Computational Physics

238 Molecular dynamics simulations

and the remaining ones are fixed by linear constraints of the form(8.119). In the case
of a planar structure we take three noncollinear atoms as a backbone. These atoms
satisfy three bond constraints and the remaining atoms are constrained linearly. In
a three-dimensional molecular structure, four backbone atoms are subject to six
bond constraints and the remaining ones to a linear vector constraint each. In the
constraint procedure, the degrees of freedom of the nonbackbone atoms are elim-
inated so that the forces they feel are transferred to the backbone. This elimination
is always possible for linear constraints such as those obeyed by the nonbackbone
atoms.
Let us now return to our CS 2 example. First we write down the equations of
motion for all three atoms, following from the extended Lagrangian (the Lagrange
parameter for the bond constraint is calledλ, that of the linear vector constraintμμμ):

mS ̈rS( 1 )=F 1 − 2 λ(rS( 1 )−rS( 2 ))−μμμ/ 2 (8.120a)

mS ̈rS( 2 )=F 2 + 2 λ(rS( 1 )−rS( 2 ))−μμμ/ 2 (8.120b)

mCr ̈C=FC+μμμ. (8.120c)

The linear constraint(8.119)is now differentiated twice with respect to time, and
using the equations of motion we obtain

FC+μμμ=

mC

2 mS

(F 1 +F 2 −μμμ). (8.121)

We can thus eliminateμμμin the equations of motion for the S-atoms and obtain,
withM= 2 mS+mC:

mSr ̈S( 1 )=

(

1 −

mC

2 M

)

F 1 −

mC

2 M

F 2 +

mS

M

FC− 2 λ(rS( 1 )−rS( 2 )); (8.122a)

mSr ̈S( 2 )=

(

1 −

mC

2 M

)

F 2 −

mC

2 M

F 1 +

mS

M

FC+ 2 λ(rS( 1 )−rS( 2 )). (8.122b)

These equations define the algorithm for the positions of the S-atoms, and the
position of the C-atom is fixed at any time by the linear constraint.
Note that we still have an unknown parameterλpresent in the resulting equations:
this parameter is fixed by demanding that the bond constraint must hold forrS( 1 )
andrS( 2 )at all times (note that we have not yet used this constraint). It is not easy
to eliminateλfrom the equations of motion as we have done withμμμ, as the bond
length constraint is quadratic. Instead, we solve forλat each time step using the
constraint equation. We outline this procedure for our example. Suppose we have
the positionsrS( 1 )andrS( 2 )at timestandt−hand that for both these time instances
the bond constraint is satisfied. According to the equations of motion (8.122) in the