Environmental Engineering Reference
In-Depth Information
�
−
�
(
)
(
)
(
)
(
)
χ +d =χ
t
t
rt
+d
t
,
pt
+d +m +d =
t
gt
t
0
€
‚
(131)
ƒ
„
(which is actually linear, since and substituting back
−
+d = d +m +d
(
)
(
)
(
)
pt
t
pt
gt
t
(132)
In pseudocode this scheme may be written:
do step
=
1, nstep
dt
(
)
p
=+
p
*f
2
p
r
=+
r
dt*
m
lambda_g
=
shake(r)
(133)
p
=+
p
lambda
lambda_g
g
r
=+
r
dt*
m
( )
f
=
force r
dt
(
)
p
=+
p
*f
2
mu_g
=
=+
rattle(r, p)
p
p
mu_g
enddo
The routine called shake here calculates the constraint forcesnecessary to ensure
that the end-of-step positionssatisfy Eq. (114). For a system of many constraints,
this calculation is usually performed in an iterative fashion, so as to satisfy each
constraint in turn until convergence; the original SHAKE algorithm was framed
in this way. These constraint forces are incorporated into both the end-of-step
positions and the mid-step momenta. The routine called rattle calculates a new set
of constraint forcesto ensure that the end-of-step momenta satisfy Eq. (115). This
also may be carried out iteratively. It is important to realize that a simulation of
a system with rigidly constrained bond lengths, is not equivalent to a simulation
with, for example, harmonic springs representing the bonds, even in the limit of
very strong springs. A subtle but crucial, difference lies in the distribution func-
tion for the other coordinates. If we obtain the conjugational distribution function
by integrating over the momenta, the difference arises because in one case a set of
momenta is set to zero, and not integrated, while in the other an integration is per-
formed, whichmay lead to an extra term depending on particle coordinates. This
is frequently called the metric tensor problem; it is explained in more detail in the
Search WWH ::
Custom Search