Molecular dynamics integration and molecular

0 downloads 0 Views 188KB Size Report
sufficiently small that we can use the dynamical molecular model. 6 ..... small enough to accurately describe the high-frequency mo- lecular vibrations. In this way ...
THE JOURNAL OF CHEMICAL PHYSICS 122, 174102 共2005兲

Molecular dynamics integration and molecular vibrational theory. II. Simulation of nonlinear molecules Matej Praprotnik and Dušanka Janežiča兲 National Institute of Chemistry, Hajdrihova 19, 1000 Ljubljana, Slovenia

共Received 17 November 2004; accepted 11 February 2005; published online 29 April 2005兲 A series of molecular dynamics 共MD兲 simulations of nonlinear molecules has been performed to test the efficiency of newly introduced semianalytical second-order symplectic time-reversible MD integrators that combine MD and the standard theory of molecular vibrations. The simulation results indicate that for the same level of accuracy, the new algorithms allow significantly longer integration time steps than the standard second-order symplectic leap-frog Verlet method. Since the computation cost per integration step using new MD integrators with longer time steps is approximately the same as for the standard method, a significant speed-up in MD simulation is achieved. © 2005 American Institute of Physics. 关DOI: 10.1063/1.1884608兴 d␩ = 兵␩,H其 = LˆH␩ , dt

I. INTRODUCTION 1

In the preceding paper, new semianalytical secondorder symplectic integrators are presented, developed by combining the molecular dynamics 共MD兲 integration2 and the standard theory of molecular vibrations.3–6 The unique feature of the new integrators is in that the standard theory of molecular vibrations, which is a very efficient tool to analyze the dynamics of the studied system from computed trajectories,7–12 is used not to analyze, but to compute trajectories of molecular systems. Information about the energy distribution of normal modes and the energy transfer between them is thus obtained without additional computations. The key property of a good MD integrator is the conservation of the system’s total energy over a long time interval. Backward error analysis13 has indicated that symplectic numerical integration methods approximately conserve the total energy of a system over time periods that are exponentially long in the size of the integration time step. Long-time conservation of the total energy by new integrators using long integration time steps is achieved by the analytical treatment of high-frequency molecular vibrations within the framework of the symplectic decomposition schemes.13–16 In this paper the new integration methods are employed to perform MD simulations of systems of nonlinear molecules with one equilibrium configuration and no internal rotation. The new integrators are superior to the standard leap-frog Verlet 共LFV兲 method17 because they allow longer integration time steps to be used for the same computational accuracy with nearly the same computational cost per integration step. II. METHODS

The Hamilton equations, which are solved for each atom of the system in MD integration, can be written in terms of Lie operators as13 a兲

Author to whom correspondence should be addressed. Electronic mail: [email protected]

0021-9606/2005/122共17兲/174102/9/$22.50

共1兲

where LˆH is the Lie operator 兵,其 is the Poisson bracket,18 and ␩ = 共q , p兲 is a vector of the coordinates and their conjugate momenta of all the particles. The formal solution of the Hamiltonian system 共1兲 is 兩␩兩tk+⌬t = exp共⌬tLˆH兲兩␩兩tk

共2兲

and it represents the exact time evolution of a trajectory in phase space composed of coordinates and momenta of all the particles from tk to tk + ⌬t, where ⌬t is the integration time step.18 A. Split integration symplectic method „SISM…

In developing new MD integration method1 we first decompose the Hamiltonian H of a system into two parts 共3兲

H = H0 + Hr ,

where H0 is the pure harmonic part of the Hamiltonian and Hr is the remaining part.19 Next, a second-order approximation for Eq. 共2兲, known as the generalized leap-frog scheme,14,15 is used

冉 冊

兩␩兩tk+1 = exp

冉 冊

⌬t ˆ ⌬t ˆ LH exp共⌬tLˆHr兲exp L H 兩 ␩ 兩 tk 2 0 2 0

+ O共⌬t3兲,

共4兲

which defines the split integration symplectic method 共SISM兲.1,20–22 The propagation by exp(共⌬t / 2兲LˆH0) is solved analytically using the normal modes of an isolated molecule,5 while the propagation by exp共⌬tLˆHr兲 is solved numerically in the same way as in the standard LFV method.17 The SISM differs from other decomposition MD integration methods in that it uses the standard theory of molecular vibrations, in particular, the concept of the Eckart frame, to define the translating and rotating internal coordinate system of a molecule for the time propagation. The

122, 174102-1

© 2005 American Institute of Physics

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-2

J. Chem. Phys. 122, 174102 共2005兲

M. Praprotnik and D. Janežič

method is described in full detail in the preceding paper.1 B. Multiple time stepping SISM „SISM-MTS…

First, we split the Hamiltonian H of the system as1 共5兲

H = H1 + H2 ,

D. Leap-frog Verlet „LFV and LFV-EQ…

To demonstrate the effectiveness of the new methods, we compared the computational performances for the same level of accuracy with the standard second-order symplectic LFV algorithm17 in which the Hamiltonian is split into the kinetic and potential energy, 共12兲

H = T + V, H1 = Vnb ,

共6兲

H2 = H0 + Vah ,

共7兲

where H0 is the pure harmonic part of the Hamiltonian, Vnb is the sum of the Coulomb and Lennard-Jones potential, and Vah is the anharmonic vibrational potential of higher terms 共cubic, quartic, etc.兲 defined in terms of the displacements of atoms from their equilibrium positions.1 The propagator exp共⌬tLˆH兲 is then approximated as

