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

Abstract

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 constrained 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 flexible 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 a toy system and in three molecules of biological relevance: methanol, N-methyl-acetamide and a tri-glycine peptide.

Figure 1. Definition of the frame of reference fixed in the system.

Figure 2. Rotation associated to a change in a bond angle.

Definition of the bond angle , 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.

Figure 3. Rotation associated to a change in a dihedral angle.

Definition of the dihedral angle , associated to atom . The positive sense of rotation is indicated in the figure, and we can distinguish between two situations regarding covalent connectivity: a) principal dihedral angle, and b) phase dihedral angle (see ref. [37]).

Figure 4. Special cases.

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

Figure 5. Molecules used in the numerical calculations in this section.

(a) Methanol, (b) N-methyl-acetamide (abbreviated NMA), and (c) 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.

Figure 6. Derivatives of some selected coordinates of methanol.

Derivatives of (a) the coordinate of atom 5 in methanol, (b) the bond length associated to it, (c) the bond angle , and (d) the dihedral angle as a function of the unconstrained coordinate . Both the results of the new algorithm and those obtained by finite differences (FD) are depicted. The key for the different types of line is the same in the four graphs.

Figure 7. Metastability of the local minima in GLY4.

(a) 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) Minimum-energy value of the constrained dihedral angle in the conformation 1044 of GLY3 for different values of the displacement in the unconstrained coordinate .

Figure 8. Dependence of the error as a function of .

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

Figure 9. Potential energy of the toy system in eq. (36).

The range of and corresponds to the one explored in this work. Contour level lines and colour level indication in the surface have been added for visual comfort. All units are arbitrary.