In general, constrained molecular dynamics generates biased statistical averages.
It can be shown that the correct average for a quantity can be obtained using the formula:
where stands for the statistical average of the quantity enclosed in angular parentheses computed for a constrained ensemble and is a mass metric tensor defined as:
It can be shown that the free energy gradient can be computed using the equation:
where is the Lagrange multiplier associated with the parameter used in the SHAKE algorithm.
The free-energy difference between states (1) and (2) can be computed by integrating the free-energy gradients over a connecting path:
Note that as the free-energy is a state quantity, the choice of path connecting (1) with (2) is irrelevant.
Constrained molecular dynamics is performed using the SHAKE algorithm..
In this algorithm, the Lagrangian for the system is extended as follows:
where the summation is over r geometric constraints, is the Lagrangian for the extended system, and λi is a Lagrange multiplier associated with a geometric constraint σi:
with ξi(q) being a geometric parameter and ξi is the value of ξi(q) fixed during the simulation.
In the SHAKE algorithm, the Lagrange multipliers λi are determined in the iterative procedure:
- Perform a standard MD step (leap-frog algorithm):
- Use the new positions q(t+Δt) to compute Lagrange multipliers for all constraints:
- Update the velocities and positions by adding a contribution due to restoring forces (proportional to λk):
- repeat steps 2-4 until either |σi(q)| are smaller than a predefined tolerance (determined by SHAKETOL), or the number of iterations exceeds SHAKEMAXITER.
- For a constrained molecular dynamics run with Andersen thermostat, one has to:
- Set the standard MD-related tags: IBRION=0, TEBEG, POTIM, and NSW.
- Set MDALGO=1, and choose an appropriate setting for ANDERSEN_PROB.
- Define geometric constraints in the ICONST-file, and set the STATUS parameter for the constrained coordinates to 0.
- When the free-energy gradient is to be computed, set LBLUEOUT=.TRUE.
- ↑ E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989).
- ↑ W. K. Den Otter and W. J. Briels, Mol. Phys. 98, 773 (2000).
- ↑ E. Darve, M. A. Wilson, and A. Pohorille, Mol. Simul. 28, 113 (2002).
- ↑ P. Fleurat-Lessard and T. Ziegler, J. Chem. Phys. 123, 084101 (2005).
- ↑ a b J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comp. Phys. 23, 327 (1977).