An exact expression to calculate the derivatives of position-dependent observables in molecular simulations with flexible constraints

In this work, we introduce an algorithm to compute the derivatives of physical observables along the constrained subspace when flexible constraints are imposed on the system (i.e., constraints in which the hard coordinates are fixed to configuration-dependent values). The presented scheme is exact, it does not contain any tunable parameter, and it only requires the calculation and inversion of a sub-block of the Hessian matrix of second derivatives of the function through which the constraints are defined. We also present a practical application to the case in which the sought observables are the Euclidean coordinates of complex molecular systems, and the function whose minimization defines the constraints is the potential energy. Finally, and in order to validate the method, which, as far as we are aware, is the first of its kind in the literature, we compare it to the natural and straightforward finite-differences approach in three molecules of biological relevance: methanol, N-methyl-acetamide and a tri-glycine peptideComment: 13 pages, 8 figures, published versio

Special cases.

<p>Special cases of atoms that do not belong to the chain connecting to atom 1, but that are nevertheless used to position .</p

Rotation associated to a change in a bond angle.

<p>Definition of the <i>bond angle </i>, associated to atom , and the unitary vector corresponding to the direction around which all atoms with chains containing rotate if is varied while the rest of internal coordinates are kept constant.</p

Dependence of the error as a function of .

<p>Average normalized error in the derivatives by finite differences as a function of (see the text for a more precise definition). (<b>a</b>) Error averaged to all conformations and all atoms of the three molecular systems studied. (<b>b</b>) Error averaged to all conformations of the -coordinate of three particular -row atoms in NMA.</p

Metastability of the local minima in GLY4.

<p>(<b>a</b>) Derivative of the constrained dihedral angle , describing a peptide bond rotation in GLY3, with respect to the unconstrained coordinate for a selected set of conformations in the working set. (<b>b</b>) Minimum-energy value of the constrained dihedral angle in the conformation 1044 of GLY3 for different values of the displacement in the unconstrained coordinate .</p

Definition of the frame of reference fixed in the system.

<p>Definition of the frame of reference fixed in the system.</p

Molecules used in the numerical calculations in this section.

<p>(<b>a</b>) Methanol, (<b>b</b>) N-methyl-acetamide (abbreviated NMA), and (<b>c</b>) the tripeptide N-acetyl-glycyl-glycyl-glycyl-amide (abbreviated GLY3). Hydrogens are conventionally white, carbons are grey, nitrogens blue and oxygens red. The unconstrained dihedral angles that span the corresponding spaces are indicated with light-blue arrows, and some internal coordinates and some atoms that appear in the discussion are specifically labeled. The constrained dihedral angle is indicated by a red arrow in GLY3.</p

Rotation associated to a change in a dihedral angle.

<p>Definition of the <i>dihedral angle </i>, associated to atom . The positive sense of rotation is indicated in the figure, and we can distinguish between two situations regarding covalent connectivity: <b>a</b>) <i>principal dihedral angle</i>, and <b>b</b>) <i>phase dihedral angle</i> (see ref. <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0024563#pone.0024563-Echenique6" target="_blank">[37]</a>).</p