A new formalism for molecular dynamics in internal coordinates

Chemical Physics 265 (2001) 63±85

www.elsevier.nl/locate/chemphys

A new formalism for molecular dynamics in internal coordinates Sang-Ho Lee, Kim Palmo, Samuel Krimm * Biophysics Research Division, Department of Physics, The University of Michigan, 930 N. University Avenue, Ann Arbor, MI 48109, USA Received 5 September 2000

Abstract Internal coordinate molecular dynamics (ICMD) has been used in the past in simulations for large molecules as an alternative way of increasing step size with a reduced operational dimension that is not achievable by MD in Cartesian coordinates. A new ICMD formalism for ¯exible molecular systems is presented, which is based on the spectroscopic Bmatrix rather than the A-matrix of previous methods. The proposed formalism does not require an inversion of a large matrix as in the recursive formulations based on robot dynamics, and takes advantage of the sparsity of the B-matrix, ensuring computational eciency for ¯exible molecules. Each moleculeÕs external rotations about an arbitrary atom center, which may di€er from its center of mass, are parameterized by the SU(2) Euler representation, giving singularity free parameterization. Although the formalism is based on the use of nonredundant generalized (internal and external) coordinates, an MD simulation in linearly dependent coordinates can be done by ®nding a transformation to a new set of independent coordinates. Based on the clear separability in the generalized coordinates between fast varying degrees of freedom and slowly varying ones, a multiple time step algorithm is introduced that avoids the previous nontrivial interaction distance classi®cation. Also presented is a recursive method for computing nonzero A-matrix elements that is much easier to apply to a general molecular structure than the previous method. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction The equation of motion of a molecular system has been described in mainly two di€erent coordinate frames. The conventional formalism for molecular dynamics (MD) simulation [1,2] has been developed in the Cartesian coordinate frame since it provides the simplest form for the systemÕs equation of motion, in which each atom is treated as a point mass whose position changes subject to the imposed external force. However, it has been

*

Corresponding author. Fax: +1-734-764-3323. E-mail address: [email protected] (S. Krimm).

noted that in this formalism the step size for numerical integrations should be kept small (0.5±1 fs) in order to maintain the systemÕs stability. This makes it almost impossible to routinely simulate in a time range longer than ns, where interesting conformational changes in biomolecules may take place. Since the instability arises mainly from rapidly varying motions like bond stretching (e.g., a typical CH stretch frequency of 3000 cm 1 corresponds roughly to a period of 10 fs) and angle bending, one way of increasing step size is to remove the fast degrees of freedom by introducing suitable constraints on interatomic distances [3±9]. In the commonly used algorithms like SHAKE [3] and RATTLE [4], the systemÕs initial states are

0301-0104/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 ( 0 1 ) 0 0 2 3 6 - 1

64

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

estimated from the unconstrained part of the equation of motion and they are corrected iteratively so as to satisfy the desired constraints to within acceptable error limits. Since the rapidly varying motions are related to bonded forces or nonbonded forces within a short distance range while nonbonded forces of long interatomic distances change slowly, another method adopts multiple time steps (MTS) [10±13] depending on the interaction distance class: a small time step for a class of fast varying forces while a long time step for a class of slowly varying forces. However, the slow force calculated at an interval of its time step may produce a resonance leading to instability [14]. In an improved algorithm, MOLLY [15], the slow forces are computed at averaged positions of the fast propagations. As a third method to increase the time step, an implicit Euler scheme is incorporated with a Langevin equation approach [16], which is introduced to e€ectively establish the systemÕs thermal equilibrium. A signi®cant improvement in this direction is made by additionally adopting a normal mode technique [17], in which in the ®rst stage the potential energy is approximated to be harmonic, which makes the concerned Langevin equation linear, and in the next stage the analytic solution of this is further corrected by the implicit Euler scheme. Alternatively, the systemÕs equation of motion can be described in internal coordinates, which are de®ned from local atomic connectivities like bond lengths, bond angles, and proper torsions. Internal coordinates have already proven to be e€ective in ab initio geometry optimizations [18,19] and in Monte Carlo simulations [20]. Although the resulting equation of motion is far more complex and has been less used in MD studies than that in Cartesian coordinates, the internal coordinate MD (ICMD) formalism [21±31] is not only important from a theoretical point of view but also has unique and interesting features, which may also prove it to be e€ective in simulations for large molecules or many-molecule systems. Since internal coordinates are clearly separable from external rotations and translations, thermodynamic quantities of a system may be better determined than in Cartesian coordinates. (At a given temperature, a systemÕs volume and pressure depend mostly on

external translations of each molecule.) In terms of internal coordinates, the fast varying degrees of freedom are clearly separable from the slowly varying ones, and a signi®cant reduction in operational dimension can be achieved by simply neglecting fast varying degrees of freedom without introducing the constraints that are necessary to freeze such coordinates, as in SHAKE [3] or RATTLE [4]. In fact, the motions related to rapidly varying degrees of freedom are highly localized [17] and their positional displacements are closer to Gaussian type ¯uctuations than to signi®cant conformational changes [32]. Thus, in the early studies of conformational energies of polypeptides, G o and Scheraga [21] used only torsions as initial variable coordinates to ®nd an approximate minimum energy conformation. Their idea has been further developed into MD formulations in torsion angle coordinates [22]. Mazur and Abagyan [23,24] extended this method to include bond stretches and angle bend coordinates, while Kneller and Hinsen [29] incorporated quaternion parameters and angular velocities for rotations of linked rigid body subunits. However, all these formulations require ®nding the inverse of a massmatrix at each time step (or at least solving a system of linear equations), which is the costliest part in the whole process. As a result, attention has been given to recursive algorithms [26±28] that are based on multi-arm robot dynamics [33,34] and avoid the direct inversion of the large mass-matrix by using equivalent local body-level equations of motion. We refer to these previous ICMD formalisms as A-matrix formulations, since they are based on the mass-matrix constructed from atomic masses and the A-matrix that transforms internal and external coordinates into Cartesian coordinates. We present another ICMD formalism, called a B-matrix formalism, which is based on the spectroscopic B-matrix [35±37]. As long as each molecule is considered to be entirely ¯exible with no rigid constraints, the spectroscopic B-matrix is well-de®ned and for each internal coordinate the corresponding B-matrix elements are nonzero only for a few (at most four) related atoms. This signi®cantly reduces the required arithmetic computations. Furthermore, the B-matrix formalism

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

allows the whole process to be executed without any inversion of a large matrix as in the fast recursive algorithms [26±28]. In this respect, the proposed B-matrix formalism may be a very ecient method for MD simulations of large or many-molecule systems. We have also developed a singularity free formulation for external rotations that uses the Euler parameters [38], equivalent to the quaternion parameters, of an SU…2† representation. The explicit separation between rapidly varying degrees of freedom and slowly varying ones also enables us to e€ectively incorporate the MTS method [12] in ICMD, the interaction distance classi®cation required in the Cartesian coordinate formulation not being necessary. The imposition of rigid constraints on some internal coordinates can also be handled within the B-matrix formalism, as we show. This requires solving a linear equation to determine the corresponding Lagrange undetermined multipliers. Although this may be less e€ective in reducing operational dimensions than the A-matrix formalisms, the overall eciency of this approach may be greater, although this remains to be determined. Alternatively, as an approximate way to avoid solving such a linear equation, we can e€ectively neglect the desired fast varying coordinates as in the A-matrix formalisms. This approximation of simply neglecting some fast varying coordinates without solving the linear equation might be an e€ective way of ®nding a path from an arbitrary initial system con®guration to an approximate equilibrium con®guration where the conformational energy of the molecule of interest is close to its minimum. Near the equilibrium con®guration we can do full B-matrix ICMD without any rigid constraint in order to derive desired physical quantities. From an arbitrary initial conformation of the known amino acid sequence of a protein, ®nding folding pathways to its observed structure or reliably predicting its stable native conformation may be one of the most interesting objectives of MD simulations. The successful accomplishment of this will ultimately depend on the full reliability of the force ®eld [39] and on the ability to carry out extended MD simulations [40]. In combination with our general type of molecular mechanics (MM)

65

energy function, called a spectroscopically determined force ®eld (SDFF) [41±43], which is developed to reproduce high level ab inito results and observed vibrational (infrared and Raman) frequencies, the B-matrix ICMD formalism presented in this paper may provide an optimal way to the ultimate goal. The conventional Cartesian coordinate MD (CCMD) formalism is summarized in Section 2. In Section 3, key elements of the conventional A-matrix ICMD formalism are summarized and extended to general dependent coordinates. A more ¯exible and systematic way of calculating the nonzero A-matrix elements is proposed, although it is associated with a few more computational steps than the previous method [22,23]. In Section 4, the B-matrix ICMD formalism is introduced with the tools for computing all the required elements. The external rotations are also represented in singularity free Euler parameters of SU…2†. In Section 5, an MTS algorithm in internal and external coordinates is introduced. 2. Dynamics in Cartesian coordinates In a general formulation, we consider a system of Nmol molecules that are all ¯exible. Let the cth molecule contain pc atoms and its kth atom have atomic mass mck and position xkc with respect to (w.r.t.) an arbitrary laboratory-®xed frame (LFF) with ^1  …1; 0; 0†, ^2  …0; 1; 0†, and ^3  …0; 0; 1† being its basis. In a classical approximation, with each bond in a molecule being considered to be massless, the systemÕs total kinetic energy T is found to be 2T ˆ

Nmol X pc X 3 X cˆ1 kˆ1 kˆ1

k

k

mck x_ kc x_ kc

…1†

where k is the index of a Cartesian coordinate component and a dot represents a derivative w.r.t. time. If the systemÕs total potential energy V is not k an explicit function of atomic velocities x_ kc , then the motion of the kth atom in the cth molecule is subject to Newton's equation: kc ˆ f kc mck x

…2†

66

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

