Rigid bonds

Rigid water is an option (rigid bond angle and bond lengths). Also bonds between hydrogens and heavy atoms can be constrained to their nominal length during integration. Three options for selecting what bonds to constrain are given in NAMD. The configuration parameter rigidBonds can be set to none, which indicates no constrained lengths, water to constrain only bonds in water molecules, or all to constrain all bonds between hydrogens and heavy atoms. Bonds between heavy atoms cannot currently be constrained, and NAMD will not work correctly (even with rigid bonds inactive) if the system contains hydrogen molecules, since no heavy (or parent) atom can be found for either hydrogen. The rigid bond code also cannot handle methane molecules, or any other molecule with four or more hydrogens bonded to a heavy atom.

The rigid bonds module sorts the atoms of interest into one of four data structures, according to how many hydrogen atoms are bonded to the parent atom. After applying the normal integration procedure, the atom positions are adjusted to bring the bond lengths back to their neutral length (ie. the length such that the bond exerts no force on either atom). These adjustments are applied in such a way that the center of mass of the bonded atoms does not move. For the one-hydrogen case, this is a straightforward calculation, but for the two- or three-hydrogen cases, an iterative procedure is followed, where each hydrogen-parent atom pair is adjusted until all the hydrogen bonds are within a certain tolerance of the desired length, or a maximum number of iterations is exceeded. These values are set by the rigidTolerance and rigidIterations parameters. A modified algorithm is applied to waters. Although waters have only two hydrogens, a constraint is also applied to the distance between the two hydrogens and the three-bond procedure followed, so that the angle between the hydrogens is also constrained.

After correcting the atom positions, the calculated velocities must also be adjusted to account for the changed positions. Without velocity correction, the kinetic energy is computed incorrectly. Velocity correction is performed by the Rattle module. The velocities calculated by the normal integration are adjusted by removing the component parallel to the bond, without changing the velocity of the center of mass of the bonded pair.

• Constraints: $\chi^{\alpha} (X)= 0, \alpha = 1, 2,\dots,$whatever
• Algorithm: given $X_{n}, V_{n-\frac{1}{2}}$
  • crunch numbers: F_n = forces
  • half-kick: $\bar{V}_{n} = V_{n-\frac{1}{2}} + \frac{1}{2} \Delta t M^{-1} F_{n}$
  • half-kick: $V_{n+\frac{1}{2}} = \bar{V}_{n} + \frac{1}{2} \Delta t M^{-1} F_{n}$
  • drift: $\bar{X}_{n+1} = X_{n} + \Delta t V_{n+\frac{1}{2}}$
  • SHAKE: $X_{n+1} = \bar{X}_{n+1} + \frac{1}{2} \Delta t^2 M^{-1} \sum_{\beta} \lambda_{\beta} \chi^{\beta}_x (X_{n})$
    • where $\{ \lambda_{\beta} \}$ solve
    • $\{\chi_{\alpha} (X_{n+1}) = 0 \}$
    • $\{V_{n+\frac{1}{2}} = (X_{n+1}-X_{n})/\Delta t\}$
  • RATTLE: Compute the new velocities using: $V_n = (I - M^{-1} \chi_x^{\rm t} ( \chi_x M^{-1} \chi_x^{\rm t} )^{-1} \chi_x) \bar{V}_{n}$

Note that the optional RATTLE step is for the purpose of computing kinetic energy. It forces velocities to satisfy velocity constraints, but these adjustments are exactly undone by the next SHAKE.

Rigid bond calculation can be performed with Verlet integration, with or without minimization or velocity rescaling. Although NAMD will run with Langevin integration and rigid bonds both active, the equilibrium temperature of the system does not converge to the correct value. Different Langevin equations are needed to produce the correct velocity distribution, which have not yet been added to NAMD.