冉 冊 冋 冉 冊 冉 冊

⌬t ˆ exp共⌬tLˆH兲 = exp LV 2 nb ⫻ exp

␦t ˆ

冉 冊册

␦t ˆ LH0 exp共␦tLˆVah兲exp LH 2 2 0

⌬t ˆ LV + O共⌬t3兲, ⫻exp 2 nb

n

共8兲

In the equilibrium SISM 共SISM-EQ兲, the numerical scheme of the SISM given by Eq. 共4兲 is used to propagate the coordinates and momenta of the atoms; the potential of the slow nonbonded forces is computed with the equilibrium positions of atoms1,23–25 共9兲 共10兲

where Vnb is the sum of the Coulomb and van der Waals potentials, Fnb = −⳵Vnb / ⳵q is the corresponding force, ⳵ / ⳵q = 共⳵ / ⳵X1 , ⳵ / ⳵Y 1 , ⳵ / ⳵Z1 , . . . , ⳵ / ⳵Xn , ⳵ / ⳵Y n , ⳵ / ⳵Zn兲, q = 共q1 , . . . , q3n兲 = 共X1 , Y 1 , Z1 , . . . , Xn , Y n , Zn兲 are the Cartesian coordinates of all atoms in the system and n is the number of atoms in the system, and d共q兲 苸 R3n are the equilibrium positions of atoms in all molecules of the system, given by the standard theory of molecular vibrations.1,3,5,6 When the scheme of the SISM-MTS, defined by Eq. 共8兲, is used to propagate the coordinates and momenta of the atoms, it gives rise to the equilibrium SISM-MTS 共SISMMTS-EQ兲 method.1 This method conserves the following quantity: H = T共p兲 + Vvib共q兲 + Vnb关d共q兲兴.

共11兲

冉 冊

⌬t ˆ ⌬t ˆ LT exp共⌬tLˆV兲exp LT 兩␩兩tk + O共⌬t3兲. 2 2 共13兲

Equation 共13兲 is explicitly written as qk⬘ = qk + M−1 · pk

⌬t , 2

⳵V 共q⬘兲, ⳵q k

qk+1 = q⬘k + M−1 · pk+1

C. Equilibrium SISM „SISM-EQ…

Fnb共q兲 → JT · Fnb关d共q兲兴,

冉 冊

兩␩兩tk+1 = exp

pk+1 = pk − ⌬t

which is used to derive the multiple time stepping SISM 共SISM-MTS兲. Here ⌬t is the integration time step and ␦t = ⌬t / n is the smaller integration time step that corresponds to the time scale of high-frequency interactions defined by Vah. The propagation by exp(共␦t / 2兲LˆH0) is performed analytically in the same way as in the SISM.1

Vnb共q兲 → Vnb关d共q兲兴,

using second-order generalized leap-frog scheme15,26

⌬t , 2

共14兲

where M 苸 R3n⫻3n is a diagonal mass matrix. The diagonal elements are M 11 = m1, M 22 = m1, M 33 = m1 , . . . , M 3n−2,3n−2 = mn, M 3n−1,3n−1 = mn, M 3n,3n = mn, where mi is the mass of the ith atom. When the numerical scheme of the LFV defined by Eq. 共13兲 is used to propagate the coordinates and momenta of the atoms, and the potential of the long-range electrostatic and van der Waals potential is calculated with the equilibrium positions of the atoms in each molecule in the same way as for the SISM-EQ and SISM-MTS-EQ,1 then this gives rise to the equilibrium LFV 共LFV-EQ兲 method. III. COMPUTATIONAL DETAILS

The applicability of the SISM for MD integration is, at present, limited to systems of molecules with one equilibrium configuration and no internal rotation, and in which the displacements of atoms from their equilibrium positions are sufficiently small that we can use the dynamical molecular model6 to describe molecular vibrations. A four-atom molecule, the hydrogen peroxide 共H2O2兲, schematically shown in Fig. 1, has been chosen as an example of a nonlinear and nonplanar molecule.

A. Model potential development

For this class of molecules we first develop an appropriate model potential to be used in MD simulations of liquid H2O2 by the SISM. The H2O2 molecule, which has no center of symmetry, is one of the simplest molecules with a hindered internal rotation of the hydrogen atoms around the bond between the oxygen atoms. It has two equivalent stable equilibrium configurations 共at ±␾0兲 and two transition states, cis 共␾0 = 0 ° 兲 and trans 共␾0 = 180° 兲.27,28 The experimental

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-3

J. Chem. Phys. 122, 174102 共2005兲

Molecular dynamics integration. II

FIG. 2. 共a兲 The torsional potential defined by Eq. 共15兲 with ␾0 = 112.46°, V0 = 2Vb / 共1 − cos ␾0兲2, and Vb = 7.28 kcal/ mol corresponding to isolated molecule. 共b兲 The torsional potential defined by Eq. 共15兲 with ␾0 = 90.2° and V0 = 140 kcal/ mol defining a molecule with one equilibrium configuration and no internal rotation. FIG. 1. Description of the positions of atoms in the equilibrium configuration of hydrogen peroxide. 共a兲 Definition of angles. 共b兲 View along the bond between atoms 2 and 3.

