Holonomic constraints : an analytical result

0.1 MB·English

by  Martial Mazars

7
0
0
2

n
Holonomic constraints : an analytical result.
a
J

0
1
] Martial MAZARS∗
h
Laboratoire de Physique Th´eorique (UMR 8627),
c
e
Universit´e de Paris Sud XI, Bˆatiment 210, 91405 Orsay Cedex, FRANCE
m
-
t
a February 6, 2008
t
s
.
t
a
m
L.P.T.-Orsay 07-01
-
d
n
o
c Abstract
[

1
v Systemssubjectedtoholonomicconstraintsfollowquitecomplicateddynamicsthatcould
1
not be described easily with Hamiltonian or Lagrangian dynamics. The influence of
2
holonomic constraints in equations of motions is taken into account by using Lagrange
2
1 multipliers. Finding the value of the Lagrange multipliers allows to compute the forces
0
induced by the constraints and therefore, to integrate the equations of motions of the
7
0 system. Computing analytically the Lagrange multipliers for a constrained system may
/ be a difficult task that is dependingon the complexity of systems. For complex systems,
t
a
it is most of the time impossible to achieve. In computer simulations, some algorithms
m
usingiterative procedures estimate numerically Lagrange multipliers or constraint forces
-
d by correcting the unconstrained trajectory. In this work, we provide an analytical com-
n putation of the Lagrange multipliers for a set of linear holonomic constraints with an
o
arbitrary number of bonds of constant length. In the appendix of the paper, one would
c
: find explicit formulas for Lagrange multipliers for systems having 1, 2, 3, 4 and 5 bonds
v
i of constant length, linearly connected.
X
r
a
∗
Electronic mail: [email protected]
1
I. Introduction.
The characteristic time scales associated with intramolecular motions are around 10 to
100 times shorter than the charecteristic time scales associated with translational and
rotational degrees of freedom of the molecule. In Molecular Dynamic computations a
convenient way to handle the multiple-time-scale in equation of motions is to treat cova-
lent bonds between atoms as rigid, this procedureintroduces constraints on equations of
motions ; the use of the Lagrange multipliers method allows to integrate the equations
of motion [1].
From a technical point of view, the Lagrange multipliers method is used in SHAKE and
RATTLE algorithms : the constrained trajectories of the systems are computed from
unconstrained trajectories, by using an iterative algorithm that allows to estimate nu-
merically the Lagrange multipliers [2, 3, 4]. Thesealgorithms are very efficient and quite
convenient. Nevertheless, a closed analytical computation of the Lagrange multipliers
is interesting on its own on a theoretical ground and can also be used to improve the
iterative estimations done in SHAKE and RATTLE algorithms.
The equations of motions of the constrained dynamics are derived as follow. For a set of
holonomic constraints, noted as σ 0, one defines a constrained Lagrangian L [q,q˙] as
α c
≡
L [q,q˙] = L[q,q˙] λ σ (q,q˙) (1)
c α α
−
α
X
where L[q,q˙] is the Lagrangian of the unconstrained systems and λ the set of Lagrange
α
multipliers. The equations of motions are given by
∂ ∂L ∂L
′ ′
= (2)
∂t ∂q˙ ∂q
which may also be written as
∂σ
α
m q¨ = F + λ = F + G (3)
i i i α i i;α
∂q
α i α
X X
where the constraint forces are defined by G . Since holonomic constraints are con-
i;α
served quantities, by requiring that the second time derivatives of all σ vanish, we
α
obtain a set of closed equations that allows to compute the Lagrange multipliers. We
have
∂σ˙ 1
α = (F + G ) σ + q˙q˙ σ
i i;β i α i j i,j α
∂t m ∇ ∇
i i β i,j
X X X
(4)
= K +C + Z λ = 0
α α αβ β
β
X
2
where we have set
1
K = F σ
α i i α
m ·∇
 i i
X
 Cα = Xi,j q˙iq˙j∇i,jσα (5)
1
Then, Lagrange multipliers areZgαivβen=byXifinmdiin∇gitσhαe·i∇nvieσrβse of the Z-matrix, as
λ = (Z 1) (K +C ) (6)
β − βα α α
−
α
X
Obviously, the structure of the Z-matrix is closely related to the geometrical features of
molecules, its computation can be quite complicated and is depending on the system.
In this work, by computing the inverse of the Z-matrix, we provide an analytical result
for holonomic constraints involved in the dynamic of linear molecules (polymers, etc.)
or having a linear sequence in its structure (branched polymers, etc). For such systems,
the holonomic constraints can be read as
σ (r ,r ) (r r )2 a2 = 0 (7)
α α α 1 α α 1
− ≡ − − −
where r are the coordinates of the atom numbered α in the linear sequence of the
α
molecule and a the length of the bond between atoms α and α 1. N is the number of
−
holonomic constraints, these constraints involve N bonds and N +1 atoms.
In the next section, we provide a closed analytical computation of the Z 1-matrix for
−
the linear holonomic constraints given by Eq.(7).
II. Analytical computation of Lagrange multipli-
ers for linear holonomic constraints.
We setr r = auˆ , whereuˆ is theunitarybondvector linkingtheatom numbered
α α 1 α α
− −
α 1 to atom α ; the linear holonomic constraints Eq.(7) gives
−
σ = 2a(δ δ )uˆ (8)
i α i,α i,α 1 α
∇ − −
whereweusetheKronecker symbolδ . Theforceacting onatom iduetoallholonomic
n,p
constraints is G = G and we found easily that
i α i;α
P
G = 0
i
i
X
3
This may be easily verified for linear holonomic constraints.
Thus, for a set of linear holonomic constraints, the Z-matrix is given by
2a2
Z = (δ δ )(δ δ )uˆ uˆ
αβ i,α i,α 1 i,β i,β 1 α β
m − − − −
i
X (9)
4a2 1 1
= ( (uˆ uˆ )δ +δ (uˆ uˆ )δ )
α 1 α α,β+1 α,β β 1 β α,β 1
m − 2 − · − 2 − · −
In the following, we set
1
γ = uˆ uˆ and Z˜ = γ δ +δ γ δ (10)
α α 1 α αβ α α,β+1 α,β β α,β 1
2 − · − − −
In appendix, computations for N = 1,2,3,4 and 5 are given explicitly. With the linear
holonomic constraints the Z-matrix is a banded symmetric matrix. To compute the
inverse matrix we take into account the properties of Z˜. The inverse matrix is built by
a sequence of multiplications on the right and the left in such a way that the identity
matrix appears progressively on the diagonal. More precisely, according to Eq.(10), if
we perform the transform
Z˜ BtAtZ˜A B = Z˜(1) (11)
−→ 1 1 1 1
with matrix A and B defined by
1 1
( 1 )
(A ) = δ +γ δ δ and (B ) = δ + 1 δ δ (12)
1 ij i,j 2 i,1 2,j 1 ij i,j i,2 2,j
1 γ2 −
− 2
q
then, the 2 2 identity matrix appears in Z˜(1) for 1 i 2 and 1 j 2. The
× ≤ ≤ ≤ ≤
matrix multiplications by A and At result in vanishing the non-diagonal coefficient (i.e.
1 1
Z˜ = Z˜ = γ ), while matrix multiplications by B and Bt scale to 1 the diagonal
12 21 − 2 1 1
coefficient (AtZ˜A ) . The matrix Z˜(1) is similar to the original matrix Z˜, with γ and
1 1 22 2
γ transformed as
3
γ γ = 0
2 −→ 2′
γ (13)
3
γ γ =
3 −→ 3′ 1 γ2
− 2
q
because of Eq.(11).
We arethenin apositionto buildtheidentity matrix fromZ˜, by transformations similar
to Eq.(11). The matrix Z˜(n) is computed from Z˜(n 1) by making
−
Z˜(n) = BtAtZ˜(n 1)A B (14)
n n − n n
With a convenient choice of the matrix A and B , we will finally obtain Z˜(N 1) = I
n n − N
and therefore, we will have built a matrix C such as
N 1
CtZ˜C = I with C = − (A B ) (15)
N n n
n=1
Y
4
and thus
Z˜ 1 = CCt (16)
−
To this end, taking into account the structure of Eq.(13) for γ , we define the sequence
3′
Ω by
n
γ2
Ω2 = 1 n (17)
n − Ω2
n 1
−
with γ = 0, Ω = 1 and, for convenient reasons, γ = 0, that gives Ω = 1. With
1 1 N+1 N+1
these notations the matrix A and B are given by
n n
γ
n+1
(A ) = δ + δ δ
n ij i,j i,n (n+1),j
Ω
n
(18)
( 1 )
(B ) = δ + 1 δ δ
n ij i,j i,(n+1) (n+1),j
Ω −
n+1
From the definition of matrix A and B and with Eq.(16), we may easily compute
n n
determinants as
1
Det(A )= 1 ; Det(B )=
n n
Ω
n+1
and
N 1
Det(Z˜ 1) = − 1 = 1 (19)
− Ω2 Det(Z˜)
n=1 n+1
Y
Det(Z˜) is closely related to the metric determinant that is used in computations of
ensemble averages [5, 6, 7].
Eqs.(16) and (17) allow to compute the Z˜ 1 matrix coefficients, the matrix is symmetric
−
Z˜ 1) = Z˜ 1) , and after some algebra we found
− ij − ji
N p
For i = j : Z˜ 1) = 1 [1+ ( 1 1)]
 − ii Ω2i p=i+1m=i+1 Ω2m −
 For i < j : Z˜ 1) = 1 [ j Xγn][Y1+ N p ( 1 1)] (20)
