Explicit factorization of external coordinates in constrained statistical mechanics models

Fixman compensating potential for general branched molecules

The technique of constraining high frequency modes of molecular motion is an effective way to increase simulation time scale and improve conformational sampling in molecular dynamics simulations. However, it has been shown that constraints on higher frequency modes such as bond lengths and bond angles stiffen the molecular model, thereby introducing systematic biases in the statistical behavior of the simulations. Fixman proposed a compensating potential to remove such biases in the thermodynamic and kinetic properties calculated from dynamics simulations. Previous implementations of the Fixman potential have been limited to only short serial chain systems. In this paper, we present a spatial operator algebra based algorithm to calculate the Fixman potential and its gradient within constrained dynamics simulations for branched topology molecules of any size. Our numerical studies on molecules of increasing complexity validate our algorithm by demonstrating recovery of the dihedral angle probability distribution function for systems that range in complexity from serial chains to protein molecules. We observe that the Fixman compensating potential recovers the free energy surface of a serial chain polymer, thus annulling the biases caused by constraining the bond lengths and bond angles. The inclusion of Fixman potential entails only a modest increase in the computational cost in these simulations. We believe that this work represents the first instance where the Fixman potential has been used for general branched systems, and establishes the viability for its use in constrained dynamics simulations of proteins and other macromolecules.

Exact and efficient calculation of Lagrange multipliers in biological polymers with constrained bond lengths and bond angles: proteins and nucleic acids as example cases

To accelerate molecular dynamics simulations, it is common to impose holonomic constraints on the hardest degrees of freedom. In this way, the time step used to integrate the equations of motion can be increased, thereby allowing longer total simulation times. The imposition of such constraints results in an aditional set of N(c) equations (the equations of constraint) and unknowns (their associated Lagrange multipliers), whose solution is closely related to any algorithm implementing the constraints in Euclidean coordinates. In this work, it is shown that, due to the essentially linear structure of typical biological polymers the algebraic equations that need to be solved involve a matrix which is not only sparse, but also banded if the constraints are indexed in a skilful way. This allows the Lagrange multipliers to be obtained through a noniterative procedure, which can be considered exact up to machine precision, and which takes O(N(c)) operations, instead of the usual O(N c3) for generic molecular systems. We develop the formalism, and describe the appropriate indexing for a number of model molecules. Finally, we provide a numerical example of the technique in a series of polyalanine peptides of different lengths. Although a use of the Lagrange multipliers without any modification in the solution of the underlying ordinary differential equations yields unstable integration algorithms, the central role of these quantities makes their efficient calculation useful for the improvement of methods that correctly enforce the exact satisfaction of the constraints at each time step. We provide several examples of this.

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 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.