value of ␾0 is 119.8° ± 3° for the gas phase.29 The value for ␾0 determined from ab initio calculations of an isolated H2O2 molecule is 112.46°.27 The torsional potential, also determined from ab initio calculations,27 can be fitted with the function of the harmonic cosine form V共␾兲 = 21 V0共cos ␾ − cos ␾0兲2 ,

共15兲

with V0 = 2Vb / 共1 − cos ␾0兲2, Vb = 7.28 kcal/ mol, ␾0 = 112.46°, and minima

␾min = 2n␲ ± ␾0 .

共16兲

From the torsional potential, shown in Fig. 2共a兲, it can be observed that there are two potential barriers surrounding the minimum at 112.46°, which corresponds to the equilibrium configuration. The high 7.28 kcal/ mol potential barrier corresponds to the cis transition state whereas the low 1.08 kcal/ mol potential barrier corresponds to the trans transition state. The equilibrium configuration of the H2O2 molecule in the gas state, however, does not correspond to the corresponding structure in the liquid state. The value of ␾0 in the liquid state, which is the most interesting physical system for

testing the efficiency of different numerical integrators for MD simulation, is therefore different from the corresponding value in the gas state.30 Because, at present state of development, the SISM is efficiently applicable only to systems of molecules with one equilibrium configuration and with no internal rotation, we have taken the experimentally determined structure in the solid state30 with ␾0 = 90.2° for the equilibrium configuration of the H2O2 molecule instead of the corresponding structure in the liquid state and we have set the height of potential barriers, which surround the minimum at ␾0 = 90.2°, artificially high at V0 = 140 kcal/ mol to ensure that the displacements of the hydrogens atoms are sufficiently small so that they can be considered as torsional vibrations. The corresponding torsional potential 共15兲 is depicted in Fig. 2共b兲 from which it can be observed that the heights of the potential barrier for the trans and cis transition states are equal. An H2O2 molecule with the equilibrium configuration determined by these parameters can therefore be considered as a molecule with one equilibrium configuration and no internal rotation. Hence, we have used it as an example of the class of nonlinear and nonplanar molecules in MD simulation by the SISM. The model Hamiltonian, which we have developed for MD simulation of the liquid H2O2, is

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-4

J. Chem. Phys. 122, 174102 共2005兲

M. Praprotnik and D. Janežič

TABLE I. Parameters of the Hamiltonian defined by Eq. 共17兲 for the H2O2 molecule. The quantity e0 is the elementary charge. Parameter

Value

bOH0 = b10 = b30 bOO0 = b20 ␪ 10 ␪ 20 ␾0 ␺0 = ␾0 − 90° kbOH = kb1 = kb3 kbOO = kb2 k ␪1 k ␪2 V0 eH eO ␴HH ␴OO ␧HH ␧OO

0.988 Å 1.453 Å 102.7° 102.7° 90.2° 0.2° 900.0 kcal/ mol/ Å2 580.0 kcal/ mol/ Å2 140.0 kcal/ mol/ radian2 140.0 kcal/ mol/ radian2 140.0 kcal/ mol 0.350 53e0 −0.350 53e0 0.40 Å 3.1507 Å 0.045 98 kcal/ mol 0.152 073 kcal/ mol

TABLE II. Experimental vibration frequencies 共Refs. 34–36兲 of liquid H2O2 and normal mode frequencies of the H2O2 molecule determined by normal mode analysis using the parameters from Table I.

a

Normal mode

1 / ␭ 共cm−1兲 共experiment兲a

1 / ␭ 共cm−1兲 共theory兲b

Antisymmetric O–H stretch Symmetric O–H stretch Symmetric angle bending Antisymmetric angle bending O–O stretch Torsional oscillation

3360 3360 1421 1350 878 635

3358 3357 1410 1386 880 1965

Experimental vibration frequencies. Normal mode frequencies.

b

B. Vibrational potential energy and internal coordinate system

The vibrational potential energy is the sum of the vibrational potential energies of all of the molecules in the system, m

Vvib = H=

p2

1

+

1 e ie j V0共cos ␾ − cos ␾0兲2 + 2 torsions i⬎j 4␲⑀0rij

+

4␧ij 兺 i⬎j



冋冉 冊 冉 冊 册 ␴ij rij

12



␴ij rij

+



6

,

共17兲

where i and j run over all atoms, mi is the mass of the ith atom, pi is the linear momentum of the ith atom, b0 and ␪0 are reference values for the bond lengths and angles, respectively, kb and k␪ are the corresponding force constants, ␾0 is the reference value for the torsional angle, V0 is the corresponding barrier height, ei denotes the charge on the ith atom, ⑀0 is the dielectric constant in vacuum, rij is the distance between the ith and jth atoms, and ␧ij and ␴ij are the corresponding constants of the Lennard-Jones potential, and the Lorentz–Berthelot mixing rules are used.2 The van der Waals and Coulomb interactions between hydrogen atoms of the same molecule are explicitly taken into account. For the reference values of the bonds, angles, and torsional angles we have taken the experimental values from Ref. 30. For the constants of the Lennard-Jones potential we have taken the corresponding values of flexible TIP3P water model31,32 and the partial charges of the atoms were calculated from the ab initio calculated dipole moment33 and the corresponding structure of an isolated H2O2 molecule.27 The force constants for bond stretching and angle bending were determined by fitting the normal mode frequencies calculated by normal mode analysis1 to the experimental frequencies in the IR and Raman spectrum of liquid H2O2.34–36 The parameters of the Hamiltonian 共17兲 are reported in Table I. The experimental and calculated normal mode frequencies for the H2O2 molecule using these parameters are reported in Table II.