WithEqs.(10) and (20), w−e mijay verΩif2iy tnh=Yai+t1Z˜Ω2n1Z˜ = Ip=X.jE+1qms.=(Y1j9+)1anΩd2m(2−0) are the main
− N
results of this paper.
According Eq.(6), we have to compute K and C to achieve the analytical computation
α α
of the Lagrange multipliers. With Eqs.(5) and (8), we found
2a
K = (F F ) uˆ (21)
α α α 1 α
m − − ·
5
and
2
C = (p p )2 (22)
α m2 α− α−1
with the momentum of atoms, p = mv = mr˙ and notations defined previously.
α α α
Lagrange multipliers are then given by applying Eq.(6).
So far, we have considered that all atoms, involved in the set of holonomic constraints,
have the same mass m and the length of each bond is equal to a. We consider now that
the sequence of atoms mass is m and the length of bonds are a . Then,
i i [0,N] i i [1,N]
the momentum of atoms is now{ p}=∈ m v = m r˙ and Eq.(8) becom{es}∈
α α α α α
σ = 2a (δ δ )uˆ (23)
i α α i,α i,α 1 α
∇ − −
Therefore, the Z-matrix is now given by
Z = ν (uˆ uˆ )δ +ω δ ν (uˆ uˆ )δ
αβ α α 1 α α,β+1 α α,β β β 1 β α,β 1
− − · − − · −

 να = −2amαα−−1a1α (24)
