- Email: [email protected]
A new approach to constrained molecular dynamics
Journal of Molecular Structure 563±564 (2001) 135±139
www.elsevier.nl/locate/molstruc
A new approach to constrained molecular dynamics Natale Neto* Department of Chemistry, University of Firenze, Via G. Capponi 9, 50121 Firenze, Italy Received 31 August 2000; revised 31 January 2001; accepted 31 January 2001
Abstract The internal motion of a semirigid molecule is described in terms of variations of intramolecular distances, valence angles, out of plane wagging and torsions. Constraints may be imposed to any selected group of these parameters and the reduction of the number of degrees of freedom is realized through a suitable basis change, splitting the space of displacements into two orthogonal subspaces, one of them spanned by the coordinates forced to vanish. The treatment differs from the previous ones, which all rely on the Lagrangian l method and thus dispense with the de®nition of unconstrained coordinates, which do not eventually appear also in the expression of the constrained Cartesian displacements given here. It turns out that the previous shake algorithm of the constrained dynamics include unnecessary calculation at each time step. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Constraints; Molecular dynamics; Vibrational dynamics
1. Introduction There is undoubtedly a close similarity between vibrational and molecular dynamics (MD) as both techniques largely use methods of classical mechanics and deal with small, but not in®nitesimal, displacements from a reference geometry. In fact in MD the motion is partitioned into a sequence of steps, each of them occurring in a very tiny time interval during which all atoms move from the initial positions. The corresponding conformational change of a polyatomic molecule may be then be described in terms of internal coordinates of the kind de®ned in the book of Wilson et al. [1], namely stretching, bending, wagging out of plane, torsion or combinations of them. These coordinates are routinely used in vibrational spectroscopy but not in constrained MD. In fact * Tel.: 139-552-757597; fax: 139-552-476961. E-mail address: [email protected]®.it (N. Neto).
in the literature only interatomic distances are explicitly considered, which are often not suf®cient to account for constraints on angles and are then complemented by linear conditions derived through vector algebra [2]. The present treatment allows a direct introduction of any kind of constraints to the internal motion and does not rely on constraint forces and Lagrange multipliers used in all previous works, see for instance the recent book by Frenkel and Smit [3]. In practice once a set of Cartesian displacements are given for a molecule not subject to kinematical restrictions, the corresponding constrained quantities are obtained, for any choice of the constraints, through a geometric procedure which does not involve forces and time step. 2. Cartesian and molecular displacements For an N atomic molecule, the variables of the motion within a given time step are 3N Cartesian
0022-2860/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0022-286 0(01)00511-7
136
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
displacements, measured in a space ®xed frame, and 3N 2 6 (or 3N 2 5; for linear molecule) independent internal coordinates, all de®ned in the following Table Time
Coordinates
Displacements
XisW
DXis Xis 2 XisW qa Ra 2 RW a
t t 1 Dt Xis
RW a Ra
p Uis mi DXis ) Us
where mi is the mass of the ith atoms, s count components and Ra are intramolecular parameters. Whenever possible we use mass weighted quantities, labelled with a unique subscript ranging from 1 to 3N. Besides, we adopt from now on the summation convention over repeated literal suf®xes. The internal and external (translational and rotational) displacements form a complete set of 3N independent molecular coordinates, all called in general q a. They are partitioned into two subsets, {q n} and {q a }, where a Greek suf®x identi®es the parameters forced to remain constant through the simulation. Actually we reserve Us and q a to indicate unconstrained displacements which, when the molecule undergoes kinematical restrictions, assume different values referred to as UsC and qaC ; with the obvious requirements that qaC ¼; UsC ; ¼ 0
a constrained coordinates 1
The molecular coordinates are characterized by a metric tensor whose contravariant and covariant components evaluated at (W), that is, at the initial time t, are de®ned by the scalar products ! a 2q 2qb 2Us 2Us ab gW gW ab 2Us W 2Us W 2qa W 2qb W W gac W gcb dab
such that
gab W
qb
qm
gW ab
qb
qm
qa
gab W
gaWm
qa
gW ab
gW am
q
gnWb
gnm W
r a qa
W a r n qn 1 g nm W gma q
q
n
gW nb
gW nm
2
nm where g nm W are elements of a submatrix g W inverse of W the block gnm : Since the transformation does not modify the q a s, the new coordinates are also divided into two subset, {r a } and {r n}. It is shown [7] that the W components g ab W and g ab ; of the metric tensor in r space, assume the following block diagonal form g ab W
rb
rm
g W ab
rb
rm
ra rn
gab W 0
0 gnm W
ra rn
g W ab 0
0 gW nm
3
ab where g W ab is the inverse of gW ; that is
W gag gb dab W g
a; b; g constrained coordinates 4
The metric tensor (3) is consistent with, or may be obtained from, a set of relations between ®rst derivatives in r and q space given [7] by a a 2r 2q 2Us 2Us 2Us W 2Us W 2r n W 2qn W 5
2Us 2r a n
2r 2Us
W
W
2r b 2Us
gW ab
gnm W
n
All derivatives here are known [1,4] and satisfy the Eckart±Sayvetz (ES) conditions [5,6]. With the above notation the components of the metric tensor in q space are written in block form as
n
where each block is labelled with the same suf®xes of one of its typical elements. The solution of the equations of the constrained motion is greatly simpli®ed by the introduction of a second set of coordinates, r a, de®ned in Ref. [7] by linear combinations of the q as as
2Us 2r m
!
W
6
W
The subset {r } formally includes the external coordinates which are, however, already factorised in q space since the ES conditions determine vanishing values of all terms gW ma coupling them with internal coordinates. Thus, the external co-ordinates are not affected by the transformation (2) and, in practice, one may forget about them as far as only intramolecular parameters are constrained. Once the constrained Cartesian displacements are determined, one may evaluate the constrained values of all r as which are also given by the direct
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
transformation a 2r rCa UC 1 2Us W s
1 2
22 r a 2Us 2Ut
! W
UsC UtC 1 ¼
7
and all rCa vanish by de®nition. The corresponding inverse expansions are ! 2 2Us 2 U s C n Us r 1 12 rn rm 1 ¼ 8 2r n W C 2r n 2r m W C C but a knowledge of these coef®cients beyond the ®rst order is not required here. We note, however, that the independence of the molecular coordinates requires that 0
22 r a 22 r a 2Us 2Ut 2r a 22 Us 1 9 2Us 2Ut 2r b 2r c 2Us 2r b 2r c 2r b 2r c
thus the second (and higher) coef®cients of the direct and inverse transformation are related. 3. Constrained Cartesian displacements Each in®nitesimal Cartesian displacement in terms of coordinates r a is split into two contributions, for the constrained and unconstrained subsets, as 2Us 2Us a dr 1 dr n 10 dUs 2r a W 2r n W The metric tensor (3) guarantees that, in the in®nitesimal limit, the equations of motion in the subspace of the r a s are completely independent from those relative to the coordinates r ns. Since the latter determine the motion of a constrained molecule, the second term in Eq. (10) represents an in®nitesimal, constrained Cartesian displacement from the initial position, that is 2Us dUsC dr n 11 2r n W
137
equations of motion, 2Us U W FmW 2r m W s 2Us 2Us 2V W W W where Fm F F 2 2qm W 2r m W s 2qm W s 13 Since the kinematical restrictions on the r a s do not affect the potential energy V ¼; qn ; ¼; the generalized components FmW of the force do not change when the accelerations U W are replaced by the s corresponding constrained quantities and one has the identities 2Us Cs W 2Us U W U 2r m W 2r m W s or, multiplying both sides by g W nm and using Eq. (6), n n 2r 2r U W U Cs W 14 2Us W 2Us W s This result and Eq. (12) establish that the generalized components of both velocity and acceleration, at the initial position and in the subspace of the r ns, are unaffected by constraints, hence the same necessarily holds for the corresponding components of the displacements from the initial position, therefore we write n n 2r 2r UsC U 15 2Us W 2Us W s This equation plays a central role in the present treatment and is cast into a more convenient form substituting the dummy index s with t, multiplying both sides by 2Us =2r n W and then writing Pst UtC Pst Ut
16
Taking the time derivatives at (W) of the displacements (10) and (11), because of the independence of the molecular coordinates the following relation holds n n 2r 2r C _ U_ W U 12 2Us W s W 2Us W s
where Pst is given by either of the two equivalent expressions n 2Us 2r Pst n 2r W 2Ut W 17 a 2Us 2r Pst dst 2 2r a W 2Ut W
A similar equation is obtained for the components of the acceleration writing, from the Newton's
These are elements of a 3N square matrix which represents a symmetric, idempotent operator whose
138
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
properties are summarized by the equations 2Ut 2Us Pst d 18 Psv Pvt Pst an 2r a W 2r n W Pst Pts
Pst
2r a 2Ut
W
dan
2r n 2Us
W
19
Eq. (18) rely on the independence of the coordinates while the factorisation of the metric tensor is essential for both Eq. (19), as it may be veri®ed using Eq. (6) and the block form (3). An equation of the kind Eq. (16) for the acceleration was obtained in Ref. [8], where the nature of the projector was recognized but not its origin from an explicit transformation (2). The properties of the projector ensure that Eq. (16) holds when UsC is written as a 2q C Us Us 2 lW Dt2 2Us W a which is the expression adopted in the literature, in the usual Verlet scheme, to derive the constrained Cartesian displacements. It is thus not surprising that the latter may be obtained, without recurring to the Lagrange multipliers, directly from Eq. (16) which is rewritten, using the second form Eq. (17) of the projector, as 2Us UsC Us 2 ga 2r a W 20 a 2r a C where g U 2 Ut 2Ut W t An expression for g a is obtained expanding UtC as done Eq. (8) and in this way the constrained Cartesian displacements become UsC Us 2Us 2 2r a 2
2r a 2Ut
! (
2r a 2Ut
W
! " W
1 2
! W
Ut
22 U t 2r n 2r m
21
!
#) W
rCn rCm 1 ¼
This equation becomes an explicit expression to the second order since rCn and r n coincide to the ®rst order, as it is evident using the identity (15) in the expansion
(7). The process may be continued but a more convenient recurrent formulation is given in the Section 4.
4. Practical algorithm The accuracy of the solution is conveniently discussed in terms of constrained displacements Ush and molecular coordinates rha calculated with them, both coincident to the hth order with the exact values and thus such that UsC Ush 1 O r h11 rCa rha ¼; Ush ; ¼ 1 O r h11 For h 1; 2 Eq. (21) gives ! 2Us 1 ra Us Us 2 2r a W
Us2
! " ! 2Us 2r a Us 2 U 2r a W 2Ut W t ! ! # 2r a 22 U t n m 1 2 2 r r 2Ut W 2r n 2r m W
where we note that from the expression of Us1 one has a a 2r 2r Us r a 1 U1 2Us W 2Us W s 2Us Us1 r n 1 O r 2 2r n W hence using the relation (9) Us2 is transformed into " ! ( ! 2Us 2r a 2 a r 1 U1 Us Us 2 2r a W 2Ut W t #) ! 22 r a 1 1 1 1 2 U U 2Ut 2Uv W t v The term in square brackets is the second order expansion of the constrained value r1a ; calculated for the a th coordinate with the approximate displacements Us1 ; therefore, letting Us0 ; Us and r0a ; r a ;
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
the three successive approximate solutions are 2Us Us1 Us0 2 ra Us0 Us 2r a W 0 2Us Us2 Us1 2 ra 2r a W 1 It may be formally shown, by induction, that this result is generalized into a recurrent relation which, using Eqs. (2), (5) and (6) and reverting to explicit summation, gives the constrained Cartesian (not mass weighted) coordinates as 1 X W 2Rb h11 h Xis Xis 2 Rh 2 R0a 22 g mi ab ab 2Xis W a Here R0a are the values kept ®xed during the simulation and in a preliminary stage known functions Ra ; ¼; Xis ; ¼ are used to calculate the ®rst Cartesian derivatives, at (W), of all constrained distances or angles. They are elements of the B matrix [1] and are also needed to evaluate the symmetric matrix W gab ab ; see Eq. (4). The iteration W and its inverse g then starts with h 0; when Xis0 are the unconstrained coordinates calculated with a suitable numerical integrator. The intramolecular parameters Rha are obtained using the constrained Cartesian coordinates Xish ; accurate to the hth order, and Eq. (22) provides with a better approximation, Xish11 : The process is iterated until the desired accuracy is reached for the constraints. Eq. (22) requires the inversion of gab W and thus resembles the matrix method of Ref. [9] which is, however, valid only for distances. In presence of several constraints nothing prevents to apply Eq. (22) in succession to each individual parameter, repeating the procedure until convergence is reached.
139
In this case Eq. (22) is considerably simpli®ed as it is aa suf®cient to computes single components gW aa 1=gW rather the inverting a matrix. This simple procedure corresponds to the SHAKE algorithm [10,11] but in the version invariably given before, from the earlier works to recent articles [12], instead of gaa W one ®nds the term 2qa 2qa Aaa 2Us 2Us W where the ®rst factor is evaluated at the current, rather than initial, position. These terms are thus computed at each cycle of the iteration, rather then only once, which is a time consuming stage, actually unnecessary according to the present treatment. References [1] E.B. Wilson, J.C. Decius, P.C. Cross, Molecular Vibrations, McGraw-Hill, New York, 1955. [2] G. Ciccotti, M. Ferrario, J.P. Ryckaert, Mol. Phys. 47 (1982) 1252. [3] D. Frenkel, B. Smit, Understanding Molecular Simulations, Academic Press, San Diego, 1996. [4] S. Califano, Vibrational States, Wiley, New York, 1976. [5] C. Eckart, Phys. Rev. 47 (1935) 552. [6] A. Sayvetz, J. Chem. Phys. 7 (1939) 383. [7] N. Neto, M. Muniz-Miranda, G. Sbrana, Spectrochim. Acta 50A (1994) 1317. [8] S.W. de Leeuw, J.W. Perram, H.G. Petersen, J. Stat. Phys. 61 (1990) 1203. [9] G. Ciccotti, J.P. Ryckaert, Comp. Phys. Rep. 4 (1986) 345. [10] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [11] J.P. Ryckaert, Mol. Phys. 55 (1985) 549. [12] G. Ciccotti, M. Ferrario, in: B.J. Berne, G. Ciccotti, D.F. Coker (Eds.), Classical and Quantum Dynamics in Condensed Phase Simulations, World Scienti®c, Singapore, 1998.
www.elsevier.nl/locate/molstruc
A new approach to constrained molecular dynamics Natale Neto* Department of Chemistry, University of Firenze, Via G. Capponi 9, 50121 Firenze, Italy Received 31 August 2000; revised 31 January 2001; accepted 31 January 2001
Abstract The internal motion of a semirigid molecule is described in terms of variations of intramolecular distances, valence angles, out of plane wagging and torsions. Constraints may be imposed to any selected group of these parameters and the reduction of the number of degrees of freedom is realized through a suitable basis change, splitting the space of displacements into two orthogonal subspaces, one of them spanned by the coordinates forced to vanish. The treatment differs from the previous ones, which all rely on the Lagrangian l method and thus dispense with the de®nition of unconstrained coordinates, which do not eventually appear also in the expression of the constrained Cartesian displacements given here. It turns out that the previous shake algorithm of the constrained dynamics include unnecessary calculation at each time step. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Constraints; Molecular dynamics; Vibrational dynamics
1. Introduction There is undoubtedly a close similarity between vibrational and molecular dynamics (MD) as both techniques largely use methods of classical mechanics and deal with small, but not in®nitesimal, displacements from a reference geometry. In fact in MD the motion is partitioned into a sequence of steps, each of them occurring in a very tiny time interval during which all atoms move from the initial positions. The corresponding conformational change of a polyatomic molecule may be then be described in terms of internal coordinates of the kind de®ned in the book of Wilson et al. [1], namely stretching, bending, wagging out of plane, torsion or combinations of them. These coordinates are routinely used in vibrational spectroscopy but not in constrained MD. In fact * Tel.: 139-552-757597; fax: 139-552-476961. E-mail address: [email protected]®.it (N. Neto).
in the literature only interatomic distances are explicitly considered, which are often not suf®cient to account for constraints on angles and are then complemented by linear conditions derived through vector algebra [2]. The present treatment allows a direct introduction of any kind of constraints to the internal motion and does not rely on constraint forces and Lagrange multipliers used in all previous works, see for instance the recent book by Frenkel and Smit [3]. In practice once a set of Cartesian displacements are given for a molecule not subject to kinematical restrictions, the corresponding constrained quantities are obtained, for any choice of the constraints, through a geometric procedure which does not involve forces and time step. 2. Cartesian and molecular displacements For an N atomic molecule, the variables of the motion within a given time step are 3N Cartesian
0022-2860/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0022-286 0(01)00511-7
136
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
displacements, measured in a space ®xed frame, and 3N 2 6 (or 3N 2 5; for linear molecule) independent internal coordinates, all de®ned in the following Table Time
Coordinates
Displacements
XisW
DXis Xis 2 XisW qa Ra 2 RW a
t t 1 Dt Xis
RW a Ra
p Uis mi DXis ) Us
where mi is the mass of the ith atoms, s count components and Ra are intramolecular parameters. Whenever possible we use mass weighted quantities, labelled with a unique subscript ranging from 1 to 3N. Besides, we adopt from now on the summation convention over repeated literal suf®xes. The internal and external (translational and rotational) displacements form a complete set of 3N independent molecular coordinates, all called in general q a. They are partitioned into two subsets, {q n} and {q a }, where a Greek suf®x identi®es the parameters forced to remain constant through the simulation. Actually we reserve Us and q a to indicate unconstrained displacements which, when the molecule undergoes kinematical restrictions, assume different values referred to as UsC and qaC ; with the obvious requirements that qaC ¼; UsC ; ¼ 0
a constrained coordinates 1
The molecular coordinates are characterized by a metric tensor whose contravariant and covariant components evaluated at (W), that is, at the initial time t, are de®ned by the scalar products ! a 2q 2qb 2Us 2Us ab gW gW ab 2Us W 2Us W 2qa W 2qb W W gac W gcb dab
such that
gab W
qb
qm
gW ab
qb
qm
qa
gab W
gaWm
qa
gW ab
gW am
q
gnWb
gnm W
r a qa
W a r n qn 1 g nm W gma q
q
n
gW nb
gW nm
2
nm where g nm W are elements of a submatrix g W inverse of W the block gnm : Since the transformation does not modify the q a s, the new coordinates are also divided into two subset, {r a } and {r n}. It is shown [7] that the W components g ab W and g ab ; of the metric tensor in r space, assume the following block diagonal form g ab W
rb
rm
g W ab
rb
rm
ra rn
gab W 0
0 gnm W
ra rn
g W ab 0
0 gW nm
3
ab where g W ab is the inverse of gW ; that is
W gag gb dab W g
a; b; g constrained coordinates 4
The metric tensor (3) is consistent with, or may be obtained from, a set of relations between ®rst derivatives in r and q space given [7] by a a 2r 2q 2Us 2Us 2Us W 2Us W 2r n W 2qn W 5
2Us 2r a n
2r 2Us
W
W
2r b 2Us
gW ab
gnm W
n
All derivatives here are known [1,4] and satisfy the Eckart±Sayvetz (ES) conditions [5,6]. With the above notation the components of the metric tensor in q space are written in block form as
n
where each block is labelled with the same suf®xes of one of its typical elements. The solution of the equations of the constrained motion is greatly simpli®ed by the introduction of a second set of coordinates, r a, de®ned in Ref. [7] by linear combinations of the q as as
2Us 2r m
!
W
6
W
The subset {r } formally includes the external coordinates which are, however, already factorised in q space since the ES conditions determine vanishing values of all terms gW ma coupling them with internal coordinates. Thus, the external co-ordinates are not affected by the transformation (2) and, in practice, one may forget about them as far as only intramolecular parameters are constrained. Once the constrained Cartesian displacements are determined, one may evaluate the constrained values of all r as which are also given by the direct
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
transformation a 2r rCa UC 1 2Us W s
1 2
22 r a 2Us 2Ut
! W
UsC UtC 1 ¼
7
and all rCa vanish by de®nition. The corresponding inverse expansions are ! 2 2Us 2 U s C n Us r 1 12 rn rm 1 ¼ 8 2r n W C 2r n 2r m W C C but a knowledge of these coef®cients beyond the ®rst order is not required here. We note, however, that the independence of the molecular coordinates requires that 0
22 r a 22 r a 2Us 2Ut 2r a 22 Us 1 9 2Us 2Ut 2r b 2r c 2Us 2r b 2r c 2r b 2r c
thus the second (and higher) coef®cients of the direct and inverse transformation are related. 3. Constrained Cartesian displacements Each in®nitesimal Cartesian displacement in terms of coordinates r a is split into two contributions, for the constrained and unconstrained subsets, as 2Us 2Us a dr 1 dr n 10 dUs 2r a W 2r n W The metric tensor (3) guarantees that, in the in®nitesimal limit, the equations of motion in the subspace of the r a s are completely independent from those relative to the coordinates r ns. Since the latter determine the motion of a constrained molecule, the second term in Eq. (10) represents an in®nitesimal, constrained Cartesian displacement from the initial position, that is 2Us dUsC dr n 11 2r n W
137
equations of motion, 2Us U W FmW 2r m W s 2Us 2Us 2V W W W where Fm F F 2 2qm W 2r m W s 2qm W s 13 Since the kinematical restrictions on the r a s do not affect the potential energy V ¼; qn ; ¼; the generalized components FmW of the force do not change when the accelerations U W are replaced by the s corresponding constrained quantities and one has the identities 2Us Cs W 2Us U W U 2r m W 2r m W s or, multiplying both sides by g W nm and using Eq. (6), n n 2r 2r U W U Cs W 14 2Us W 2Us W s This result and Eq. (12) establish that the generalized components of both velocity and acceleration, at the initial position and in the subspace of the r ns, are unaffected by constraints, hence the same necessarily holds for the corresponding components of the displacements from the initial position, therefore we write n n 2r 2r UsC U 15 2Us W 2Us W s This equation plays a central role in the present treatment and is cast into a more convenient form substituting the dummy index s with t, multiplying both sides by 2Us =2r n W and then writing Pst UtC Pst Ut
16
Taking the time derivatives at (W) of the displacements (10) and (11), because of the independence of the molecular coordinates the following relation holds n n 2r 2r C _ U_ W U 12 2Us W s W 2Us W s
where Pst is given by either of the two equivalent expressions n 2Us 2r Pst n 2r W 2Ut W 17 a 2Us 2r Pst dst 2 2r a W 2Ut W
A similar equation is obtained for the components of the acceleration writing, from the Newton's
These are elements of a 3N square matrix which represents a symmetric, idempotent operator whose
138
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
properties are summarized by the equations 2Ut 2Us Pst d 18 Psv Pvt Pst an 2r a W 2r n W Pst Pts
Pst
2r a 2Ut
W
dan
2r n 2Us
W
19
Eq. (18) rely on the independence of the coordinates while the factorisation of the metric tensor is essential for both Eq. (19), as it may be veri®ed using Eq. (6) and the block form (3). An equation of the kind Eq. (16) for the acceleration was obtained in Ref. [8], where the nature of the projector was recognized but not its origin from an explicit transformation (2). The properties of the projector ensure that Eq. (16) holds when UsC is written as a 2q C Us Us 2 lW Dt2 2Us W a which is the expression adopted in the literature, in the usual Verlet scheme, to derive the constrained Cartesian displacements. It is thus not surprising that the latter may be obtained, without recurring to the Lagrange multipliers, directly from Eq. (16) which is rewritten, using the second form Eq. (17) of the projector, as 2Us UsC Us 2 ga 2r a W 20 a 2r a C where g U 2 Ut 2Ut W t An expression for g a is obtained expanding UtC as done Eq. (8) and in this way the constrained Cartesian displacements become UsC Us 2Us 2 2r a 2
2r a 2Ut
! (
2r a 2Ut
W
! " W
1 2
! W
Ut
22 U t 2r n 2r m
21
!
#) W
rCn rCm 1 ¼
This equation becomes an explicit expression to the second order since rCn and r n coincide to the ®rst order, as it is evident using the identity (15) in the expansion
(7). The process may be continued but a more convenient recurrent formulation is given in the Section 4.
4. Practical algorithm The accuracy of the solution is conveniently discussed in terms of constrained displacements Ush and molecular coordinates rha calculated with them, both coincident to the hth order with the exact values and thus such that UsC Ush 1 O r h11 rCa rha ¼; Ush ; ¼ 1 O r h11 For h 1; 2 Eq. (21) gives ! 2Us 1 ra Us Us 2 2r a W
Us2
! " ! 2Us 2r a Us 2 U 2r a W 2Ut W t ! ! # 2r a 22 U t n m 1 2 2 r r 2Ut W 2r n 2r m W
where we note that from the expression of Us1 one has a a 2r 2r Us r a 1 U1 2Us W 2Us W s 2Us Us1 r n 1 O r 2 2r n W hence using the relation (9) Us2 is transformed into " ! ( ! 2Us 2r a 2 a r 1 U1 Us Us 2 2r a W 2Ut W t #) ! 22 r a 1 1 1 1 2 U U 2Ut 2Uv W t v The term in square brackets is the second order expansion of the constrained value r1a ; calculated for the a th coordinate with the approximate displacements Us1 ; therefore, letting Us0 ; Us and r0a ; r a ;
N. Neto / Journal of Molecular Structure 563±564 (2001) 135±139
the three successive approximate solutions are 2Us Us1 Us0 2 ra Us0 Us 2r a W 0 2Us Us2 Us1 2 ra 2r a W 1 It may be formally shown, by induction, that this result is generalized into a recurrent relation which, using Eqs. (2), (5) and (6) and reverting to explicit summation, gives the constrained Cartesian (not mass weighted) coordinates as 1 X W 2Rb h11 h Xis Xis 2 Rh 2 R0a 22 g mi ab ab 2Xis W a Here R0a are the values kept ®xed during the simulation and in a preliminary stage known functions Ra ; ¼; Xis ; ¼ are used to calculate the ®rst Cartesian derivatives, at (W), of all constrained distances or angles. They are elements of the B matrix [1] and are also needed to evaluate the symmetric matrix W gab ab ; see Eq. (4). The iteration W and its inverse g then starts with h 0; when Xis0 are the unconstrained coordinates calculated with a suitable numerical integrator. The intramolecular parameters Rha are obtained using the constrained Cartesian coordinates Xish ; accurate to the hth order, and Eq. (22) provides with a better approximation, Xish11 : The process is iterated until the desired accuracy is reached for the constraints. Eq. (22) requires the inversion of gab W and thus resembles the matrix method of Ref. [9] which is, however, valid only for distances. In presence of several constraints nothing prevents to apply Eq. (22) in succession to each individual parameter, repeating the procedure until convergence is reached.
139
In this case Eq. (22) is considerably simpli®ed as it is aa suf®cient to computes single components gW aa 1=gW rather the inverting a matrix. This simple procedure corresponds to the SHAKE algorithm [10,11] but in the version invariably given before, from the earlier works to recent articles [12], instead of gaa W one ®nds the term 2qa 2qa Aaa 2Us 2Us W where the ®rst factor is evaluated at the current, rather than initial, position. These terms are thus computed at each cycle of the iteration, rather then only once, which is a time consuming stage, actually unnecessary according to the present treatment. References [1] E.B. Wilson, J.C. Decius, P.C. Cross, Molecular Vibrations, McGraw-Hill, New York, 1955. [2] G. Ciccotti, M. Ferrario, J.P. Ryckaert, Mol. Phys. 47 (1982) 1252. [3] D. Frenkel, B. Smit, Understanding Molecular Simulations, Academic Press, San Diego, 1996. [4] S. Califano, Vibrational States, Wiley, New York, 1976. [5] C. Eckart, Phys. Rev. 47 (1935) 552. [6] A. Sayvetz, J. Chem. Phys. 7 (1939) 383. [7] N. Neto, M. Muniz-Miranda, G. Sbrana, Spectrochim. Acta 50A (1994) 1317. [8] S.W. de Leeuw, J.W. Perram, H.G. Petersen, J. Stat. Phys. 61 (1990) 1203. [9] G. Ciccotti, J.P. Ryckaert, Comp. Phys. Rep. 4 (1986) 345. [10] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [11] J.P. Ryckaert, Mol. Phys. 55 (1985) 549. [12] G. Ciccotti, M. Ferrario, in: B.J. Berne, G. Ciccotti, D.F. Coker (Eds.), Classical and Quantum Dynamics in Condensed Phase Simulations, World Scienti®c, Singapore, 1998.