1

k

1

兺i 2mi i + 2 bonds 兺 kb共b − b0兲2 + 2 angles 兺 k ␪共 ␪ − ␪ 0兲 2

1

Vvib = 兺 kb共b − b0兲2 + 兺 k␪共␪ − ␪0兲2 兺 2 bonds 2 angles k=1 1 V0共cos ␾ − cos ␾0兲2 , 2 torsions



共18兲

where Vvibk is the vibrational potential energy of the kth molecule in the system and m is the number of molecules in the system. In the SISM, the propagation by exp(共⌬t / 2兲LˆH0) is integrated analytically using the normal coordinates to describe the vibrational, rotational, and translational degrees of freedom of each molecule in the system. For the transformation of Cartesian coordinates and momenta into the normal coordinates, the relative Cartesian displacement coordinates are required. To determine the vibrational frequencies and normal modes of vibration of the kth molecule in the system, the mass-weighted Hessian M−1/2 · Hk · M−1/2 苸 R3N⫻3N has to be diagonalized. The matrix Hk is a symmetric matrix with the elements Hkij = Hk ji =



⳵2Vharmk ⳵⌬qi⳵⌬q j



共19兲 0

and M is a diagonal mass matrix with the elements M 11 = m1, M 22 = m1, M 33 = m1 , . . . , M 3N−2,3N−2 = mN, M 3N−1,3N−1 = mN, M 3N,3N = mN, where N is the number of all of the atoms in the kth molecule. The harmonic vibrational potential energy Vharmk for the kth molecule is defined as

冉 兺冉

3N ⳵2Vvibk 1 Vharmk = 2 i,j=1 ⳵⌬qi⳵⌬q j



=

3N ⳵2Vharmk 1 2 i,j=1 ⳵⌬qi⳵⌬q j

冊 冊

⌬qi⌬q j 0

⌬qi⌬q j 0

3N

=

1 1 Hk ⌬qi⌬q j = ⌬q · Hk · ⌬q, 2 i,j=1 ij 2



共20兲

where ⌬q = 共⌬x1 , ⌬y 1 , ⌬z1 , . . . , ⌬xN , ⌬y N , ⌬zN兲 is a vector of

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-5

J. Chem. Phys. 122, 174102 共2005兲

Molecular dynamics integration. II

the relative Cartesian displacement coordinates.1 The Hessian Hk as well as the functional form of Vharmk are equal for every molecule in the system and therefore the index k is omitted. The harmonic potential Vharm is then expressed as 共21兲

Vharm = Vstretch + Vbend + Vtorsion ,

where Vstretch, Vbend, and Vtorsion are the bond stretching, angle bending, and torsional potentials, respectively, expressed as a quadratic forms in terms of relative Cartesian displacement coordinates. The equilibrium configuration of the H2O2 molecule is shown in Fig. 1 where atom 1 is the first hydrogen atom,

atoms 2 and 3 are oxygen atoms, and atom 4 is the second hydrogen atom, respectively. The definitions of the angles ␪10, ␪20, and ␺0 are evident from Fig. 1共a兲, and the distances are b10 = b30 = bOH0 and b20 = bOO0. The unit vector pointing from atom 2 to atom 1 is then 共cos ␪10 , −sin ␺0 sin ␪10 , cos ␺0 sin ␪10兲, the unit vector pointing from atom 2 to atom 3 is 共1, 0, 0兲, and the unit vector pointing from atom 3 to atom 4 is 共−cos ␪20 , sin ␪20 , 0兲. The expression for Vstretch in terms of the relative Cartesian coordinates is obtained by projecting the difference of the displacements of the atoms onto the unit vectors, which point along the bonds between atoms.37 Then

Vstretch = 21 kb1关共⌬x1 − ⌬x2,⌬y 1 − ⌬y 2,⌬z1 − ⌬z2兲 · 共cos ␪10,− sin ␺0 sin ␪10,cos ␺0 sin ␪10兲T兴2 + 21 kb2关共⌬x3 − ⌬x2,⌬y 3 − ⌬y 2,⌬z3 − ⌬z2兲 · 共1,0,0兲T兴2 + 21 kb3关共⌬x4 − ⌬x3,⌬y 4 − ⌬y 3,⌬z4 − ⌬z3兲 · 共− cos ␪20,sin ␪20,0兲T兴2 = 21 kb1关共⌬x1 − ⌬x2兲cos ␪10 − 共⌬y 1 − ⌬y 2兲sin ␺0 sin ␪10 + 共⌬z1 − ⌬z2兲cos ␺0 sin ␪10兴2 + 21 kb2共⌬x3 − ⌬x2兲2 + 21 kb3关− 共⌬x4 − ⌬x3兲cos ␪20 + 共⌬y 4 − ⌬y 3兲sin ␪20兴2 ,

共22兲

where · denotes the dot product of two vectors, kb1 = kb3 = kbOH, and kb2 = kbOO. Similarly, the expression for Vbend in terms of the relative Cartesian coordinates are obtained by taking the components of the difference of the displacements of the atoms perpendicular to the bonds between the atoms.37 Therefore



