Constrained Simulation

In the method of constraints, a force of the form AV is applied at each step such that remains constant throughout the simulation. To see how A can be related to dA/d , we recall that the free energy A is the effective or average potential acting on . From physical intuition, it should be true that -dA/d is the average of the force acting on . In a constraint simulation, this force is equal to -A. Therefore we can expect to have d/l/d - - (A). We now make this statement more rigorous. 

In constrained simulations, the Hamiltonian M is supplemented by a Lagrange multiplier 

The additional term A( — (x)) is needed to enforce the constraint (x) = constant. This corresponds to an additional force equal to —V(A( — (x))). A is chosen such that (x) = . Therefore the force can be expressed more simply as AV . Then, the interpretation is quite natural. In order to enforce the constraint we are applying a force parallel to V which opposes the mechanical force acting on . 

The condition (x) = constant implies in particular that = 0. With this condition it is possible to derive an expression for A as a function of positions and velocities  

See (4.10) for the definition of Z/=. Therefore A is in general a function of x and x. With this expression for A, we see that if we start a simulation with = (x) and = 0 then = (x) at all times. Note that different methods, such as SHAKE or rattle [40], are employed numerically to avoid a drift of . Those methods are not based on a direct application of (4.20). 

---

The free energy differences obtained from our constrained simulations refer to strictly specified states, defined by single points in the 14-dimensional dihedral space. Standard concepts of a molecular conformation include some region, or volume in that space, explored by thermal fluctuations around a transient equilibrium structure. To obtain the free energy differences between conformers of the unconstrained peptide, a correction for the thermodynamic state is needed. The volume of explored conformational space may be estimated from the covariance matrix of the coordinates of interest, = ((Ci [13, lOj. For each of the four selected conform-... 

Fig. 4.2. Free energy computation using constraint forces. It may be difficult to sample the surface (x) = using a constrained simulation because of the presence of energy barriers separating different reaction pathways. Left a barrier is shown in the middle of the pathway from reactant A to product B. Right two barriers are shown at B...

Importantly, in contrast to constrained simulations, the system is allowed to evolve freely and in particular to explore the various pathways connecting A and B, see Fig. 4.2. This is one of the reasons why ABF can converge much faster than the method of constraints. 

Before we derive the appropriate expressions to calculate cL4/d from constrained simulations, we note an important difference between sampling in constrained and unconstrained simulations. There are two ways to gather statistics at (x) = . In unconstrained simulations, the positions are sampled according to exp —iiU while the momenta are sampled according to exp —j3K. If a constraint force is applied to keep fixed the positions are sampled according to A( (x) — x) exp —iiU. The momenta, however, are sampled according to a more complex statistical ensemble. Recall that... 

In constrained simulations, = 0 so that is not an independent variable but rather a function of q and pg. Let us discuss the implications of this fact. Consider an arbitrary function /(x) and the following average ... 

This means that the average of / (f) can be computed using a constrained simulation for a modified Hamiltonian... 

We are now going to establish a connection between the free energy A( ) and constrained simulations. Specifically, there is a very simple relation between A(f) and... 

This equation is a very important result because it shows that constrained simulations can be used to calculate dA/d . A possible algorithm would consist in running a constrained simulation with the Hamiltonian [which contains the extra Fixman potential 1/(2/ ) In Zc(q)], calculate the rate of change of -Z/J with at each step. Finally by averaging this rate of change the derivative of A can be computed. 

Compared with (4.15), we see that this new expression for constrained simulations (4.28) is somewhat similar but a striking difference is that A is a function of the velocity x whereas (4.15) involves only an average in configurational space. Those two results are linked and we show in Appendix B that (4.28) is actually a special case of (4.15). 

This expression involves only the computation of the gradient of . A comparison with (4.28) which requires the same constrained simulation but a different equation for dA/d underscores the simplicity of this new expression. A similar expression was derived by Schlitter et al. 141 -44]. 

We now describe a different approach which is simpler than the method of constraints and also very efficient. It does not require running a constrained simulation and can be performed entirely with a single molecular dynamics run. 

One reason for the inefficiencies of constraint methods is that they may prevent an efficient sampling of the set (x) = . This is illustrated by Fig. 4.2. It is common that many pathways separated by high energy barriers exist to go from A to B. In constrained simulation, the system can get trapped in one of the pathways. In the most serious cases, this leads to quasi-nonergodic effect where only a part of the set (x) is effectively explored. In less serious cases, the convergence is quite slow. 

Fig. 4.8. Free energy computation using ABF and a constrained simulation. Reprinted with permission from Darve et al. 2001 [28], Copyright 2001, American Institute of Physics...

Equations 8-34 and 8-35 are more suited for calculations based on simulation data as require a constrained simulation, considering only the vibrational ground state condition, where the energy shift does not include the short-range term which is likely to be not properly sampled in a constrained simulation with a fixed subset of classical coordinates. 

S. M. Morrill, R.G. Lane, LI. Rosen, Constrained simulated annealing for optimized radiation therapy treatment planning. Computer Methods and Programs in Biomedicine, 33 (1990) 135. 

Constrained simulations rely on the calculation of the first term of (23) from a series of simulations conducted at different values of This average force is then corrected by adding the second term of (23), and then numerically integrated to give a potential of mean force for the desired range of... 

In the other constrained simulations (constraint respectively at 2.75, 3.0, 3.5, 3.75, 3.9 A) the formation of this kind of defects has been observed in the whole simulation times, and defects transform each other in a very short time ( 100 fs). However for such values of the constraint, the reacting O atom was still bound to the framework. The reactive event occurred only when the constraint was set to an NO distance of 4.0 A. After few fs, the oxygen previously trapped in the framework defects, left the four ring region and diffused in the adjacent cage colliding with the second NOj. Such collision first led to the transient species [N02---02] , that appeared in the unconstrained nitrite sodalite -I- O2 simulation. Then, such species reached rapid equilibrium with the separated NOj and O2 compounds. 

Several molecular dynamics simulation approaches are available, including regular simulation, constrained simulation, and simulated annealing techniques. In addition, the simulated annealing dynamics simulation also was popular in computational chemistry studies of proteins and drug molecules (34.44.45). A common experimental protocol has been widely used in molecular dynamics simulation as following (34.46.47) ... 

Kleijnen, J. Response surfece methodology for constrained simulation optimization An overview. Simulation Modelling Practice and Theory, 16(l) 50-64, 2008. 

In constrained simulations near the transition state, it is difficult to sample from the delta function constraint (3.41), especially if a Monte Carlo sampling is used. It is therefore convenient to approximate the constrained weight function (3.41) by using the Gaussian approximation of the delta function in Eq. (3.15). This Gaussian constraint corresponds to a harmonic constraint potential Vconstr,y(r) [30,32] ... 

In simulations, F z ) has been interpreted as the force that needs to be applied to keep the solute at the fixed position 2 along 2. The calculation of F z ) is carried out accordingly after each molecular dynamics step the 2 coordinate of the center of mass of the solute is reset to its constrained value of 2. F (2 ) is proportional to the distance by which the solute is moved to reset its position. Averaging over the constrained molecular dynamics trajectory produces Fz z ). If constrained simulations are repeated for several different locations of the solute along 2, Afiej c z) can be calculated by integrating equation (6). 

A third area where continuum solvent models are useful is in highly constrained simulations. These include X-ray crystallographic and 2D-NMR structural refinements. In these situations, with the additional restraints (and additional computation) arising from the experimental data, the extra expense of explicit solvent models would be inappropriate. 

The structure of chymotrypsin inhibitor 2 (CI2) is particularly well-characterized, but the crystal structure and NMR solution structure have small but distinct differences. Early NMR results suggested that there is an additional pair of anti-parallel -strands in comparison with the crystal structure. Subsequent experiment and constrained simulation concluded that the difference between these experimental techniques is caused by the refinement methodology used. 

---