Journal of Chemical Physics, Vol.100, No.12, 9129-9139, 1994
Change of Bond-Length in Free-Energy Simulations - Algorithmic Improvements, But When Is It Necessary
Various methods are compared for calculating the constraint free energy in simulations with change of bond lengths. Bond forces can be evaluated as restraint forces with the use of flexible bonds, and also, when bend lengths are held fixed, as constraint forces calculated either from SHAKE displacements, or else from Gauss’ principle of least constraint, and these forces can then be used to calculate the potential of mean force or free energy for the bond length changes. An alternative is to use so-called mean-potential forces. Analysis of the equations derived from Gauss’ principle shows that the constraint force (and also the restraint force) contain a "dynamic" term, that depends on the molecular motion and for simple molecules averages to 2 kT/r, where r is the bond length. Integration of this term, which is not present with use of the mean-potential force, or in a Monte Carlo simulation, gives a free energy that corresponds to the change in conformational entropy of two particles maintained at a certain distance, when this distance is changed. It is unnecessary to compute this free energy term by simulation, as was done in earlier work. The accuracy and precision of these methods, and of methods in which the bond lengths do not change, are tested in a series of model calculations of increasing complexity. First, the equivalence of these methods is tested in a model system with small molecules whose motion is maintained by Brownian dynamics, and which can be subjected to very simple external forces, in order to permit very long simulations. Second, a comparison is made of techniques for calculating free energy change with and without change of bond lengths, applied to the transformation of ethane into ethanol in water. A third test system is the conversion of an alanine residue into a serine residue in an ct-helix in water. It is found that results are much improved if the "vanishing" atoms (in this case hydroxyl O and H) are coupled to a heat bath via Brownian forces. Statistical errors are estimated from four successive slow-growth simulation cycles. Calculations with twofold change of C-O and O-H bond lengths and linear coupling of the nonbonded energy terms are more precise than calculations without change of bond lengths and nonlinear coupling. The free energy change for the system with flexible bonds contains a significant contribution (circa 3 kJ mol(-1)) from elongation of the bond lengths as a result of interaction of the polar OH group with water. The constraint forces calculated from the SHAKE displacements and from Gauss’ principle turn out to be almost identical. The three methods for calculating the constraint forces achieve similar precision in the ethane/ethanol system as well as in the alpha-helix system. We also discuss the application of flee energy calculations with changing bond lengths to changes of a molecule’s chemistry in the middle of a chain.