1 1 1 Vbend = k␪1 共⌬x1 − ⌬x2,⌬y 1 − ⌬y 2,⌬z1 − ⌬z2兲 · 共− sin ␪10,− sin ␺0 cos ␪10,cos ␺0 cos ␪10兲T + 共⌬x3 − ⌬x2,⌬y 3 2 b 10 b 20 − ⌬y 2,⌬z3 − ⌬z2兲 · 共0,sin ␺0,− cos ␺0兲T +

册 冋 2

1 1 + k ␪2 共⌬x2 − ⌬x3,⌬y 2 − ⌬y 3,⌬z2 − ⌬z3兲 · 共0,− 1,0兲T 2 b 20

1 共⌬x4 − ⌬x3,⌬y 4 − ⌬y 3,⌬z4 − ⌬z3兲 · 共sin ␪20,cos ␪20,0兲T b 30





2

1 1 = k ␪1 关− 共⌬x1 − ⌬x2兲sin ␪10 − 共⌬y 1 − ⌬y 2兲sin ␺0 cos ␪10 + 共⌬z1 − ⌬z2兲cos ␺0 cos ␪10兴 2 b 10 +

1 关共⌬y 3 − ⌬y 2兲sin ␺0 − 共⌬z3 − ⌬z2兲cos ␺0兴 b 20





2

1 1 1 + k ␪2 − 共⌬y 2 − ⌬y 3兲 + 关共⌬x4 − ⌬x3兲sin ␪20 + 共⌬y 4 − ⌬y 3兲cos ␪20兴 2 b 20 b 30

where k␪1 = k␪2 = k␪OOH. To express the torsional angle ␾ in terms of the relative Cartesian coordinates we define



2

,

共23兲

a34 = − 共− b30 cos ␪20 + ⌬x4 − ⌬x3,b30 sin ␪20 + ⌬y 4 − ⌬y 3,⌬z4 − ⌬z3兲,

共26兲

where a␣␤ = a␣ + ␳␣ − 共a␤ + ␳␤兲, ␳␣ = 共⌬x␣ , ⌬y ␣ , ⌬z␣兲, ␣ , ␤ = 1 , . . . , 4, and

a12 = − 共− b10 cos ␪10 + ⌬x2 − ⌬x1,b10 sin ␺0 sin ␪10 + ⌬y 2 − ⌬y 1,− b10 cos ␺0 sin ␪10 + ⌬z2 − ⌬z1兲, 共24兲

a1 = 共b10 cos ␪10,− b10 sin ␺0 sin ␪10,b10 cos ␺0 sin ␪10兲, 共27兲

a23 = − 共b20 + ⌬x3 − ⌬x2,⌬y 3 − ⌬x2,⌬z3 − ⌬z2兲,

共25兲

a2 = 共0,0,0兲,

共28兲

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-6

J. Chem. Phys. 122, 174102 共2005兲

M. Praprotnik and D. Janežič

a3 = 共b20,0,0兲,

共29兲

a4 = 共b20 − b30 cos ␪20,b30 sin ␪20,0兲.

共30兲

cos ␾ is then obtained as

2

cos ␾ =

共a12 ⫻ a23兲 · 共a23 ⫻ a34兲 , 兩a12 ⫻ a23兩兩a23 ⫻ a34兩

共31兲

where 兩·兩 denotes the vector norm. In Eq. 共31兲 we only keep the linear and quadratic terms in relative Cartesian displacement coordinates. Then



1 1 Vtorsion = V0 − 共⌬y 2 − ⌬y 1兲 2 b10 sin ␪10 +

1 关− 2 sin ␺0共⌬x3 − ⌬x2兲 − 共cot ␪10 b 20

+ sin ␺0 cot ␪20兲共⌬y 3 − ⌬y 2兲 + cos ␺0 ⫻cot ␪20共⌬z3 − ⌬z2兲兴 +

1 关cos ␺0共⌬z4 − ⌬z3兲 b30 sin ␪20

− sin ␺0共⌬y 4 − ⌬y 3兲兴



2

共32兲

.

The elements of the Hessian H are determined by Eq. 共19兲. The dimension of the mass-weighted Hessian M−1/2 · H · M−1/2 is 12⫻ 12 in the case of the H2O2 molecule. The 3N − 6 = 6 vibrational normal mode frequencies for the H2O2 molecule are given in Table II. The mass-weighted Hessian M−1/2 · H · M−1/2 was diagonalized using subroutines TRED2 and TQLI taken from Ref. 38. These subroutines were also used to diagonalize the symmetric positive definite Gram matrix F.1 The translating and rotating internal coordinate system is defined by the right-handed triad of unit vectors fi, i = 1 , 2 , 3, where f j · fk = ␦ jk, with the origin in the center of mass of a molecule. The constant equilibrium distances of the atoms from the molecule’s center of mass ci␣, i = 1 , 2 , 3.1 which are required for defining of the internal coordinate system1 are then obtained from vectors c␣, ␣ = 1 , . . . , 4, determined as 3

c␣ = a␣ − R =

c␣i fi = 共c1␣,c2␣,c3␣兲, 兺 i=1

共33兲

where R is

兺 ␣ m ␣a ␣ 1 R= = 关m1b1 m + m + 兺 ␣ m␣ 1 2 m3 + m4

0

cos ␪10

