Computational Physics

8.6 Molecular systems 239

Verlet scheme, predictions for the positions att+hare given by

rS( 1 )(t+h)= 2 rS( 1 )(t)−rS( 1 )(t−h)+h^2

(

1 −

mC

M

)

F 1 (t)

−h^2

mC

M

F 2 +h^2

mS

M

FC(t)− 2 h^2 λ(t)[rS( 1 )(t)−rS( 2 )(t)]; (8.123a)

rS( 2 )(t+h)= 2 rS( 2 )(t)−rS( 2 )(t−h)+h^2

(

1 −

mC

M

)

F 2 (t)

−h^2

mC

M

F 1 +h^2

mS

M

FC(t)+ 2 h^2 λ(t)[rS( 1 )(t)−rS( 2 )(t)]. (8.123b)

The predictions for the positions att+hare linear functions ofλand if we substitute
them into the bond constraint (8.118), we obtain a quadratic equation forλwhich
can be solved trivially. Of the two solutions, we keep the smallest value ofλ. This
means that the bond constraint is now satisfied to computer precision for all times.
It should be noted that the value ofλin this procedure will deviate slightly from its
value in the exact solution of the continuum case, but the deviation remains within
the overall orderh^4 error of the integration scheme[52].
We have given the CS 2 example here because it illustrates the general procedure
involving linear constraints which are all eliminated from the equations of motion,
thereby reducing the degrees of freedom to those of the backbone atoms (the two
sulphur atoms in our example). These are confined by quadratic bond constraints.
The Lagrange multipliers associated with the latter are kept in the problem and
fixed by the bond constraints themselves. After solving for the backbone, the linear
constraints fix the positions of the remaining atoms uniquely.
The nitrogen molecule which was discussed in the previous subsection using a
direct approach can be treated with the method of constraints. It is a simple problem
because there are no linear constraints which have to be used to remove redundant
degrees of freedom from the equations of motion, and we are left with the following
equations:

m 1 r ̈ 1 =F 1 +λ(r 1 −r 2 ) (8.124a)

m 2 r ̈ 2 =F 2 −λ(r 1 −r 2 ). (8.124b)

The Verlet equations lead again to linear predictions forr 1 andr 2 at the next time
step and substituting these into the bond constraint leads to a quadratic equation
which fixes the Lagrange multiplierλ. For an implementation, see Problem 8.9.

8.6.3 General procedure: partial constraints

In the previous section we have considered systems consisting of completely rigid
molecules. Now we discuss partially rigid molecules, consisting of rigid fragments
which can move with respect to each other. The motion of two fragments attached