with f kc  oV =oxkc being the force acting on the atom at xkc due to V. The potential energy is usually expandable in the neighborhood of an arbitrary system con®guration by " # Nmol X pc X 3 X k k oV V ˆ Vo ‡ …xkc xkoc † k oxkc o cˆ1 kˆ1 kˆ1 3 c X b X mol X 1X j …xmb 2 b;cˆ1 mˆ1 kˆ1 j;kˆ1 " # o2 V ‡


 k j oxmb oxkc o N

p

p

‡

2

k

j

k

xmob †…xkc

xkoc † …3†

k j =oxmb oxkc Šo

where the ‰o V de®ne the systemÕs Hessian matrix. In the cases of MM or MD simulations, the total potential energy is usually separable into two parts: V ˆVr ‡VS

a X b b mol X mol X 1X 1X lm Vablm …rab †‡ V lm …rlm † 2 a6ˆb lˆ1 mˆ1 2 bˆ1 l6ˆm bb bb

N

p

p

N

p

…5†

where Vablm …ˆ Vbaml † and Vbblm represent the inter- and intra-molecular nonbonded interaction potentials between atom l in the ath molecule and atom m in the bth molecule, respectively. Since Vablm is a lm function of the distance rab de®ned by m lm rab  jxlm ab j  xb

xl

…6†

a

where xlm atom l to atom m, the ab is the vector from k derivatives of V r w.r.t. xkc are found to be

ˆ

!

 lk Nmol X pa X o oVac

ˆ

k elk lk ac orac

lkj a6ˆc lˆ1 oxac

‡

 lk pc X o oVcc

k elk lk cc orcc

j

l6ˆk

oxlk cc

 …dbc dkm



…dkm

p

Nmol X pa X oV lk ac lk or ac a6ˆc lˆ1

p

k

elk ac ‡

pc X oV lk cc lk or cc l6ˆk

c Scc ˆ Soc ‡

‡

pc X 3 X kˆ1 kˆ1

k

k

…7†

dlm †dbc

…8†

k

Bcckk …xkc

pc X 3 1X

xkoc † j

‰Bc Šc k …xm 2 m;kˆ1j;kˆ1 2 mj k c

j

k

xmoc †…xkc

k

xkoc † ‡


…9†

c ‰Bc2 Šmj kk

Bcckk

and are called the where the coecients ®rst- and the second-order B-matrix elements, respectively. It is implicitly assumed in the expression of Eq. (9) that changes in Scc are due only to changes in atomic coordinates in the cth molecule. When there are strong intermolecular (e.g., hydrogen bond) interactions and some internal coordinates must be de®ned based on atomic coordinates that belong to di€erent molecules, we should consider the group of related molecules as a super-molecule that is formed by virtual bonds between the related molecules. The derivatives of k V S w.r.t. xkc are found to be mol X oV S X oV S oSbb X oV S b ˆ ˆ B ; k b b ckk kk oxkc bˆ1 b oSb oxc b oSc

N

o2 V S k ˆ j oxmb oxkc

X o2 V S b;c

oScc oScb

X oV S b

elk cc ;

dba dlm †

b lk lk where elk ac  xac =rac and dc ˆ 0 for b 6ˆ c and b dc ˆ 1 for b ˆ c. In the neighborhood of an arbitrary expansion center xkoc , the cth internal coordinate Scc can be expanded by

‡

a c mol X oV r X oVaclk X oVcclk ˆ ‡ k lkk lkk oxkc a6ˆc lˆ1 oxac l6ˆk oxcc

N

oV r k oxkc

…4†

where V r is due to nonbonded (van der Waals or electrostatic) interactions and V S represents bonded interaction terms. Typically, V S is expressed by each moleculeÕs internal coordinates Scc , while V r can be put in the form of Vr ˆ

o2 V r o j k  mj k oxmb oxb oxc

…10†

Bccmj Bbckk

b ‰Bc2 Šmj kk oScb

! dbc :

…11†

Once f kc is computed from Eqs. (7) and (10) at some instant t, which is the costliest part in CCMD, the new positions and velocities after a time step Dt in the trajectory (following the ve-

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

locity Verlet propagation scheme) are determined by xkc …t ‡ Dt† ˆ xkc …t† ‡ Dtx_ kc …t† ‡ x_ kc …t ‡ Dt† ˆ x_ kc …t† ‡

Dt2 k f …t†; 2mck c

…12†

Dt  k f c …t† ‡ f kc …t ‡ Dt† : c 2mk …13†

3. Internal coordinate dynamics with A-matrix 3.1. Basic formulation A molecule's internal coordinates S c can be constructed from arbitrary functions of the primitive internals nk , which are bond lengths, bond angles, or proper torsion angles based on the atomic connectivity in the molecule. (Since the entire systemÕs equations of motion can be determined from that of each molecule, we concentrate on a speci®c molecule, viz., the cth molecule, and drop the molecular index if it causes no confusion.) The number, n, of all the primitive internals may exceed the molecule's internal degrees of freedom f, in which case some of them are connected through n±f redundancy (linear or nonlinear) relations. The changes in internal coordinates due to in®nitesimal displacements in the primitive internals can be expressed by n n X 1X c DS c ˆ Ukc Dnk ‡ Umk Dnm Dnk ‡


2 kˆ1 m;kˆ1

…14†

c  o2 S c =onm onk . Although there are many with Umk ways of de®ning a molecule's internal coordinates, we consider only such internal coordinates that give a nonsingular inverse relation of Eq. (14) as

Dnk ˆ

n n X 1X Vck DS c ‡ Vmck DS m DS c ‡


: 2 cˆ1 m;cˆ1

…15†

This enables us to calculate values in the primitive internals nk from the speci®ed values in internal coordinates S c , which may involve iterative processes if we only use terms to the ®rst order. The atomic coordinates in an arbitrary coordinate system ®xed to the molecule can be obtained from

67

the primitive internals by constructing a Z-matrix, which is commonly used in setting up ab initio geometry optimization in internal coordinates [18,19]. Not all the primitive internals are necessary to construct a Z-matrix but only a suitable set of f independent primitive internals are required in order to de®ne all atomic coordinates in a molecule-®xed frame (MFF). We refer to the primitive internals in this independent subset as Z-coordinates Rk , which are orthogonal to each other. For convenience, we put the primitive internals in the order of  k R ; for 1 6 k 6 f ; k n ˆ …16† nk ; for f < k 6 n: In order to determine atomic coordinates in an MFF, we consider an arbitrary molecule (which may be constructed from some virtual atoms or virtual bonds) as in Fig. 1, where each atom is indexed according to the order of the step in constructing a desired Z-matrix. Since a molecule's center of mass is obtained only after determining its atomic coordinates, instead of using the center of mass frame (Appendix B), we select the atom f…ˆ 1† as the origin (or base) of the MFF with its basis being f^1m ; ^2m ; ^3m g, which is arbitrarily de®ned by ^1m  e12m ˆ ^3m 

e21m ;

u13m 

…17†

…e21m  e23m †= sin h13

ˆ …e12m  e23m †= sin h13 ˆ u31m ; ^2m  ^3m  ^1m ˆ e12m  u13m :

…18† …19†

The ®rst three atomic coordinates in f^1m ; ^2m ; ^3m g are easily found to be z1m ˆ …0; 0; 0†m ;

…20†

z2m ˆ r12 e12m ˆ …r12 ; 0; 0†m ;

…21†

z3m ˆ z2m

r23 cosh13 e12m ‡ r23 sin h13 e12m  u13m : …22†

Thus, the coordinates of atom d, which can be de®ned from the Z-matrix input representation of

68

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

Fig. 1. Orientation of an arbitrary molecule with respect to a laboratory-®xed coordinate frame.

d, k, rkd , m, hmd , l, and sld , are to be obtained by the trigonometric relation zdm ˆ zkm

rkd cos hmd emkm ‡ rkd sin hmd

 cossld emkm  ulkm


sin sld ulkm ;

…23†

where emkm , ulkm , and emkm  ulkm can be computed from the already determined zlm , zmm , and zkm in the previous steps. In this way, all atomic coordinates in the MFF can be obtained recursively from the ones determined in the previous steps. Since the position of the fth atom is taken as the origin of the MFF, it is evident that the position vector of atom f in f^ 1; ^ 2; ^ 3g itself is the translation vector from the LFF to the MFF. Therefore, the atomic coordinates, xk , in f^ 1; ^ 2; ^ 3g can be determined from zkm by xk ˆ xf ‡ Nzkm

…24†

where N is the relative rotation matrix between the two frames. Just for convenience, we take the external rotation matrix N to be parameterized as [44] N…~ /†  N…/1 ; /2 ; /3 †  exp …/3 D3 † exp…/2 D2 † exp …/1 D1 † with

…25†

0

0

0

0

1

B C 1 A; D1  @ 0 0 0 1 0 0 1 0 0 1 B C D2  @ 0 0 0 A; 1 0 0 0 1 0 1 0 B C D3  @ 1 0 0 A: 0 0 0

…26†

Thus, the matrix N…0† of zero rotational parameters is just an identity matrix, which implies that the MFF is related to the LFF by translations only. The coordinate vector xd in f^1; ^2; ^3g, which is equivalent to Eq. (23), is found to be xd ˆ xk rkd cos hmd emk ‡ rkd sin hmd cossld emk  ulk


sin sld ulk : …27†

Eqs. (24)±(27) enable us to naturally derive a relation that is complementary to Eq. (9), X k xkk ˆ xkok ‡ Akc …S c Soc † c

1X kk ‡ ‰A2 Šbc …S b 2 b;c

Sob †…S c

Soc † ‡


;

…28†

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85 k

k

k