+ m3b20 + m4共b20 − b30 cos ␪20兲, − m1b10 sin ␺0 sin ␪10 + m4b30 sin ␪20, m1b10 cos ␺0 sin ␪10兴.

共34兲

The matrix F−1/2, which is required to determine the unit vectors of the internal coordinate system of a molecule in the SISM1, is calculated as F−1/2 = P · D−1/2 · PT ,

共35兲

where D−1/2 is a diagonal matrix with the elements D−1/2 ii = 1 / 冑␭i and ␭i are the eigenvalues of the symmetric positive definite Gram matrix F. The columns of the transition matrix P are the eigenvectors of F and PT is the transpose of P. C. Simulation protocol

We have carried out MD simulation of a system of 256 H2O2 molecules with the density ␳ = 1.4425 g / cm3 at T = 298 K corresponding to the liquid state.39 The corresponding size of the simulation box was a = 21.6 Å. Periodic boundary conditions were imposed to overcome the problem of surface effects; the minimum image convention was used.2 The Coulomb interactions were truncated using the force-shifted potential40 with a cutoff distance rof f = 8.5 Å.41 The Lennard-Jones interactions were shifted by adding the term Cijr6ij + Dij to the potential, where Cij and Dij were chosen such that the potential and force are zero at rij = rof f .32 The initial positions and velocities of the atoms were chosen at random. The system was then equilibrated for 50 ps where the velocities were scaled every 500 integration time steps, followed by an additional 50 ps of equilibration at constant energy of the system to ensure that the velocities assume the Maxwell distribution at T = 298 K. To obtain physically and numerically relevant initial conditions to perform the MD simulation of a system of flexible molecules, the equilibration was also monitored using the Vieillard–Baron rotational order parameter.22,42 IV. RESULTS AND DISCUSSION

To demonstrate the effectiveness of the SISM, in all our numerical experiments we compared the computational performances for the same level of accuracy with the standard second-order LFV algorithm using an integration time step small enough to accurately describe the high-frequency molecular vibrations. In this way it is assured that the physical properties of the system determined from the trajectories computed by new integrators using long integration time steps are reliable. For that purpose the error in total energy ⌬E / E defined as M





E0 − Ek ⌬E 1 = , E M k=1 E0



共36兲

where E0 is the initial energy, Ek is the total energy of the system at the integration step k, and M is the total number of integration steps, was monitored for all methods and was used as a measure of the efficiency and accuracy of numerical integrators for MD simulation. The speed-up of the new methods due to a prolongation of the integration time step can be determined from the error in total energy, which is depicted in Fig. 3共a兲 for the system of 256 H2O2 molecules for the LFV, SISM, and SISM-MTS.

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-7

Molecular dynamics integration. II

FIG. 3. 共a兲 The error in the total energy of the system of 256H2O2 molecules with ␳ = 1.4425g / cm3 at T = 298 K using the LFV, SISM, and SISM-MTS for M = 1000. 共b兲 The error in the total energy of the system of 256H2O2 molecules with ␳ = 1.4425 g / cm3 at T = 298 K using the LFV-EQ, SISM-EQ, and SISM-MTS-EQ for M = 1000.

The period for the antisymmetric stretching of the bond between the oxygen and the hydrogen atom in the H2O2 molecule is 9.9 fs 共Table II兲. We estimate that the maximal acceptable size of the integration time step for the LFV to be 0.5 fs. From results in Fig. 3共a兲 we conclude that the error in total energy for a 1.25 fs integration time step in the case of the SISM corresponds to the error in the total energy using a 0.5 fs integration time step in the case of the LFV. This means that the SISM allows the use of an up to a two and a half times longer time step than the LFV for the same level of accuracy. For the system of H2O2 molecules, a high amount of anharmonic forces derived from Vah is expected due to the strong anharmonic potential describing the interactions in the system. Surprisingly enough the SISM-MTS, which employs ␦t = 0.5 fs for the integration of motions generated by highfrequency anharmonic interactions is less accurate than the SISM 关Fig. 3共a兲兴. This can be explained by the fact that the large amount of anharmonic forces stems from the strong electrostatic and van der Waals interactions and not from the high-frequency anharmonic interactions defined by Vah. Due to strong intermolecular forces we expect good performance of the SISM-EQ and SISM-MTS-EQ integrators. This is confirmed by the results in Fig. 3共b兲 in which the error in total energy for the LFV-EQ, SISM-EQ, and SISM-MTS-EQ is shown. The total energy in this case equals expression

J. Chem. Phys. 122, 174102 共2005兲