Eq.(24) is simωilaαr to=E2qa.(2α9().mWα1−e1m+amy1αtr)ansform Eq.(24) to Eq.(10) with the help of a
diagonal matrix D, according to
Z˜ = DZD (25)
with
1
D = δ (26)
ij ij
√ω
i
Then, Z˜ is again given by Eq.(10), but with
m m
α α 2
γα = − (uˆα 1 uˆα) (27)
s(mα 1+mα)(mα 1 +mα 2) − ·
− − −
where inertia parameters of bonds appear explicitly (the canonical partition function
of freely jointed chains is depending on inertia parameters of the bonds as defined in
ref.[8]). One may note also that if m = m, Eq.(27) corresponds to Eq.(10).
α
The inverse of the Z-matrix is now given by
Z 1 = DZ˜ 1D = D(CCt)D (28)
− −
with Eq.(27), (17) and (20) ; and K and C are given by
α α
F F
Kα = 2aα( α α−1) uˆα (29)
m − m ·
α α 1
−
and
Cα = 2(pα pα−1 )2 (30)
m − m
α α 1
−
6
III. Discussion.
From an analytical point of view, an interesting result is given by the computation of
the determinant of Z˜ (Eq.(19)). For holonomic constraints, this determinant is known
as the metric determinant ; and for non-hamiltonian dynamical systems [5, 6, 9] it can
be related to the Jacobian determinant, some ensemble averages can be formulated with
the help of this determinant [5, 6, 10]. According to computations done in ref.[11], we
obtain
1 N 1 N
2N 1 2 dγ Det(Z˜ 1) = 2N 1 2 dγ 1
− n − − n
Z0 n=2 q Z0 n=2 Det(Z˜)
Y Y (31)
q
1 3 1 1
= [N 1]He[N]( ; ; ,..., )
−
2 2 4 4
where[N 1]He[N] is themultiplehypergeometric function definedinref.[8], related to the
−
canonical partition function of freely jointed chains. In particular, for N = 2, Eq.(31) is
equivalent to the elementary relation
1
2 dγ2 1 1 3 1 1
2 = F ( , ; ; )= 2arcsin( ) (32)
2 1
Z0 1−γ22 2 2 2 4 2
q
In ref.[8], we have shown that themultiple hypergeometric function [N 1]He[N] is defined
−
for any mass sequence, Eq.(31) may therefore be extended to any mass sequence, and
the multiple hypergeometric function [N 1]He[N] must be evaluated at a point defined
−
by the inertia parameters of the coupling between bonds [8]. The inertia parameter of
the coupling between bonds α and α 1 is given by
−
m m
α α 2
xα = − (33)
(m +m )(m +m )
α 1 α α 1 α 2
− − −
this parameter appears explicitly in Eq.(27), so in the integral too (one may note also
that if m = m, for all α, then we have x = 1/4). For any mass sequence, the relation
α α
equivalent to Eq.(31) is
√x2dγ ... √xN dγ Det(Z˜ 1) = ( N x1/2) [N 1]He[N](1;3; x ) (34)
Z0 2 Z0 Nq − n=2 α − 2 2 { i}i∈[2,N]
Y
with γ given by Eq.(27). Eq.(34) can be considered as an integral definition of the
α
multiple hypergeometric function [N 1]He[N] ; this multiple hypergeometric function is
−
quite complicated, an iterative approximation scheme, called Independent Motions Ap-
proximation (IMA), is described in ref.[8].
From the point of view of the physical chemistry, the main purpose of this paper was to
computeanalytically theLagrangemultipliers for aset oflinear holonomic constraints in
view of applications to molecular simulations or statistical physics of complex systems.
7
For applications oftheseanalytical resultsinmolecular simulation, themaindrawback is
that the number of coefficients in the Z˜ 1-matrix, that have to be computed and stored,
−
scales as N2 (the matrix is fully filled, see appendix for some examples), N being the
number of bonds linearly connected in the molecule. Therefore, for practical implemen-
tations, this brute force method of computing Lagrange multipliers is certainly not as
efficient as SHAKE or RATTLE algorithms, at least when N is rather large. Another
point that one should kept in mind in trying to implement these analytical results in a
computer code is that, because of C , velocities (or momentum) contribute to Lagrange
α
multipliers. This may rise to some technical problems in the Verlet-velocity algorithm
and in constant temperature Molecular Dynamics computations, since forces at t+δt
have to be computed to obtain the velocity at t+δt [1, 2, 4]. Nevertheless, clever uses
of these analytical results may certainly be used to improve the efficiency and accuracy
of RATTLE algorithms ; for instance, when many bonds verify uˆ uˆ 1.
α 1 α
On general grounds, despite the complexity of analytical formu−las·of Z≃˜ 1 given by
−
Eq.(20), we believe that the analytical results presented in this work are of some in-
terest, since matrix similar to Z˜ appears frequently in many models or theories in one
dimension or with nearest neighbors interactions (spin models, harmonic chains, etc.).
8
Appendix : Lagrange multipliers for N = 1,2,3,4 and 5.
In this appendix, we give explicitly the analytical formulas for linear holonomic con-
straints that correspond to N=1, 2, 3, 4 and 5 bonds. To obtain analytical results for
the matrix Z˜ 1, we compute firstly Ω2 from the recursive relation in Eq.(17),
− n
Ω2 = 1 γ2
2 − 2

 ΩΩ2324 == 11−−1γγ−2222γ1−−22γγ3322γ2−γ42γ(21−γ22) (A.1)