where the coecients Akc and ‰A2 Šbc are called the ®rst- and the second-order A-matrix elements, respectively. The nonzero A-matrix elements for internal coordinates can be computed from the direct derivatives of Eq. (27) w.r.t. the Z-coordinates, or from simple relations related to vector rotations [22,23], which will be further investigated in this section. It should be noted that Eq. (28) is well-de®ned also for the external rotations and translations, with related A-matrix elements being computed from the direct derivatives of Eq. (24). When S c represent not only the internal coordinates but also the molecule's external rotations and translations, we will designate the S c as generalized coordinates. From Eq. (28) the velocity of the kth atom is found to be X k x_ kk ˆ Akc S_ c : …29† c

Substituting this into Eq. (1), we have for the systemÕs kinetic energy 2T ˆ

Nmol X X cˆ1 b;c

gcbc S_ cb S_ cc ˆ

Nmol X cˆ1

S_ Tc gc S_ c

…30†

with gcbc 

p X 3 X kˆ1 kˆ1

  k k mck Akcb Akcc ˆ ATc uc 1 Ac bc ;

…31†

where uc is a diagonal matrix containing triads of the inverse atomic masses and ATc is the transpose matrix of Ac . The systemÕs Euler±Lagrange equation of motion w.r.t. the generalized coordinates Scc leads to X c

gcbc Scc ‡

X gcbcd S_ cc S_ cd ˆ c;d

oV oScb

…32†

with  ogcbd ogcbc ogccd ‡ oScc oScd oScb pc X 3 n k k k 1X k mck Akcb ‰Ac2 Šk‡cd ‡ Akcc ‰Ac2 Šk db ˆ 2 kˆ1 kˆ1 o k kk ‡ Akcd ‰Ac2 Š cb ;

gcbcd 

1 2



…33†

69

where kk

kk

kk

‰Ac2 Š‡cd  ‰Ac2 Šcd ‡ ‰Ac2 Šdc k ‰Ac2 Šk cd



k ‰Ac2 Škcd

and

k ‰Ac2 Škdc :

kk

Note that the values of ‰Ac2 Š cd are zero for all internal coordinates but are nonzero for external rotations due to the noncommutativity of external rotations. The derivative of V w.r.t. Scb is found to be oV oScb

ˆ

Nmol X pa X pc X 3 X oV lm ac lmj a6ˆc lˆ1 mˆ1 jˆ1 oxac

‡

pc X 3 X oV lm cc lmj l6ˆm jˆ1 oxcc

j

j

Amcb

Amcb ‡

oV S oScb

:

…34†

A usual way of solving Eq. (32) requires inverting the matrix gc , which is also called a massmatrix since the conjugate momentum of Scc is given by X oT Jcc  c ˆ gccd S_ cd : …35† oS_ c d Although gc is positive de®nite for linearly independent generalized coordinates, it is singular when there is a redundancy relation in Scc , which requires a special treatment to solve Eq. (32). Since gc  ATc uc 1 Ac is symmetric, there is an orthogonal
matrix PTc  KTc LTc such that  
Kc T Pc gc Pc  …ATc uc 1 Ac † KTc LTc Lc   Cc 0 ; …36† ˆ 0 0 where Cc is a nonsingular diagonal matrix of order 3pc . The corresponding generalized (Moore±Penrose) inverse of gc is de®ned so that  1  0 Cc T T T 1 Pc gc Pc  Pc …Ac uc Ac † Pc ˆ : …37† 0 0 In fact, it is of no use to apply gc directly to the left-hand side of Eq. (32). Note that Lc Sc de®nes a set of redundant coordinates that remain constant under an in®nitesimal displacement in the molecule, while Kc Sc de®nes a set of nonredundant coordinates. If we de®ne Sc  Kc Sc , Ac  Ac KTc ,

70

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

gc  ATc uc 1 Ac …ˆ Kc gc KTc †, and gc  Kc gc KTc , then from Eqs. (36) and (37) we have gc gc ˆ 13pc 3pc ˆ gc gc :

…38†

This shows that gc is the inverse matrix of gc , and we denote its elements by gbc with raised indices. c With the consideration of  ˆ PT Pc g PT Pc S  ˆ KT g Kc S;  gc S c c c c c

…39†

if we apply gc 1 Kc to the left-hand side of Eq. (32)   Kc S  c , and Ac  and de®ne S_ c  Kc S_ c , S c 2 c T Kc A2 Kc , then the systemÕs equation of motion reduces to a Sc ‡

3pc X c;dˆ1

c d

gaccd S_ c S_ c ˆ f ac

…40†

with gaccd 

f ac 

3pc X gab gc ; c bcd

…41†

bˆ1

3pc h X bˆ1

gc 1 Kc

iab oV oScb

ˆ

3pc X oV gab : b c oS c bˆ1

…42†

Thus, the equivalent to Eqs. (12) and (13) can be formulated as 2

Dt a a S ac …t ‡ Dt† ˆ S ac …t† ‡ DtS_ c …t† ‡ S …t†; 2 c

…43†

o Dt n a a a a S_ c …t ‡ Dt† ˆ S_ c …t† ‡ S c …t† ‡ Sc …t ‡ Dt† : 2

…44†

At each time step, once new values of all Sc , S_ c ,  are determined, values in redundant coorand S c dinates are obtained by Sc ˆ KTc Sc , S_ c ˆ KTc S_ c , and   c ˆ KT S S c c. For nonredundant coordinates, the inverse matrix of gc is well-de®ned and the transformation to a new set of nonredundant coordinates is not necessary, and we can simply remove the underlines in Eqs. (40)±(42). We assume the generalized coordinates Scc to be nonredundant if not otherwise noted. In order to advance Scc along the time tra-

jectory by Eq. (40), however, it is necessary at each time step to compute not only oV =oScb and gcbcd , which involve computations of the ®rst- and the second-order A-matrices, but also gc , which is known to be the costliest part. The whole process is so complex and time consuming that ICMD, compared to CCMD, has rarely been used thus far and the attention that has been given to it has been devoted to developing fast algorithms to compute A-matrices and gc . Although a signi®cant achievement was made in this direction by adopting, instead of the direct matrix inversion of gc , the recursive formulation [26±28], it will not be pursued here. One of the big advantages obtainable from an ICMD simulation is the fact that we can signi®cantly reduce the operational dimension of each molecule by simply neglecting sti€ and rapidly varying degrees of freedom such as valence bond length and bond angle coordinates in Eqs. (29) and (32), without introducing any constraint equation as is necessary in the equivalent treatments like SHAKE [3] or RATTLE [4] in a CCMD formulation. In real molecular systems, valence bond length and bond angle coordinates of each molecule do not deviate much from their average values during a ®nite time, which may be much longer than the time step, while torsion angles may undergo signi®cant changes [21]. Therefore, the early formulation of ICMD to ®nd the minimum energy conformations of macromolecules was developed by taking only torsional degrees of freedom into consideration with other degrees of freedom being kept rigid [22,25], which simpli®es the calculations of A, A2 , and gc . 3.2. Computation of A-matrix elements for internal coordinates The method of computing A-matrix elements in the previous ICMD formalisms is based on the observation that (for a speci®c illustration see Fig. 1) the changes in x6 due to in®nitesimal changes in h25 and s15 can be viewed as in®nitesimal rotations of x36 about u25 and e23 , respectively [22,23]: ox6 ˆ u25  …x6 oh25

x3 †;

…45†

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

ox6 ˆ e23  …x6 os15

x3 †;

…46†

and the derivative of x6 w.r.t. r23 can be obtained from ox6 o ˆ …x2 ‡ x23 ‡ x35 ‡ x56 † ˆ e23 : or23 or23

…47†

Values of x36 , u25 , and e23 are easily computable from the de®ned atomic coordinates up to the stage for atom 6. Similarly, with Rc involving all Zcoordinates corresponding to stages from atom 6 to atom 1 (the base), values of ox6 =oRc can be determined. Otherwise, ox6 =oRc are all zero. (The derivatives of x6 w.r.t. the primitive internals not yet de®ned to stage 6 or belonging to di€erent branches from the base are all zero.) In this way all the nonzero ®rst- and second-order A-matrix elements in Z-coordinates can be computed once a suitable atomic numeration of the molecular tree topology is done. But great care, especially in the case of a molecule containing an internal ring structure, should be taken in the atomic ordering for constructing a proper Z-matrix in order to keep relations like Eqs. (45)±(47) valid [24]. (Note that in this formulation even the atomic ordering given in Fig. 1 cannot be used but should be modi®ed by setting the base to one of the end atoms, e.g., atom 6.) An alternative way of recursively computing Amatrix elements in Z-coordinates can be found from Eq. (27). Direct di€erentiations of this equation w.r.t. the newly speci®ed Z-coordinates, rkd , hmd , and sld , respectively, give oxd ˆ cos hmd emk orkd
‡ sin hmd cos sld emk  ulk sin sld ulk ; …48† oxd ˆ rkd sin hmd emk ohmd ‡ rkd cos hmd cossld emk  ulk oxd ˆ osld


sin sld ulk ; …49†


rkd sin hmd sin sld emk  ulk ‡ cossld ulk : …50†

71

The derivatives of xd w.r.t. Z-coordinates that are not de®ned to the current stage of atom d are all zero. The derivatives of xd w.r.t. all the previously de®ned Rc (di€ering from rkd , hmd , and sld ) are obtained from oxd oxk ˆ oRc oRc

oemk oRc
o ‡ rkd sin hmd cos sld c emk  ulk oR oulk rkd sin hmd sin sld c ; oR rkd cos hmd

…51†

where values of oemk =oRc and oulk =oRc are found to be, respectively, oemk ˆ oRc



0; 1 rmk

for Rc ˆ
any bond stretch coordinate; oxm ; otherwise; c oR

oxk oRc

…52†

8 > < 0;

oulk ˆ > oRc :

for Rc ˆ hlk or any bond stretch  coordinate;  oeml 1 mk  emk ‡ eml  oe ; otherwise: oRc oRc sin hlk …53†