共11兲. Since the vibrational potential Vvib taken into account in the LFV-EQ is the same as in the LFV, the estimated value of the maximal acceptable time step for the LFV-EQ is 0.5 fs. From Fig. 3共b兲 we can determine that the error in total energy for the LFV-EQ with a 0.5 fs integration time step corresponds to the error in the case of the SISM-EQ using a 1.5 fs integration time step or 2.0 fs in the case of the SISMMTS-EQ. The SISM-EQ therefore allows the use of a three times longer integration time step than the LFV-EQ for the same level of accuracy, whereas the SISM-MTS-EQ allows the use of even up to a four times longer time step as the LFV-EQ. We can also conclude that the SISM becomes unstable for integration time steps longer than 3.75 fs whereas no drift occurs in total energy using the SISM-MTS-EQ with integration time steps shorter than 5.0 fs. Using longer time steps results in a drift in the total energy, which is consistent with the conclusions in Ref. 25 where a linear numerical instability is predicted for the integration time step size corresponding to around half of the period of the fastest normal mode. The prolongation of the maximal acceptable integration time step by the SISM-EQ and SISM-MTS-EQ in comparison to the LFV-EQ comes from the fact that the maximal acceptable integration time step is limited by the atoms’ motion generated by the intermolecular forces in the case of the SISM-EQ and SISM-MTS-EQ.22 The maximal integration time step allowed by the LFV-EQ method, however, is limited by the intramolecular high-frequency vibrations that are on the considerably smaller time scale. We can also report resonances occurring in the error in total energy when the size of the integration time step corresponds to a multiple of the period of the fastest normal mode of a molecule 共not shown in Figs. 3共a兲 and 3共b兲 where only physically meaningful integration time step sizes are considered兲 as in the case of the Verlet-I/r-RESPA method.26,43 The Lennard-Jones and electrostatic interactions represent the external driving forces on the internal motion of the molecules and resonance always occurs if the frequency of the driving force corresponds to a multiple of the oscillator frequency, regardless of whether the high-frequency molecular vibrations are integrated numerically, as in the case of the VerletI/r-RESPA method, or analytically, as in the case of the SISM and its derivatives, which are also not stable in these resonances. Since the period of the fastest motion is 9.9 fs, we can draw the conclusion that the size of the maximal integration time step of the new methods is also limited in this case by the nonlinear instabilities at a third or fourth of the period of the fastest motion and by linear instabilities at half of the period of the fastest motion, as found for the Verlet-I/rRESPA method.25,44 The actual speed-up of the new methods is confirmed by measuring the CPU time spent by the methods per integration step. The CPU times for the three methods 共the SISM, SISM-MTS, and LFV兲 for 1000 MD steps measured on an AMD Athlon XP 1600+ processor for different system sizes 共m兲 and equal time steps 共1 fs兲 are given in Table III. The results indicate that the computation cost per integration step is slightly larger for the SISM and SISM-MTS than for the LFV. However, for larger systems consisting of

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-8

J. Chem. Phys. 122, 174102 共2005兲

M. Praprotnik and D. Janežič

TABLE III. CPU time 共s兲 for SISM and LFV for 1000 MD steps measured on an AMD Athlon XP 1600+ processor for different system sizes 共m兲 and equal integration time steps 共1 fs兲. m

t共LFV兲 共s兲

t共SISM兲 共s兲

t共SISM-MTS兲 共s兲

64 128 256

14.21 49.57 170.52

16.39 54.51 182.65

18.72 58.51 188.64

more than 100 molecules the cost of computation per integration step becomes approximately the same for all methods because the computation of long-range forces, which is the same for all methods, prevails over the computation of extra transformations required by the SISM and SISM-MTS. Therefore the significant speed-up of the SISM and SISMMTS over the LFV is due to the larger integration time step allowed by the SISM and SISM-MTS. The applied rotational barrier of 140 kcal/ mol is much higher than the estimated experimental value of ⬇4 kcal/ mol for the potential barrier in the trans transition state.34 For efficient performance of the SISM in the case of the realistic rotational barrier, the normal mode associated with the torsional oscillation should be treated as the internal rotation and should be therefore excluded from the description by the normal coordinates.4,5 In addition, an additional internal coordinate system should be introduced to describe the internal rotation of the molecule.45–49 Further verification that the new integrators yield correct trajectories using long integration time step can be obtained by computing various statistical properties of a molecular system, e.g., radial distribution functions, diffusion coefficients, orientational correlation times, vibrational spectra, for different integration time step sizes. We have performed this for a system of liquid water at ambient conditions. The results are presented in the next paper of this series50 where we give further evidence that correct dynamics is achieved by the new integrators.

V. CONCLUSIONS

In the present paper we have tested the efficiency of newly introduced semianalytical symplectic integration methods for MD simulation on systems of nonlinear and nonplanar hydrogen peroxide molecules. The new integration methods are also stable even when using integration time steps several times longer than the standard LFV method. The numerical results show that the new integration methods allow the use of up to four times longer integration time step than the standard LFV method for the same level of accuracy. The results also indicate that the computation cost per integration time step of newly developed MD integrators is basically the same as that of the standard method. Therefore the up to a fourfold simulation speed up of new symplectic integrators is due to the larger time step they allow to use. However, much work remains to be done in the development of this approach to explore its advantages and limitations.

t共SISM兲

t共SISM-MTS兲

t共LFV兲

t共LFV兲

1.15 1.10 1.07

1.32 1.18 1.11

ACKNOWLEDGMENTS

The authors thank Dr. M. Hodošček, and Dr. F. Merzel for helpful discussions, and U. Borštnik for careful reading of the manuscript. This work was supported by the Ministry of Education, Science and Sports of Slovenia under Grant Nos. P1-0002 and J1-6331. 1