− 2 − 3
The computationofΩZ25˜ 1=fo1r−Nγ=22−1,γ1232,−3−,γ4γ2242a−(n1γd−325−γi22sγ)d42−(o1nγe−52(bγ1y22−)usγi22ng−Eγq32.)(20).
−
N = 1
•
For N = 1, Z˜ = Z˜ 1 = (1) and, following Eqs.(20) and (21), we have
−
2a 2
K = (F F ) uˆ and C = (p p )2 (A.1.1)
1 m 1− 0 · 1 1 m2 1− 0
Then, the Lagrange multiplier is
1 1
λ = [(F F ) uˆ + (p p )2] (A.1.2)
1 2a 1− 0 · 1 ma 1− 0
the constraint force acting on atom numbered 1 is given by
1
G = λ ∇ σ = [(F F ) uˆ + (p p )2]uˆ (A.1.3)
1 − 1 1 1 − 1− 0 · 1 ma 1− 0 1
and the force acting on atom 0 is G = G .
0 1
−
N = 2
•
For N = 2, we have
1 γ 1 1 γ
Z˜ = −γ2 −12 ! and Z˜−1 = 1−γ22 γ2 12 ! (A.2.1)
9
therefore, with γ = uˆ uˆ /2, Lagrange multipliers are given by
2 1 2
·
λ = 1 [(F F ) uˆ + 1 (p p )2
 1 2a(1 (uˆ1·uˆ2)2) 1− 0 · 1 ma 1− 0
− 4
 λ = 1 [+121((uˆuˆ1·uˆuˆ2)[)([F(F2−FF1))·uˆuˆ2++ m11a((pp2−pp1))22]]] (A.2.2)
The for2ces ac2tai(n1g−on(uˆa1t·4ouˆm2)s2)are+2(F21−· F21)·uˆ12−+ m01a·(p21−pm1)a2] 1− 0
G = 2aλ uˆ
0 1 1
−

 G1 = 2aλ1uˆ1−2aλ2uˆ2 (A.2.3)
G = 2aλ uˆ
2 2 2
N = 3 
•
With N = 3, matrix Z˜ and Z˜ 1 are given by
−
1 γ 0
2
−
Z˜ = γ 1 γ (A.3.1)
 2 3 
− −
0 γ 1
3
 − 
 
and
1 γ2 γ γ γ
− 3 2 2 3
1  
Z˜ 1 = γ 1 γ (A.3.2)
− 1−γ22−γ32  2 3 
 γ γ γ 1 γ2 
 2 3 3 − 2 
 
As for N = 1, following Eqs. (20) and (21), for 1 i 3, K and C are given by
i i
≤ ≤
2a 2
K = (F F ) uˆ and C = (p p )2 (A.3.3)
i m i− i−1 · i i m2 i− i−1
10