Similarly, the second order A-matrix elements in Z-coordinates can also be computed by direct di€erentiation. It should be noted that, with predetermined values of oxl =oRc , oxm =oRc , and oxk =oRc in the previous stages, the new method requires only values of rmk , rkd , hlk , hmd , and sld , which are all Z-coordinates de®ned from a local orientation of only four atoms, and is much easier to apply since it does not require as strict atomic ordering in constructing a proper Z-matrix as the method using Eqs. (45)±(47), although it requires a few more computational steps. Expansion to A-matrix elements in all primitive internals nk can be achieved by using redundancy relations, hd …n† ˆ 0 …d ˆ 1; . . . ; n f †. Analytic redundancy relations in primitive internals corresponding to some important local molecular substructures are shown in Appendix A. Each redundancy relation can be expanded to the second order by

72

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

f X bˆ1

adb DRb ‡

n f f X 1X d bdb Dnf ‡b ‡ abc DRb DRc 2 bˆ1 b;cˆ1

1X d b Dnf ‡b Dnf ‡c 2 b;cˆ1 bc

oxkk ˆ dkj ; oxjf

…60†

o2 xkk o2 xkk o2 xkk ˆ 0 ˆ ˆ j j oxif oxf oS c oxf oxjf oS c

…61†

n f

‡ ‡

f X n f X

1 cbcd DRb Dnf ‡c ˆ 0 2 bˆ1 cˆ1

…54†

with adbc  o2 hd =oRb oRc , bdbc  o2 hd =onf ‡b onf ‡c , and cbcd  o2 hd =oRb onf ‡c , where Rb is not involved in the dth redundancy relation if adb ˆ 0. If we de®ne Ud  fRb jadb 6ˆ 0g and arrange so that nf ‡d is contained in the dth redundancy relation, then the derivatives of xk w.r.t. the excluded primitive internals are found to be oxk ˆ onf ‡d 2

X oxk bd d ; oRb adb b2U

o xk ˆ onf ‡d oRc

…55†

X o2 xk bd d ; b oRc ad oR b b2U

…56†

d

X X o2 xk bcc bd o2 xk d ˆ f ‡d oRa oRb aca adb on on a2Uc b2Ud

gijtran 

p X oxk oxk mk i
j oxf oxf kˆ1

ˆ mf dij

…1 6 i; j 6 3†:

…63†

Considering the following two relations: ‰Di ; Dj Š ˆ eijk Dk ;

X oxk bdcd : oRb adb b2U d

…57† Therefore, from Eq. (15) the A-matrix elements in internal coordinates S c can be obtained by n oxd X oxd k ˆ V ; k c c oS kˆ1 on

…58†

n n X o2 x k o2 xk x c X oxk c ˆ V V ‡ V : x c b a oS a oS b c;xˆ1 on on onc ab cˆ1

…59†

3.3. Computation of A-matrix elements for external coordinates For any changes in external translations, the second term on the right-hand side of Eq. (24) remains constant, so A-matrix elements for the external translations are simply found to be

…62†

The A-matrix elements for the external rotations given by Eq. (25) can be obtained by derivatives of Eq. (24) w.r.t. external rotation parameters. Since values of xf and zkm stay constant for any external rotations at atom f, we have oxk oN ˆ zkm : o/j o/j

d

f ‡c

with S c being any generalized coordinate. The corresponding g-matrix elements are found to be

…64†

exp… /k Dk †Dj exp …/k Dk † ˆ Dj ‡ /k ‰Dj ; Dk Š ‡

1 k 2 …/ † ‰‰Dj ; Dk Š; Dk Š ‡


; 2! …65†

with eijk being the Levi±Civita density, we can derive oN…~ /† ˆ … cos/2 cos/3 D1 ‡ cos /2 sin /3 D2 o/1 sin /2 D3 †N…~ /†; oN…~ /† ˆ … cos/3 D2 o/2 oN…~ /† ˆ D3 N…~ /†: o/3

sin /3 D1 †N…~ /†;

…66† …67†

…68†

Substituting these into Eq. (63), we have with zk  Nzkm

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

oxk o~ h  zk ; j ˆ o/j o/

…69†

where ohk =o/j  Wjk are de®ned by the matrix elements of the transformation 10 1 1 0 11 0 dh sin /3 0 d/ cos/3 cos /2 @ dh2 A  @ sin /3 cos /2 cos /3 0 A@ d/2 A sin /2 dh3 0 1 d/3 0 11 d/  W…~ /†@ d/2 A; d/3 …70†

Ijk f  ˆ

p X

 mk …zk
zk †djk

kˆ1 p X

mk

kˆ1

73

zjk zkk



oxk oxk
: ohj ohk

…77†

4. Internal coordinate dynamics with B-matrix 4.1. Basic formulation

oxk ^ ˆ j  zk ; ohj

…73†

  o2 xk ^j  ^i  zk ; ˆ ohi ohj

…74†

We present an ICMD method that avoids direct inversion of the mass matrix of Eq. (31). Note that its inverse can be obtained by ®nding the inverse of the A-matrix. However, instead of inverting the A-matrix of Eq. (28), we directly use the easily computable spectroscopic B-matrix [35±37], which was originally devised to explain observed vibrations of molecules. Each atom constantly vibrates or changes its relative coordinates w.r.t. the local molecular orientation, and a molecule is considered to be entirely ¯exible. As long as a molecule is considered to be entirely ¯exible, the spectroscopic B-matrix elements for an internal coordinate, viz., the changes in the internal coordinate due to in®nitesimal changes in each atomic coordinate, are nonzero only for the atoms that are involved in de®ning the internal coordinate. We now show that the spectroscopic B-matrix computed in nonredundant coordinates of a completely ¯exible molecule is the exact inverse of the A-matrix. Since for a set of generalized coordinates S c we have

…75†

oxkk X oxkk oS c m k ˆ j ˆ dk dj ; c oS oxjm ox m c

which also explicitly de®nes the molecule's angular velocity due to the external rotation at atom f by _ _ ~ ~ X h  W…~ /†~ /. The second-order A-matrix elements for the external rotations can be obtained from the derivatives of Eq. (69): h o2 xk o2~ o~ h  i j ˆ i j  zk ‡ o/ o/ o/ o/ o/j

! o~ h  zk ; o/i

o2 xk o~ h oxk o2 xk ˆ  ˆ j j oS c o/j oS c o/ oS c o/

…71†

…72†

with S c being an internal coordinate. Thus, the corresponding elements in hj become, respectively,

2

2

o xk oxk o xk ˆ ^j  c ˆ j c : oS oS c ohj oh oS

The g-matrix elements for external rotations of ~ / are found to be gijrot 

p X kˆ1

mk

 oxk oxk  T i
j ˆ W If W ij o/ o/

…1 6 i; j 6 3† …76†

where If is the molecule's inertia tensor at atom f de®ned by

…78†

the B-matrix elements for external rotations, in combination with spectroscopic B-matrix elements for the internal coordinates, can be determined to give AB ˆ 13px3p

…79†

with the A-matrix being obtained by the methods described in Section 3. (More details will follow in Section 4.2.) If the set of generalized coordinates

74

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

S c is linearly independent, then both A and B are square nonsingular matrices and BA ˆ 13px3p :

…80†

Thus, the spectroscopic B-matrix calculated in linearly independent coordinates of an entirely ¯exible molecule is the exact inverse of the A-matrix computed by the methods described in Section 3. In this case, moreover, the spectroscopic matrix G  BuBT is the inverse of the mass-matrix g  AT u 1 A, and instead of Eqs. (41) and (42), we have 8 p P 3 P > kk > Bakk ‰A2 Šrs ; > gars ˆ > > kˆ1kˆ1 > > > p P 3 > a k P > > ˆ Bakk ‰A2 Škrj g > rj > > kˆ1kˆ1 >  > > 3 P > kk > 1 kk al > ‡ 2 mk Ar G ‰A2 Š jl ; <  lˆ1 gacd ˆ p P 3 P > kk > a 1 > Bakk ‰A2 Š‡ij gij ˆ 2 > > > kˆ1kˆ1 > >  k  3 > P > kk kk al k kk > > ‡ m G A ‰A Š ‡ A ‰A Š ; k 2 jl 2 il > i j > > lˆ1 > > p 3 > P P > kn > : gajn ˆ 12 mk Gal ‰A2 Š jl ; kˆ1

lˆ1

…81†

fa ˆ

p X 3 X 1 a oV Bmj j ; m oxm m mˆ1 jˆ1

…82†

which allows the entire process to be executed very eciently (where r and s are internal coordinate indices while i, j, and l represent external rotations and n represents an external translation). Values of all nonzero Bakk are analytically calculated without any inversion of a large matrix. If the considered generalized coordinates are subject to some redundancy relations, a spectroscopic nonsquare matrix B satisfying Eq. (79) can still be obtained with the corresponding row vectors of B not being independent. In this case Eq. (80) no longer holds, but we can construct a suitable set of nonredundant coordinates in a way similar to that used in Eqs. (36)±(42), viz., from an orthogonal matrix PT  …KT LT † that diagonalizes the spectroscopic symmetric matrix G  BuBT . Since LB ˆ 0, inserting PT P into the left-hand side

of Eq. (79), we have AB ˆ 13px3p with A  AKT and B  KB. Note that although K is a rectangular matrix, both A and B are square matrices of order 3p and are inverses of each other so that BA ˆ 13px3p . Let us de®ne G  BuBT and g  AT u 1 A, which are inverses of each other. Then, in terms of the new nonredundant coordinates, each molecule's equation of motion is now found to be of the same form as Eq. (40) with gacd and f a being also given by Eqs. (81) and (82), respectively, neglecting the underlines. Thus, the actual dynamics can eciently be executed in the transformed nonredundant coordinates. Whenever new values _ and S   KS  are determined, of S  KS, S_  KS, the corresponding values in the original redundant coordinates can be obtained by KT . It should be noted that an ICMD simulation can always be done using a suitable set of nonredundant internal coordinates, where the transformation matrix K to a new set of linearly independent coordinates need not be computed. Thus, for nonredundant generalized coordinates Eqs. (40), (81) and (82) constitute an ecient MD simulation method without any inversion of a large matrix. Moreover, if we use a set of Zcoordinates for internal degrees of freedom, extra transformations as in Eqs. (58) and (59) can be avoided. Even in the case that the V S (the bonded part of the potential) is expressed in redundant coordinates, an MD simulation can always be performed by using a set of Z-coordinates and external coordinates, since the required Cartesian derivatives oV =oxjm in Eq. (82) are easily obtained from Eqs. (7) and (10). As previously assumed, we consider only nonredundant generalized coordinates and remove all the underlines if not otherwise noted. In applying the proposed new MD formalism, nonzero elements of A and A2 are to be computed by one of the two methods described in Section 3 while those of B are to be directly computed as in vibrational spectroscopy. In addition to the advantage of avoiding any direct inversion of a large matrix, as in the recursive algorithms [26±28] in A-matrix formulations, the new ICMD formalism based on Eqs. (40), (81), and (82) has other features. Since only a few (at most four) atoms are involved in de®ning an internal coordinate, the full B and A2 need not be

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

stored if we keep related atom indices for each internal coordinate. Thus, in computing gars elements for the ath internal coordinate the actual multiplications in the summation of Eq. (81) need be done only for nonzero B-matrix elements, which signi®cantly reduces the number of required arithmetic operations compared to using Eq. (41) with multiplications for all coordinate indices b (where even the required computation of all gbcd is nontrivial). By analogy, considering only nonzero j B-matrix elements in Eq. (82) with oV =oxmc computed from Eqs. (7) and (10) gives an e€ective way of computing f a compared to using Eq. (42) with oV =oScb computed from Eq. (34). (As an example, the computation of f a for the CCCC torsional coordinate in butane, which consists of 14 atoms, can be done only for 12 backbone Cartesian components rather than all 42 Cartesian components.) We now consider a method to freeze some (nonredundant) internal coordinates in the course of a time trajectory. When an internal coordinate is constrained to a constant value, a movement of an atom directly related to this internal may result in a change in an adjacent internal coordinate that is not de®ned from this atom but from others, making the spectroscopic B-matrix physically ill de®ned in the sense that it is not the exact inverse of the A-matrix [30]. However, we may ®rst construct the systemÕs equation of motion from the spectroscopic B-matrix computed by regarding each molecule to be completely ¯exible, and next introduce the desired constraints in the equation, which simply become S r …t† ˆ S r …0† ˆ constant

…83†

with r being the indices of internals to be frozen. The resulting systemÕs equation of motion is found to be 3p X cˆ1

gbc Sc ‡

3p X

gbcd S_ c S_ d ˆ

c;dˆ1

X s oV ‡ ks db ; oS b s

…84†

where each additional parameter ks (Lagrange undetermined multiplier) represents the required constraint force to maintain the imposed constraint [38]. When each coordinate to be frozen is nonredundant, its related B-matrix elements are

75

still nonzero while its time derivatives are all zero. Thus, Eq. (84) can be solved to give X X guvw S_ v S_ w ˆ f u ‡ Gus ks ; …85† Su ‡ v;w

s

X X grvw S_ v S_ w ˆ f r ‡ Grs ks ; v;w

…86†

s

where u, v, and w are indices for unconstrained coordinates, while r and s are those for constrained ones. With grvw and f r given by Eqs. (81) and (82), respectively, the parameters, ks , can be determined by solving Eq. (86). These can be substituted into Eq. (85) in order to determine new accelerations for unconstrained coordinates. Since this method requires solving a linear equation, it may be less e€ective in reducing operational dimensions than the A-matrix formalism, in which the concerned coordinates can be simply neglected since the systemÕs kinetic energy is directly determined from the A-matrix [30]. Alternatively, since it is more realistic to regard each internal coordinate as ¯exible in a real molecule and the spectroscopic B-matrix is still a good approximation to the inverse of the A-matrix, as a useful approximation we may directly use the spectroscopic B-matrix elements in order to ®nd an initial equilibrium con®guration in the early stages of MD simulations by simply neglecting the fast varying internals without solving the constrained equation as in the A-matrix formalism. However, near the equilibrium con®guration we can eciently do full B-matrix ICMD without any rigid constraints in order to derive desired physical quantities. We note that it is ecient to compute also the related nonzero spectroscopic B-matrix elements in Z-coordinates after de®ning new atomic coordinates from Eq. (27) at each Z-matrix stage. B-matrix elements in a set of redundant internal coordinates can alternatively be computed from these. Once all B-matrix elements in Z-coordinates are determined, those in the primitive internals can be obtained from Eq. (54) by 8 b oR > < oxk ; for 1 6 b 6 f ; b on k P adc oRc ˆ …87† k ; for b ˆ f ‡ d; oxk > : bd oxk c2Ud

d

k

76

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

8 o2 Rb ; for 1 6 b 6 f ; j > > < oxm oxkk d o2 nb P ac o2 Rc P adac oRa ˆ j k j d k b bd oxj ox ox > oxm oxk > : c2Ud d m k a;c d m for b ˆ f ‡ d:

oRc oxkk

;

…88†

Using Eq. (14), B-matrix elements in (redundant) internal coordinates S a can be alternatively obtained by n c oS a X a on ˆ U ; …89† c oxkk oxkk cˆ1 n n b c 2 c X X o2 S a a o n a on on ˆ U ‡ U : c bc j j j k oxm oxkk oxm oxkk b;cˆ1 oxm oxk cˆ1

…90†

4.2. B-matrix elements for external coordinates The coordinates of atom f are translational parameters between the LFF and the MFF, and the B-matrix elements for the external translations are simply found to be oxkf oxjm

ˆ dmf dkj ;

o2 xkf

oxil oxjm

p X oS c zk  ˆ0 oxk kˆ1

If~ h mj ˆ

oxk mk zk  c ˆ 0; oS kˆ1

…94†

mk

oxjk m …dk o~ h

dmf †

f p X X oxk oxk mk
o~ h oS c cˆ1 kˆ1 p X kˆ1

…92†

Thus, this term is guaranteed to be small when the term in parenthesis is small, viz.,

p X kˆ1

ˆ ˆ 0:

…95†

derived from the invariance of all internal coordinates under external rotations, viz., oS c =o/j ˆ 0 …j ˆ 1; 2; 3†. It should be noted that the A-matrix elements oxkk =oS c for internal coordinates do not satisfy Eq. (94) in general, while the B-matrix elements for external rotations based on the Eckart condition of Eq. (94) do not satisfy Eqs. (79) and (80) with the A-matrix. Without any conditions, the correct B-matrix elements in ~ h are found to satisfy

…91†

Since internal vibrations are the main concerns in conventional spectroscopy, B-matrix elements for external rotations are usually determined under the Eckart condition [45,46], which is imposed to make the interaction terms between internal vibrations and external rotations small enough to be neglected in a molecule's kinetic energy. Using Eq. (69), the interaction term between an external rotation of /_ j and an internal motion S_ c in the kinetic energy expression of Eq. (30) is found to be ! p p X X oxk oxk o~ h oxk mk j
c ˆ j
mk zk  c : …93† oS o/ o/ oS kˆ1 kˆ1

p X

which is the Eckart condition for external rotations at atom f. This should not be confused with the relation

mk …zk  ^j†…dmk

!

oS c oxjm

dmf †

f p X X oxk m k zk  c oS cˆ1 kˆ1

!

oS c ; oxjm

…96†

where S c represent only internal coordinates and ~ hmj  o~ h=oxjm . Values of ~ hmj can be obtained by solving Eq. (96). The second derivatives ~ h=oxil oxjm can be obtained from the dehli mj  o2~ rivative of Eq. (96) w.r.t. xil . However, since these contain A2 elements for internal coordinates (due to the second term in the right-hand side of Eq. (96)), they can alternatively be computed directly by using the relation ~ hli mj ˆ

3p p X 3 X X o2 x k ~ hk k a k b oS oS a;bˆ1 kˆ1 kˆ1

!

oS a oS b oxil oxjm

…97†

obtained from a Cartesian derivative of Eq. (80). From the inverse relation of Eq. (70), the B-matrix elements for external rotations of ~ /, viz.,

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

~ /mj  o~ /=oxil oxjm , can be de/=oxjm and ~ /li mj  o2~ termined by ~ /mj ˆ W 1 …~ /†~ h mj ;

~ /li mj

…98†

( ˆ W …~ /† ~ hli mj 1

8 > < 1 ~ ˆ W …/† ~ hli mj > :

B /† ˆ @ W 1 …~

cos/3

e2 ‡ ie e0 ‡ ie3 2 1 …e ie † e0 ie3   1 0 ˆ e0 ‡ ie
r 0 1

Qˆ

1 cos /2 3

sin / sin /2 cos/3 cos /2

sin /3

1 cos /2 3

cos / sin /2 sin /3 cos /2

0

1

C 0 A: 1

Note that W …~ /† of Eq. (100) is not well-de®ned when /2 ˆ p=2. This is also a problem in the A-matrix formulation, where grot of Eq. (76) is singular and we cannot solve Eq. (32) properly. This kind of singularity may happen frequently during MM or MD simulations of a system containing many molecules. Thus, for a system of many molecules it is much better to use the Euler representation in SU…2† for the external rotations [38,47,29].

A vector vm ˆ v1m ^ 1m ‡ v2m ^ 2m ‡ v3m ^ 3m in the MFF can be represented as vm ˆ

e0 ˆ cos

e1 ˆ cos

e2 ˆ cos

 …101†

1 2 in complex two dimensions with v‡ m  vm ‡ ivm and 1 2 vm  vm ivm . Thus, a rotation N of vm with v ˆ Nvm in f^ 1; ^ 2; ^ 3g can be regarded as a representation Q with v  Qy vm Q in SU(2) such that Q is parameterized as

…103†

In particular, for the rotation given by Eq. (25) the corresponding new parameters …e0 ; e†, which are di€erent from the ones derived from the Euler angles [38,47,29], are found to be

4.3. Euler representation of rotations in SU(2)

vm v3m

…99†

det…Q† ˆ e0 e0 ‡ e
e ˆ 1:

1

v3m v‡ m

…102†

where rj …j ˆ 1; 2; 3† are the Pauli spin matrices and the vector e determines the direction of the rotation axis, which satis®es

…100†



77

 1

) oW…~ /† ~ /mj oxil 0 2 1 19 /li /mj cos /3 sin /2 ‡ /3li …/1mj sin /3 cos/2 ‡ /2mj cos/3 † > = B C ‡ @ /2li /1mj sin /3 sin /2 /3li …/1mj cos /3 cos /2 /2mj sin /3 † A ; > ; /2li /1mj cos /2

where W 1 …~ /† is found to be 0



e3 ˆ cos

/3 /2 /1 cos cos 2 2 2

sin

/3 /2 /1 sin sin ; 2 2 2 …104†

/3 /2 /1 /3 /2 /1 cos sin ‡ sin sin cos ; 2 2 2 2 2 2 …105† /3 /2 /1 sin cos 2 2 2

sin

/3 /2 /1 cos sin ; 2 2 2 …106†

/3 /2 /1 /3 /2 /1 sin sin ‡ sin cos cos ; 2 2 2 2 2 2 …107†

and the rotation matrix N in real three-dimension in terms of …e0 ; e† is found to be

78

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

0

e0 e0 e
e ‡ 2e1 e1 0 @ N…e ; e† ˆ 2e1 e2 ‡ 2e0 e3 2e1 e3 2e0 e2

2e1 e2 2e0 e3 0 0 ee e
e ‡ 2e2 e2 2 3 2e e ‡ 2e0 e1



Considering the following two relations: oN N oe0 oN N oej

1

1

0

ˆ 2‰e 1

e
DŠ;

ˆ 2‰ej 1 ‡ e0 Dj

…109† j

…e  D† Š;

1 2e1 e3 ‡ 2e0 e2 2e2 e3 2e0 e1 A: 0 0 ee e
e ‡ 2e3 e3

…110†



e0mj e mj

e0li mj eli mj

the time derivative of Eq. (24) in terms of …e0 ; e† becomes

with

~  zk ‡ N…e0 ; e†_zkm ; x_ k ˆ x_ f ‡ X

and

…111†

_ ~ ~ h, where the molecule's angular momentum, X is de®ned by !   0 e_ 0 0 …112† ~ ˆ W…e ; e† e_ X with

0

e0 B e1 W…e0 ; e†  2B @ e2 e3

e1 e0 e3 e2

e2 e3 e0 e1

1 e3 e2 C C: e1 A e0

…113†

The A-matrix elements in …e0 ; e† are found to be: oxk o~ h ˆ m  zk ; m oe oe o2 xk o ˆ oel oem oel

o~ h oem

…114† !

o~ h  zk ‡ m  oe

! o~ h  zk : oel …115†

The ®rst row of Eq. (112) corresponds just to the time derivative of Eq. (103). It is evident that, for all …e0 ; e† satisfying Eq. (103), W…e0 ; e† has the well-de®ned unique inverse: 1 W 1 …e0 ; e† ˆ WT …e0 ; e†: …116† 4 Once values of ~ hmj and ~ hli mj are calculated from Eqs. (96) and (97), the B-matrix elements in …e0 ; e† can be obtained from the inverse relation of Eq. (112):



 ˆW

1

 0 ; ~ h mj

 ˆW

1

…108†

~ hli mj

…117† 1~ h i
~ h mj 2 l 1~ i h ~ h mj 2 l

! …118†

e0mj  oe0 =oxjm ; emj  oe=oxjm ; e0li mj  o2 e0 =oxil oxjm ; eli mj  o2 e=oxil oxjm : Following Eqs. (104)±(107), the case of ~ /ˆ0 corresponds to that of …e0 ; e† ˆ …1; 0†, and in the limit of …e0 ; e† ! …1; 0† the direction of e_ coincides ~ with that of X. Although this treatment for external rotations not only removes the singularity problem involved in Eq. (100) but also avoids the computations of trigonometric functions inherent in using Eq. (98), we have to solve the problem arising from the constraint of Eq. (103) for each molecule. This can be solved by a treatment similar to that used in the case of dependent coordinates. Considering that grot in …e0 ; e† is found to be   0 0 grot ˆ WT W …119† 0 If and there exists an orthogonal matrix Q such that QIf QT ˆ D

…120†

with D being diagonal, it is evident that grot is diagonalized by    
1 rot 1 0 0 T 1 1 g W W 0 QT 0 QT   0 0 ˆ : …121† 0 D Therefore, a desired transformation matrix K to new nonredundant rotation coordinates e can be taken as

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

K  …0 Q†W 0 1 e ˆ 2Q@ e2 e3

0

e e3 e2

3

e e0 e1

2

1

e e1 A e0

…122†

where Q has to be constructed from eigenvectors of nonzero eigenvalues of If . The time derivative of e becomes ! e_ 0 _ ~ h ˆ QX; …123† e_  K ˆ Q~ e_ which is the angular velocity in the frame of principal axes of If , while the new rotation coordinates are the null coordinates since  0 e eK ˆ 0: …124† e For nonlinear molecules, it is more ecient to h, by QT : transform e to new rotation parameters, ~ _  ~ h  QT e_ , and ~ h  QTe. Thus, the transh  QT e, ~ formation matrix from Euler parameters …e0 ; e† to ~ h is found to be 0 1 1 e e0 e3 e2 N  Q T K ˆ 2@ e 2 e 3 …125† e0 e 1 A; 3 2 1 e e e e0 which is the submatrix taken from the second row to the fourth row of W…e0 ; e†. Then, the corresponding B- and A-matrix elements become, respectively,  0  e ~ h mj  N mj ˆ ~ …126† h mj ; e mj  0  e ~ hli mj  N li mj ˆ ~ h l i mj eli mj

79

which are computed from Eqs. (73)±(75) and (96)± (97) with S c being internal coordinates. Finally, we brie¯y summarize how the above results can be used in applying the proposed new ICMD formalism to model systems. Initial Euler parameters …e0 ; e† specify the relative rotation matrix N between the MFF and LFF by Eq. (108), while those of xf determine the translations between the two frames. The dynamics equation, Eq. (40), with Eqs. (81) and (82) is most eciently executed in a nonredundant set of Z-coordinates, xf (translations), and ~ h (rotations), for a nonlinear molecule. From the given initial values of Zcoordinates, we generate each molecule's atomic coordinates xk in the LFF by using Eqs. (24) and (27). At each stage of atomic coordinate generation we also compute nonzero elements of A, A2 , and B in Z-coordinates directly. (The method of computing spectroscopic B-matrix elements can be found elsewhere [35±37].) Next, we compute Aand B-matrix elements for external coordinates. Those for xf are obtained from Eqs. (60), (61), and h are obtained from Eqs. (91), while those for ~ (128)±(130) and (126) with Eqs. (73)±(75) and (96) being used. The time trajectory of …e0 ; e† is obtained in the following way. Initial values of …_e0 ; e_ † _ h by N, and from Eq. are transformed to those of ~  (40) new values of ~ h are computed and transformed back to values of …e0 ; e† by NT =4. These are then used to determine new values of …e0 ; e† by Eq. (43). 5. Multiple time step algorithm

…127†

With Jcc being the conjugate momenta of the generalized coordinates Scc as in Eq. (35), the systemÕs Hamiltonian is expressed by

3 X 3 oxk X oxk ohk oel oxk ˆ ˆ ; k oel ohj ohj ohj kˆ1 lˆ0 oh

…128†

H …S; J† ˆ

o2 xk o2 xk ; i j ˆ oh oh ohi ohj

…129†

which provides for the systemÕs Hamiltonian equations of motion

o2 xk o2 xk o2 xk ; j ˆ j ˆ c c oS oh oS oh ohj oS c

…130†

oH J_ ca ˆ a ˆ oSc

1~ hli  ~ h mj ; 2

mol 1X JT g 1 Jc ‡ V …S† 2 cˆ1 c c

N

1 _ T ogc _ oV Sc a Sc ‡ a ; 2 oSc oSc

…131†

…132†

80

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

a oH  S_ ca ˆ a ˆ gc 1 Jc : oJc

…133†

While Eq. (133) is equivalent to the de®nition of Jcc , Eq. (132) is equivalent to Eq. (32), the systemÕs Euler±Lagrange equation of motion. Let us concentrate on the cth molecule of the system and remove the molecular index. Then, the molecule's Liouville operator L can be de®ned by  X o o iL  S_ b b ‡ J_ b b ; …134† oS oJ b which is Hermitian and propagates the initial values of fS…0†; J…0†g to those of fS…t†; J…t†g ˆ exp…iLt†fS…0†; J…0†g

…135†

at time t. Note that in the generalized coordinates S c the fast varying degrees of freedom, e.g., bond stretches or angle bend coordinates, can be explicitly separated from the slowly varying degrees of freedom such as torsions, external translations, or external rotations. So, we can decompose the Liouville operator into two parts for simplicity: iL ˆ iLs ‡ iLf  X o r r o _ _  S ‡J oS r oJ r r   X o o ‡ S_ u u ‡ J_ u u ; oS oJ u

…136†

where r and u are coordinate indices for slowly varying and fast varying coordinates, respectively. This enables us to write, to the order of O…Dt3 † [12],     Dt Dt exp …iLDt†  exp iLs exp …iLf Dt† exp iLs 2 2   Dt n ˆ exp iLs f exp …iLf dt†g 2   Dt  exp iLs …137† 2 where the inner fast varying part is further resolved by a smaller time step of dt  Dt=n (n can be chosen to give an optimized process). The state at time Dt is determined by applying this propagator to the initial state of fS…0†; J…0†; 0g:

(a) exp …iLs Dt=2†fS…0†; J…0†; 0g provides the state of fS r …Dt=2†; J r …Dt=2†; Dt=2g and fS u …0†; J u …0†; Dt=2g, where S r …Dt=2† and J r …Dt=2† are found to be, respectively,   Dt r Dt S …138† ˆ S r …0† ‡ S_ r …0†; 2 2  J

r

Dt 2

 ˆ J r …0† ‡

Dt _ r J …0†: 2

…139†

Instead of Eq. (139), with Eq. (40) we can equivalently use   _S r Dt ˆ S_ r …0† ‡ Dt Sr …0†: …140† 2 2 (b) exp …iLf dt†fS…Dt=2†; J…Dt=2†; Dt=2g provides the state of fS r …Dt=2†; J r …Dt=2†; Dt=2 ‡ dtg and fS u …dt†; J r …dt†; Dt=2 ‡ dtg. Applying this operation a total of n times, which is equivalent to using the velocity Verlet integrator iteratively n times on the system with only fast varying degrees of freedom, generates the state of fS r …Dt=2†; J r …Dt=2†; Dt=2g and fS u …Dt†; J u …Dt†; Dtg. (c) Applying exp…iLs Dt=2† to the ®nal state of step (b) provides the state of fS…Dt†; J…Dt†; Dtg. Thus, we may use MTS depending on the nature of the degrees of freedom: a longer time step for slowly varying degrees of freedom while a smaller time step for fast varying degrees of freedom. This is di€erent from the conventional MTS algorithm [10±13], where the Liouville operator is decomposed depending on the classi®cation of interaction distances, since the forces acting on an atom due to atoms located at far distances stay roughly constant while those arising from near atoms change fast. In fact, each molecule changes its position and the classes of interaction distances may also change in time, which reduces the e€ectiveness of the conventional method. It should be mentioned that, if we follow the spectroscopic Bmatrix formulation, at each time step, whether it is for fast or slow moving degrees of freedom, we have to solve a linear equation to determine the Lagrange undetermined multipliers introduced to freeze internal degrees of freedom. However, as a way to obtain an initial rough equilibrium con®guration, we may also explore the trajectories

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

that result from simply neglecting the frozen coordinates, as in the A-matrix formalism. Such con®gurations could then be used for further unconstrained dynamics. 6. Concluding remarks Based on the spectroscopic B-matrix in nonredundant generalized coordinates, we propose a new ICMD formalism that does not require any matrix inversion except for the 3  3 inertia tensor, which makes it useful for MD simulations on systems of large molecules. The inversion of the mass-matrix is avoided by directly computing the desired spectroscopic B-matrix that gives the exact inverse of the conventional A-matrix. It turns out that the B-matrix elements for external rotations under the Eckart conditions [45,46] do not provide such an exact inverse of the A-matrix. Instead, the correct B-matrix elements for external rotations are obtained (by solving Eq. (96)) without any conditions. The sparsity in the spectroscopic Bmatrix enhances the eciency of the proposed formalism by signi®cantly reducing the number of arithmetic operations required. This approach is ideal for MD simulations of many-molecule systems in combination with the proposed dynamics method of singularity free external rotations in a representation of SU…2† algebra [38]. Based on the explicit separability between fast varying degrees of freedom and slowly varying ones, we introduce an MTS method in generalized coordinates in which the nontrivial classi®cation of interaction distances [10±12] is not required. An MD simulation in redundant internal coordinates can be done by ®nding a transformation to a new set of independent coordinates, in which the actual dynamics equation is solved. Since this may involve nontrivial computations of eigenvectors of the spectroscopic matrix G, the proposed B-matrix formalism is then not optimal. In every case, however, we can easily ®nd a suitable set of nonredundant internal coordinates for MD simulations. In particular, if we use a proper set of Zcoordinates, which are linearly independent and orthogonal to each other, the transformations to other (redundant or nonredundant) internals are

81

not necessary, giving an ecient set of internal coordinates for the proposed B-matrix formalism. Even in the case that the bonded part of the potential energy function contains terms expressed in redundant internal coordinates, its Cartesian derivatives are easily computed and substituted into Eq. (82) so that the actual dynamics simulation can be performed in Z-coordinates. Although the spectroscopic B-matrix is not well-de®ned for rigid constraints on internal coordinates [30], it is still a good approximation to the inverse of the A-matrix for real macromolecules. Thus, as an e€ective way of ®nding an approximate initial equilibrium con®guration before any full rigorous MD process, this approach can be used by simply neglecting the constrained coordinates to reduce the operational dimension, as in the A-matrix ICMD formalism. This may also be a useful way of ®nding a folding path from an arbitrary extended conformation of a macromolecule like a protein or high polymer to its stable (native) structure in a solvent. The proposed B-matrix ICMD formalism is under implementation in our laboratory and it will be applied to realistic model systems and compared with other MD simulation methods.

Acknowledgements We are indebted to Weili Qian for many helpful discussions. This research was supported by NSF grants MCB-9903991 and DMR-9902727. Additional support was provided by the Air Force Research Laboratory/Materials & Manufacturing Directorate and by the Common High Performance Software Support Initiative of the Department of Defense High Performance Computing Program. Appendix A. Redundancy relations in primitive internals Most redundancy relations in primitive internals of a molecule can be derived from a basic equation that is related to a local structure of ®ve atoms. This local structure can be best represented by the one from atom 1 to atom 5 in Fig. 1. Using

82

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

Eq. (27), from the atomic ordering for this local structure we have

Substitution of Eq. (A.8) into Eq. (A.5) provides another redundancy relation

e34 ˆ

cos…s79

cosh24 e23 ‡ sin h24 … coss14 e23  u13

e35 ˆ

cosh25 e23 ‡ sin h25 … coss15 e23  u13

sin s14 u13 †; …A:1†

s210 †:

…A:9†

Similar considerations to cosh810 and cosh910 lead to, respectively, sin s15 u13 †: …A:2†

Since u13
…e23  u13 † ˆ 0, we obtain a relation between the excluded (from Z-coordinates) bond angle h45 … cos h45  e34
e35 † and Z-coordinates h24 , h25 , s14 , and s15 [48]: cos h45 ˆ cosh24 cosh25 ‡ sin h24 sin h25 cos…s14

s15 †:

…A:3†

When s14 s15 ˆ p, viz., a planar local structure, h45 is independent of s14 and s15 : h45 ‡ h24 ‡ h25 ˆ 2p:

…A:4†

Redundancy relations for a locally tetrahedral structure can be easily derived by applying Eq. (A.3). As an illustration, we consider the left-end local structure in Fig. 1. Assume that the coordinates of atoms 9, 10, and 11 are de®ned starting not from atom 7 but from atom 2, viz., the primitive internals h19 , h110 , h111 , s29 , s210 , and s211 belong to the Z-coordinates, while h910 , h911 , h1011 , s79 , s710 , and s711 do not. The redundancy relations for h910 , h911 , and h1011 in terms of Z-coordinates are found to be, respectively,

cos…s79 cos…s710

s711 † ˆ cos …s29 s711 † ˆ cos …s210

s211 †; s211 †:

s210 †;

…A:5†

‡ sin h912 sin h92 cos…s712

…A:6†

‡ sin h13 sin h113 cos…s15 s211 †: …A:7†

Note that cosh910 also satis®es cos h910 ˆ cosh19 cosh110 s710 †:

…A:8†

s121 †;

…A:13†

s93 †;

…A:14†

s115 †;

…A:15†

cosh111 ˆ cos h13 cos h113

cos h1011 ˆ cos h110 cos h111 ‡ sin h110 sin h111 cos…s210

…A:12†

cosh13 ˆ cos h111 cos h113 ‡ sin h111 sin h113 cos…s91

s211 †;

s72 †;

cosh912 ˆ cos h92 cos h122

cos h911 ˆ cosh19 cosh111 ‡ sin h19 sin h111 cos …s29

…A:11†

cosh122 ˆ cos h912 cosh92

‡ sin h92 sin h122 cos…s91

‡ sin h19 sin h110 cos …s29

…A:10†

Next, we consider the redundancy relations in benzene, which is a basic cyclic structure. Each atom is numbered as in Fig. 2, and the unnecessary bond length and bond angle internal coordinates for the Z-matrix construction are designated by dashed lines while those of Z-coordinates are shown by solid lines. Related to the excluded bond vector x112 , thirteen primitive internals excluded from the Z-matrix construction are listed to be not only h122 , h111 , s91 , s102 , s114 , s121 , and s123 , but also the six coordinates r112 , h92 , h113 , s72 , s93 , and s115 coming from the inside hexagonal carbon ring. Ten analytic redundancy relations can be obtained in ways similar to the above:

cos h910 ˆ cosh19 cosh110

‡ sin h19 sin h10 cos…s79

s710 † ˆ cos …s29

cos…s72

s712 † ˆ cos …s102

s1012 †;

…A:16†

cos…s72

s102 † ˆ cos …s712

s1012 †;

…A:17†

cos…s93

s91 † ˆ cos…s123

s121 †;

…A:18†

cos…s93

s123 † ˆ cos …s91

s121 †;

…A:19†

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

83

Fig. 2. Schematic drawing of the benzene molecule with atomic numbering consistent with a Z-matrix construction.

cos …s115

s114 † ˆ cos…s15

s14 †;

…A:20†

cos …s115

s15 † ˆ cos…s114

s14 †:

…A:21†

The remaining three redundancy relations can be obtained from a closure condition of the internal carbon ring: P  x23 ‡ x35 ‡ x57 ‡ x79 ‡ x911 ‡ x112 ˆ 0: …A:22† For example, redundancy relations belonging to planar symmetry species A1g , B2u , and E2g can be found, respectively, from [36] P
…x23 ‡ x35 ‡ x57 ‡ x79 ‡ x911 ‡ x112 † ˆ 0; …A:23† P
…x23

x35 ‡ x57

x79 ‡ x911

x112 † ˆ 0; …A:24†