D. Janežič, M. Praprotnik, and F. Merzel, J. Chem. Phys. 122, 174101 共2005兲, preceding paper. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids 共Clarendon, Oxford, 1987兲. 3 C. Eckart, Phys. Rev. 47, 552 共1935兲. 4 A. Sayvetz, J. Chem. Phys. 7, 383 共1939兲. 5 E. B. Wilson, J. C. Decius, and P. C. Cross, Molecular Vibrations 共McGraw-Hill, New York, 1955兲. 6 J. D. Louck and H. W. Galbraith, Rev. Mod. Phys. 48, 69 共1976兲. 7 R. Rey, J. Chem. Phys. 104, 2356 共1996兲. 8 R. Rey, Chem. Phys. 229, 217 共1998兲. 9 R. Rey, J. Chem. Phys. 108, 142 共1998兲. 10 B. R. Brooks, D. Janežič, and M. Karplus, J. Comput. Chem. 16, 1522 共1995兲. 11 D. Janežič and B. R. Brooks, J. Comput. Chem. 16, 1543 共1995兲. 12 D. Janežič, R. M. Venable, and B. R. Brooks, J. Comput. Chem. 16, 1554 共1995兲. 13 J. M. Sanz-Serna and M. P. Calvo, Numerical Hamiltonian Problems 共Chapman & Hall, London, 1994兲. 14 H. F. Trotter, Proc. Am. Math. Soc. 10, 545 共1959兲. 15 G. Strang, SIAM 共Soc. Ind. Appl. Math.兲 J. Numer. Anal. 5, 506 共1968兲. 16 H. Yoshida, Phys. Lett. A 150, 262 共1990兲. 17 L. Verlet, Phys. Rev. 159, 98 共1967兲. 18 H. Goldstein, Classical Mechanics, 2nd ed. 共Addison-Wesley, New York, 1980兲. 19 D. Janežič and F. Merzel, J. Chem. Inf. Comput. Sci. 35, 321 共1995兲. 20 D. Janežič and M. Praprotnik, Int. J. Quantum Chem. 84, 2 共2001兲. 21 M. Praprotnik and D. Janežič, Cell. Mol. Biol. Lett. 7, 147 共2002兲. 22 D. Janežič and M. Praprotnik, J. Chem. Inf. Comput. Sci. 43, 1922 共2003兲. 23 B. Garcia-Archilla, J. M. Sanz-Serna, and R. D. Skeel, SIAM J. Sci. Comput. 共USA兲 20, 930 共1998兲. 24 J. A. Izaguirre, S. Reich, and R. D. Skeel, J. Chem. Phys. 110, 9853 共1999兲. 25 R. D. Skeel and J. A. Izaguirre, in Computational Molecular Dynamics: Challenges, Methods, Ideas, Lecture Notes in Computational Science and Engineering Vol. 4, edited by P. Deuflhard, J. Hermans, B. Leimkuhler, A. E. Mark, S. Reich, R. D. Skeel 共Springer, Berlin, 1999兲, pp. 318–331. 26 H. Grubmüller, H. Heller, A. Windemuth, and K. Schulten, Mol. Simul. 6, 121 共1991兲. 27 J. Koput, Chem. Phys. Lett. 236, 516 共1995兲. 28 V. A. Benderskii, I. S. Irgibaeva, E. V. Vetoshkin, and H. P. Trommsdorff, Chem. Phys. 262, 369 共2000兲. 29 R. L. Redington, W. B. Olson, and P. C. Cross, J. Chem. Phys. 36, 1311 共1962兲. 30 W. R. Busing and H. A. Levy, J. Chem. Phys. 42, 3054 共1965兲. 31 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 共1983兲. 32 P. J. Steinbach and B. R. Brooks, J. Comput. Chem. 15, 667 共1994兲. 33 J. Koput, Chem. Phys. Lett. 257, 36 共1996兲. 34 P. A. Giguere, J. Chem. Phys. 18, 88 共1950兲. 35 R. C. Taylor, J. Chem. Phys. 18, 898 共1950兲. 2

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp

174102-9

J. Chem. Phys. 122, 174102 共2005兲

Molecular dynamics integration. II

O. Bain and P. A. Giguere, Can. J. Chem. 33, 527 共1955兲. L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics: Mechanics, 3rd ed. 共Pergamon, Oxford, 1976兲, vol. 1. 38 W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing 共Cambridge University Press, Cambridge, 1987兲. 39 CRC Handbook of Chemistry and Physics, 66th ed., edited by R. C. Weast 共CRC, Cleveland, 1986兲. 40 C. L. Brooks III, B. M. Pettitt, and M. Karplus, J. Chem. Phys. 83, 5897 共1985兲. 41 M. Prevost, D. van Belle, G. Lippens, and S. Wodak, Mol. Phys. 71, 587 共1990兲.

J. E. Vieillard-Baron, J. Chem. Phys. 56, 4729 共1967兲. M. E. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. 97, 1990 共1992兲. 44 Q. Ma, J. A. Izaguirre, and R. D. Skeel, SIAM J. Sci. Comput. 共USA兲 24, 1951 共2003兲. 45 J. B. Howard, J. Chem. Phys. 5, 442 共1937兲. 46 J. B. Howard, J. Chem. Phys. 5, 451 共1937兲. 47 B. Kirtman, J. Chem. Phys. 37, 2516 共1962兲. 48 B. Kirtman, J. Chem. Phys. 41, 775 共1964兲. 49 B. Kirtman, J. Chem. Phys. 49, 2257 共1968兲. 50 M. Praprotnik and D. Janežič, J. Chem. Phys. 122, 174103 共2005兲, following paper.

36

42

37

43

Downloaded 09 Aug 2005 to 194.95.63.241. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp