Molecular Dynamics - Constraints

The RATTLE Algorithm

Constraints are applied during dynamics runs via the RATTLE algorithm. RATTLE is the velocity version of SHAKE (Ryckaert et al. 1977) which is a popular algorithm for introducing distance constraints during molecular dynamics simulations. The RATTLE procedure is to go through the constraints one by one, adjusting the coordinates so as to satisfy each in turn. The procedure may be iterated until all the constraints are satisfied to within a given tolerance. Unlike SHAKE, RATTLE (Andersen 1983) makes sure that velocities are adjusted also, to satisfy the constraints, and is suitable for use with the velocity version of Verlet integrators.

RATTLE can be used to constrain bonds, angles spanned by two constrained bonds, as well as any pair of atoms in periodic and nonperiodic systems. It can be used with the constant-volume ensembles. However, for the time being, RATTLE cannot be used with constant-pressure dynamics, since that involves calculating the pressure on a molecule basis, which is not currently available.

Bond vibrations constitute the highest frequencies in the system and thus determine the largest timestep that can be used during dynamics. If the bonds are constrained, longer timesteps can be used during dynamics, because of the absence of these high-frequency bond vibrations. In a test run on crambin, a timestep as long as 3 fs could be used with the RATTLE algorithm, but dynamics without RATTLE can have a timestep of no more than 1 fs. Furthermore, energy conservation for the 3 fs timestep with RATTLE was as good as that at 1 fs without RATTLE.

Although constraining bonds allows the timestep to be increased without incurring significant overhead for carrying out RATTLE iteratively, constraining angles does not really help to reduce the computational cost, because it takes a long time for the iterative procedure to converge. In some circumstances, it actually takes much longer than without RATTLE. It is not recommended that rattle angles be used even though it is available. If angle constraints have to be used, you should use a larger tolerance than with bonds.

Another major use of RATTLE is to enable the use of a fixed-geometry water model. With RATTLE, you can use the SPC, TIP3P, or current water models. The initial configurations of the water molecules are set to the equilibrium geometry of the SPC or TIP3P model or to the current geometry of the water model.

Tests of RATTLE Implementation

RATTLE has been tested on two systems--crambin (a nonperiodic system) and a box of 216 water molecules with periodic boundary conditions.

Crambin

NVE dynamics was performed on crambin with and without RATTLE. No cutoff was used for the nonbond interactions. When RATTLE was used, all the bonds were constrained with a tolerance of 1e-7 and then 1e-5. Dynamics was run for 100 fs with different timesteps. The goal was to check how large the timestep could be before the job would crash due to failure of the integrators.

The standard deviation (SDV) of the total energy in an NVE simulation is commonly used to indicate how accurate the integrator is. In Table 5-2, the SDV in total energy and the CPU time per step are listed as functions of the timestep.

With RATTLE, a much longer timestep may be used before the dynamics run crashes (5.5 fs vs 3 fs). In addition, the standard deviation is a lot smaller for the same timestep. For similar standard deviations (accuracies), as long as 3.5 fs can be simulated with RATTLE, compared to 1 fs without RATTLE. Comparing the RATTLE runs with tolerances of 1e-5 and 1e-7, the SDV in total energy was virtually the same; however, using a larger tolerance saved CPU time by about 4%.

Table 5-3 shows the timing of a dynamics simulation with RATTLE constraining i) only bonds lor ii) both bonds and angles. With angle constraints, the tolerance must be much looser than 1e-7. The calculation converges only when the tolerance is increased to 1e-3. The run takes much longer (more than twice as long) to converge when angle constraints are used. Hence, angle constraints are not recommended for the purpose of using a longer timestep.

Table 5-2. Relation of standard deviation in total energy and the CPU time per step with and without RATTLE

               No RATTLE                    RATTLE
               ---------        ------------------------------------
                                 tol. = 1e-7          tol. = 1e-5
                                ----------------    ----------------
   Timestep  SDV     Time/step   SDV   Time/step    SDV     Time/Step
    (fs)  (kcal mol^-1) (sec)    (kcal mol^-1) (sec)  (kcal mol^-1)  (sec)

    0.5     0.740    0.889      0.058   0.971      0.058     0.937
    1.0     2.941    0.893      0.232   0.985      0.232     0.944
    1.5     6.352    0.884      0.525   0.988      0.526     0.953
    2.0     9.832    0.886      0.939   0.993      0.939     0.957
    2.5    10.217    0.887      1.494   0.998      1.494     0.967
    3.0            crashed      2.224   1.004      2.224     0.961
    3.5                         3.239   1.008      3.239     0.966
    4.0                         4.859   1.006      4.858     0.967
    4.5                         6.821   1.011      6.820     0.970
    5.0                        11.126   1.012     11.125     0.973
    5.5                           crashed             crashed

(``tol'' = tolerance)

Table 5-3. CPU time used per step when using i) bond constraints and ii) both bond and angle constraints

 Condition                                  CPU time/step
 ---------                                  -------------

No RATTLE with timestep of 1fs               0.893 sec
RATTLE with bonds (tol 1e-5) 3fs             0.961 sec
RATTLE with bonds and angles (tol 1e-3)      2.325 sec
RATTLE with bonds and angles (tol 1e-5)      job not finished in 15 min

Water

A box of 216 water molecules under periodic boundary conditions was used for a study of the constant-temperature constant-volume dynamics (NVT ensemble) of water with a group-based cutoff (using the default value). A timestep of 1 fs was used. The CPU time required for each step with different tolerances are shown in Table 5-4.

For water, one great advantage of RATTLE is that fixed-geometry water models like SPC and TIP3P can be used. Both models are now available in the CVFF and AMBER forcefields. SPC and TIP3P cannot be used with the CFF93 forcefield, because the van der Waals term is 6-12 for these water models but 6-9 for CFF93. The SPC and TIP3P water models are also not available with ESFF. Thus, to use RATTLE with the latter two forcefields, you have to set the water type to be Current.

Table 5-4. Simulation of water with and without RATTLE

         Tolerance (Å)              Time/step (sec)
         -------------              ---------------

         1e-4                       0.842
         1e-5                       0.894
         1e-6                       0.931
         1e-7                       0.964

        no RATTLE                   0.781

The intermolecular interactions for SPC and TIP3P water were -10.24 kcal mol^-1 and -9.87 kcal mol^-1, respectively. This compares well with the experimental value of -9.92 kcal mol^-1, as well as with other calculations on the SPC and TIP3P models (-10.18 and -9.82 kcal mol^-1 in Jorgensen et al. 1983).

In conclusion, RATTLE not only allows the use of longer timesteps (3 fs with reasonable accuracy) but the use of rigid water models. For crambin, RATTLE enables dynamics to run 2.79 times faster by allowing the use of timestep of 3 fs vs 1 fs without RATTLE, with the overhead for the RATTLE algorithm taken into account.

Main access page

Theory/Methodology access.

Dynamics access Dynamics - Pressure & Stress Dynamics - General Methodology