P
…x23

2x35 ‡ x57 ‡ x79

2x911 ‡ x112 † ˆ 0: …A:25†

In order to compute all the dot products in Eqs. (A.23)±(A.25), which will provide relations in terms of primitive internals, we need relations similar to 0

1 0 1 e57 e35 @ e57  u37 A ˆ T35 …h25 ; s14 †@ e35  u25 A ; u37 u25 5 3 …A:26† where fe35 ; e35  u25 ; u25 g can be considered as the basis of a local coordinate system at atom 3 and the transformation matrix T35 …h25 ; s14 †, which can be obtained from local geometric relations and Eq. (27), is found to be

T35 …h25 ; s14 † ˆ exp … i‰p h24 ŠD3 † exp … is14 D1 † 0 1 cos h24 sin h24 cos s14 sin h24 sin s14 ˆ @ sin h24 cos h24 coss14 cos h24 sin s14 A: 0 sin s14 coss14

…A:27†

84

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

From the ring closureness, we should have T35 T57 T79 T911 T112 T23 ˆ 1;

…A:28†

which contains no bond length dependency. The underlying essence in Eqs. (A.22) and (A.28) can be applied to any molecular ring structure. It should be noted that Eq. (A.28) provides three additional redundancy relations. For simplicity, we can choose …T35 T57 T79 T911 T112 T23 †^ 1ˆ^ 1;

…A:29†

since bringing the vector e35 to the initial state through cyclic symmetry transformations is equivalent to bringing the other two vectors e35  u25 and u25 to the starting position. For the case of the benzene molecule, where six hydrogen atoms are sticking out from the internal carbon ring, the three redundancy relations from Eq. (A.28) are not necessary, since Eqs. (A.12)±(A.21) provide all the requirements. Appendix B. External coordinates in the center of mass frame From a molecule's atomic coordinates xk in f^ 1; ^ 2; ^ 3g, the center of mass (c.m.) position, c, is de®ned by 1X mk xk : M kˆ1 p

c

…B:1†

If we de®ne yk as the distance vector from the c.m. to the kth atom, then the equivalent to Eq. (24) is expressed by xk ˆ c ‡ yk

…B:2†

with yk  Nykm , where N represents the relative rotation between the LFF and the c.m. frame (CMF). The de®nition of Eq. (B.1) in combination with Eq. (B.2) provides p X kˆ1

mk yk ˆ 0:

…B:3†

If we regard the components of c as translational parameters, the corresponding A- and Bmatrix elements are found to be

oxkk ˆ dkj ; ocj

…B:4†

o2 xkk o2 xkk o2 xkk ˆ 0 ˆ ˆ ; oci ocj oS c ocj ocj oS c

…B:5†

oc mm ^ j; j ˆ oxm M

…B:6†

o2 c ˆ0 oxil oxjm

…B:7†

with S c being internal coordinates. Since for rotations at the c.m., c and the internal coordinates S c remain constant, the A-matrix elements in ~ h are found to be oxk ^ ˆ j  yk ; ohj

…B:8†

  o2 x k ^ ^ i j ˆ j  i  yk ; oh oh

…B:9†

2 o2 xk ^j  oxk ˆ o xk : ˆ j oS c ohj oS c oS c oh

…B:10†

It should be noted that, in general for ¯exible molecules, a slight change in xk or in internal coordinates S c will result in changes in the c.m. position, giving 3 ocj XX ocj oxkk ˆ 6 0: ˆ oS c oxkk oS c kˆ1 kˆ1 p

…B:11†

Thus, Eq. (80) no longer holds. However, if we take the reference frame as the CMF making c ˆ 0, the corresponding A and B satisfy the required Eqs. (79) and (80) for B-matrix ICMD simulations. Since oykk =oxjm ˆ …dmk mm =M†dkj and oxjm =oykk ˆ djk , the desired elements of B and A for internal coordinates S c in CMF are obtained by, respectively, oS c oS c ˆ ; oymj oxjm oymj oxjm ˆ oS c oS c

…B:12† 1 X oxjk mk : M kˆ1 oS c p

…B:13†

S.-H. Lee et al. / Chemical Physics 265 (2001) 63±85

Eq. (B.13) holds also for external coordinates with oxjm =oS c being given by Eqs. (B.4) and (B.8). The elements of A2 in CMF are obtained from derivatives of Eq. (B.13). Thus, the desired elements of B for external rotations, ~ hmj  o~ h=oymj , can be obtained by solving the relation equivalent to Eq. (96): ! f p j X X oy oy oy oS c k k Icm~ h mj ˆ m m m mk
c oymj o~ h o~ h oS cˆ1 kˆ1 ! f p X X oyk oS c ^ ˆ m m ym  j m k yk  c oS oymj cˆ1 kˆ1 …B:14† with Ijk cm , the molecule's inertia tensor w.r.t. its c.m., being given by Ijk cm  ˆ

p X kˆ1 p X

 mk …yk
yk †djk mk

kˆ1

oyk oyk
: ohj ohk

yjk ykk



…B:15†

Note that the quantity in parenthesis on the righthand side of Eq. (B.14) is nonzero in general. References [1] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [2] J.A. McCammon, S.C. Harvey, Dynamics of Proteins and Nucleic Acids, Cambridge University Press, Cambridge, 1987. [3] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [4] H.C. Andersen, J. Comput. Phys. 52 (1983) 24. [5] J.P. Ryckaert, J. Comput. Phys. 55 (1985) 549. [6] R. Edberg, D.J. Evans, G.P. Morriss, J. Chem. Phys. 84 (1986) 6933. [7] D.J. Tobias, C.L. Brooks III, J. Chem. Phys. 89 (1988) 5115. [8] E. Barth, K. Kuczera, B. Leimkuhler, R.D. Skeel, J. Comput. Chem. 16 (1995) 1192. [9] R. Kutteh, J. Chem. Phys. 111 (1999) 1394. [10] W.B. Streett, D.J. Tildesley, G. Saville, Mol. Phys. 35 (1978) 639. [11] H. Grubm uller, H. Heller, A. Windemuth, K. Schulten, Mol. Simul. 6 (1991) 121. [12] M. Tuckerman, B.J. Berne, G.J. Martyna, J. Chem. Phys. 97 (1992) 1990.

85

[13] D.D. Humphreys, R.A. Friesner, B.J. Berne, J. Phys. Chem. 98 (1994) 6885. [14] J.J. Biesiadecki, R.D. Skeel, J. Comput. Phys. 109 (1993) 318. [15] J.A. Izaguirre, S. Reich, R.D. Skeel, J. Chem. Phys. 110 (1999) 9853. [16] C.S. Peskin, T. Schlick, Commun. Pure Appl. Math. 42 (1989) 1001. [17] G. Zhang, T. Schlick, J. Chem. Phys. 101 (1994) 4995. [18] B. Paizs, G. Fogarasi, P. Pulay, J. Chem. Phys. 109 (1998) 6571. [19] J. Baker, D. Kinghorn, P. Pulay, J. Chem. Phys. 110 (1999) 4986. [20] W.L. Jorgensen, J. Tirado-Rives, J. Phys. Chem. 100 (1996) 14508. [21] N. G o, H.A. Scheraga, Macromolecules 3 (1970) 178. [22] T. Noguti, N. G o, J. Phys. Soc. Jpn. 52 (1983) 3283. [23] A.K. Mazur, R.A. Abagyan, J. Biomol. Struct. Dyn. 6 (1989) 815. [24] A.K. Mazur, R.A. Abagyan, J. Biomol. Struct. Dyn. 6 (1989) 833. [25] K.D. Gibson, H.A. Scheraga, J. Comput. Chem. 11 (1990) 468. [26] A. Jain, N. Vaidehi, G. Rodriguez, J. Comput. Phys. 106 (1993) 258. [27] L.M. Rice, A.T. Br unger, Proteins 19 (1994) 277. [28] J. Turner, P. Weiner, B. Robson, R. Venugopal, H. Schubele III, R. Singh, J. Comput. Chem. 16 (1995) 1271. [29] G.R. Kneller, K. Hinsen, Phys. Rev. E 50 (1994) 1559. [30] S. He, H.A. Scheraga, J. Chem. Phys. 108 (1998) 271. [31] A.K. Mazur, J. Chem. Phys. 111 (1999) 1407. [32] A. Amadei, A.B.M. Linssen, H.J.C. Berendsen, Proteins 17 (1993) 412. [33] D.-S. Bae, E.J. Haug, Mech. Struct. Mach. 15 (3) (1987) 359. [34] G. Rodriguez, A. Jain, K. Kreutz-Delgado, J. Astronaut. Sci. 40 (1992) 27. [35] E.B. Wilson, J.C. Decius, P.C. Cross, Molecular Vibrations, McGraw-Hill, New York, 1955. [36] S. Califano, Vibrational States, Wiley, New York, 1976. [37] S.-H. Lee, K. Palmo, S. Krimm, J. Comput. Chem. 20 (1999) 1067. [38] H. Goldstein, Classical Mechanics, Addison-Wesley, Reading, MA, 1971. [39] H.J.C. Berendsen, Science 282 (1998) 642. [40] Y. Duan, P.A. Kollman, Science 282 (1998) 740. [41] K. Palmo, L.-O. Pietila, S. Krimm, J. Comput. Chem. 12 (1991) 385. [42] K. Palmo, L.-O. Pietila, S. Krimm, Comput. Chem. 17 (1993) 67. [43] K. Palmo, N.G. Mirkin, S. Krimm, J. Phys. Chem. A 102 (1998) 6448. [44] W.D. Allen, A.G. Csaszar, J. Chem. Phys. 98 (1993) 2983. [45] C. Eckart, Phys. Rev. 47 (1935) 552. [46] A. Sayvetz, J. Chem. Phys. 7 (1939) 383. [47] D.J. Evans, Mol. Phys. 34 (1977) 317. [48] T. Schlick, J. Comput. Chem. 9 (1988) 861.