Molecular dynamics of polymers with explicit but

M OLECULAR P HYSICS, 2001, VOL. 99, N O. 3, 155 –165

Molecular dynamics of polymers with explicit but frozen hydrogens JEAN -PAU L R YCK AER T 1 *, G IAN LU IG I AR IALD I 1 and SIM ON E M ELCHION N A 2 1 Laboratoire de Physique des Polyme`res, CP 223, U niversite´ Libre de Bruxelles, Bd du Triomphe, B-1050 Brussels, Belgium 2 D epartment of Chemistry, Lens eld R oad, Cambridge CB2 1EW, U K (R eceived 18 M ay 2000; accepted 4 S eptember 2000) A new form of holonomic constraint, called a conic constraint, is introduced for the purpose of eliminating the fast vibrations of hydrogen atoms in molecular simulations of systems of aliphatic chains. It can easily be combined with bond constraints in SH AK E/R ATTLE algorithms for which a uni ed tolerance criterion is de ned. The new form of constraint allows the use of rather large time steps (in the 2–3 fs range). The procedure is illustrated for a full atomic model of polyprop ylene.

1. Introduction M otions of the hydrogen atoms are responsible for the highest frequency vibrations in hydrogen ated molecules, including hydrocarbon chains and biological macromolecules. These motions place an upper limit of about 1 fs on the size of the time step which can be used in molecular dynamics (M D ) simulations. When the aim of the calculations is to probe dynamical events over much larger time scales, an important example being the local dynamics relaxation of polymer melts well above T g which takes place in the 0.1–10 ns time window, it is advantageous to eliminate the fast vibrating motions of hydrogens. In the rst M D studies of hydrocarbon chains, the high-frequency hydrogen motions were eliminated by treating methyl (CH 3 ), methylene (CH 2 ) and methine (CH ) groups as ‘united atoms’ and the bond C–C stretching vibrations were eliminated by use of bondlength constraints [1, 2]. Such models, in which explicit account is taken of C–C–C bending and torsion modes, continue to be used in the study of saturated hydrocarbon melts [3, 4]. H owever, models with explicit hydrogen s are also widely used not only for hydrocarbon solid phases, where the hydrogen centres of force are needed to stabilize the experimental crystalline structure [5], but also in the liquid phase, where a fully atomistic model gives better agreement with experimental data [6, 7]. In order to introduce explicit hydrogen atoms in semi exible chains in a way which copes e ciently with their associated fast vibrations, di erent strategies have been

* Author for correspondence. e-mail: [email protected]

adopted. The rst approach historically was to use holonomic constraints to eliminate the fast motions. F or saturated hydrocarbon chains, the four atoms of a methyl end chain group and the carbon atom directly adjacent to it were treated as a rigid unit of ve particles linked by nine or ten bond constraints (for the latter option, see the discussion on the redundant constraints in the test section 5.1 of the present paper and in [8]). F or methylene hydrogens, it was required that the bisector of the H –C–H rigid triangle should lie along the bisector of the adjacent C–C–C bending angle while an additional constraint forces the H –H bond to remain permanently orthogonal to the C–C–C plane [9]. To our knowledge, no constraint scheme has so far been proposed to eliminate vibrations associated with methine hydrogen s in –CH R – groups (see gure 1) while retaining the bending vibrations for the pending group R , e.g. a methyl group in polypropylene (PP), a benzene ring in polystyrene (PS), a chlorine atom in polyvinylchloride (PVC), or a side chain in a branched alkane. In simulations of PVC, for example [10, 11], the fast vibrations associated with methine hydrogen s were arti cially slowed down by redistributing the masses of the –CH Cl– fragment without changing the sum of the masses and the moment of inertia around the C–C bonds. R ecently, Antoniadis et al. [12] has suggested treating the methylene and methine hydrogens in atactic polypropylene (aPP) only as centres of forces, the forces on these geometrical points being redistributed on the adjacent carbon nuclei. A disadvantage of this procedure is that it eliminates the inertial in uence of the hydrogen nuclei on the skeletal dynamics. M oreover, the absence of the methine hydrogen from the model leads

M olecular Physics ISSN 0026–8976 print /ISSN 1362–3028 online # 2001 Taylor & Francis Ltd http://www.tandf.co.uk/journals D OI: 10.1080/0026897001000728 0

156

J.-P. R yckaert et al.

R

H d CH

a w

C i v-

u b /2 b /2 d CC

d CC C

Ci - 1

i+1

b F igure 1. Geometry of a hydrocarbon chain fragment for the de nition of the conic constraints associated with one skeletal hydrogen atom, the R particle being either an isolated atom (hydrogen, chlorine,. . .) or the rst atom of a side group (carbon of a methyl group, carbon of a phenyl ring,. . .). The angle ¬ between the C i –H bond and the bisector b of the variable bending angle is constant, the C i –H bond and the bisector b de ning a plane perpendicular to the C i¡1 –C i –C i‡1 triangle (see text). The basis of the conic constraint is the fact that the C i –H and C i –C i¡1 (or C i –C i‡1 ) bond s de ne an angle which is determined by the instantaneous bending angle (see text).

to accidental chiral inversions of the methine carbon, in agreement with what we independently observed for a united atom PP model [13]. The problem is that no realistic force constant associated with bending C–C–C modes in PP is large enough to prevent the rst atom of the pendent side group to exchange its position with the ‘virtual’ hydrogen . F orcing a speci c chirality to a methine carbon atom requires an extra, ad hoc term in the potential [12, 13] or the use of a so-called improper dihedral angle potential contribution [14] which requires a speci c reparametrization of the model. Another route which can be followed is to retain all vibrations but use a multiple time-step method to update the atomic interactions once they have been divided into appropriate classes, such as ‘soft’ (e.g. van der Waals) and ‘hard’ (e.g. chemical bond) potentials [15]. The method uses small time steps for the hard, vibrational, degrees of freedom while the soft forces (which are the most time consuming in CPU time) need to be recalculated only once every n steps, where n is an integer parameter to be chosen in an optimal way. This elegant and rigorous approach is certainly

of value when short time dynamics (º 0:1 ps time scale) are of interest, but it does have some drawbacks. Equipartition of kinetic energy over all degrees of freedom can be rather slow to establish because of the weak coupling between soft and hard degrees of freedom, which in turn gives rise to an ill-de ned temperature. The problem can be cured at the price of adding a N ose´ thermostat, or even a chain of N ose´ thermostats [16], but this and the multiple time-step procedure itself make the method rather tricky to program while the choice of parameters controlling the e ciency of the integration algorithm requires experience. In the present paper we propose a new scheme which makes use of holonomic constraints and is computationally simple to implement. The integration of the molecular dynamics trajectory of the constrained system can be implemented within the traditional SH AK E/ R ATTLE procedure combined with the Verlet algorithm [17, 18]. We aim at freezing all fast degrees of freedom associated with the hydrogen atoms in methylene and methine carbon groups (e.g. in vinyl polymers) for a semi- exible chain model with explicit skeletal bending and torsional modes. In such a model, because the distance between the hydrogen atom and its second neighbour main chain atom is not constant, the desired reduction in the number of degrees of freedom cannot be achieved by means of bond constraints alone. Our new procedure is based on the combination of two types of constraint only. The rst of these is the traditional bond constraint, which is applied to the C–H and C–C bonds. The second is a new, scalar constraint, the origin of which is the fact that the angle between the C–H bond and one adjacent C–C bond is xed by the instantaneous value of the adjacent C–C–C bending angle. U se of the new constraint, which we call a conic constraint, has two advantages. F irst, its treatment is only marginally more intricate mathematically than that of the usual bond constraint. Secondly, such a unique additional scalar constraint treats all skeletal hydrogens, both methylene and methine hydrogens, on an equal footing, making its inclusion in the SH AK E and R ATTLE routines potentially useful for a wide range of hydrogen ated polymers. As an illustration of the new method, we have applied it to a model of an atactic polypropylene melt in which the saturated hydrocarbon molecule is represented as a semi- exible chain having exible bending angles in the backbone. In this particular case, we nd a gain by a factor two to three in the M D step size when the conic constraint is used. As a result, time steps of 2–3 fs are usable in realistic polymer simulations in which hydrogen atoms are explicit point masses carrying centres of force. We also propose a unique relative tolerance criterion for constraint satisfaction, valid for all types of constraints, in the control of the convergence of both the

157

M olecular dynamics of polymers with ex plicit but frozen hydrogens SH AK E and R ATTLE iterative procedures. We stress that in certain cases, such as methyl groups made fully rigid by bond constraints, the use of redundant constraints can improve dramatically the rate of convergence of SH AK E and R ATTLE routines. 2. General equations of motion Consider a set of N classical particles (atoms or united atoms) interacting via an arbitrary intra and inter-molecular force eld U … frg† subject to certain holonomic constraints. The equation of motion of particle i is d 2 ri 1 ˆ … Fi ‡ Gi† ; 2 mi dt

… 1†

where Fi is the force on i coming from the potential U and Gi is the total force of constraint acting on i. A typical simulation treats a set of molecules in interaction. F or our purposes, it is su cient to focus on a single molecule consisting of N atoms linked by a set of L holonomic (scalar) constraints of the form ¼k… frgnk † ˆ 0;

k ˆ 1;. . . ;L ;

… 2†

where frgnk represents a subset of nk atomic coordinates. The time derivative of equation (2) may be written as X X Ñ i ¼k… frgnk † ¢ r_ i ² Ak ;i ¢ r_ i ˆ 0; … 3† ¼_ k… frgnk † ˆ i

i2nk

where Ak ;i is the gradient of constraint k with respect to the coordinates of atom i; the symbol i 2 nk means that the sum on i is restricted to the nk atoms involved in the constraint k. The total constraining force on atom i can be written formally as Gi ˆ ¡

L X

¶k Ak ;i ;

… 4†

kˆ1

where the sum on k runs over all constraints in which atom i is involved. The parameters ¶ k are a set of L Lagrange multipliers which vary with time through the imposition of the L constraints on the mechanical trajectory. The constraining forces associated with the molecule under consideration make a contribution Sc to the atomic stress tensor S of the total system given by Sc ˆ ¡

N L N X 1X 1X ri Gi ˆ ri Ak ;i ; ¶k V iˆ1 V kˆ1 iˆ1

…

†

… 5†

where V is the volume of the system. H ence the contribution W c of the constraints to the total virial of the system is W c ˆ ¡V Tr Sc ˆ

N X iˆ1

ri ¢ Gi ˆ ¡

L X kˆ1

¶k

N X

…

iˆ1

†

ri ¢ Ak ;i :

… 6†

P The total virial W ˆ i ri ¢ Fi ‡ W c is related to the atomic pressure p by pV ˆ … 2hK i ‡ hW i† =3, where K is the total kinetic energy. 3. General algorithms In this section we review the treatment of holonomic constraints in the case when the velocity Verlet (VV) algorithm is used to propagate the molecular trajectory in time. The VV algorithm is also customarily employed in the multiple time-step formulation of molecular dynamics [15].

3.1. T he velocity version of the V erlet algorithm The velocity version of the Verlet algorithm consists of two steps in which the atomic coordinates and velocities are separately advanced [18]. F or the coordinates, propagation by a time step D t gives ri… t ‡ D t † ˆ ri… t † ‡ D t r_ i… t † ‡ ˆ rui … t ‡ D t † ‡

D t

D t2 F t ‡ G t … i… † i… †† 2mi

2

2mi

Gi… t † ;

… 7†

where rui , a quantity we shall need later on, is the result of an unconstrained advancement during the time interval … t ; t ‡ D t †. The atomic velocities are then advanced as r_ i… t ‡ D t † ˆ r_ i… t† ‡

D t F t ‡G t ‰ …† i… † 2mi i

‡ Fi… t ‡ D t † ‡ Gi… t ‡ D t †Š

… D 2t† ‡ 2mD t ‰F … t ‡ D

² r_ i t ‡

² r_ ui … t ‡ D t † ‡

i

i

t † ‡ Gi… t ‡ D t †Š

D tG t‡ t … D †; 2mi i

… 8†

where r_ ui … t ‡ D t † is the velocity predicted at time t ‡ D t in the absence of constraining forces during the time interval … t ‡ D t =2;t ‡ D t †. The implementation of the method of constraints that we adopt [2, 17] is to evaluate the Lagrange multipliers ¶ k involved in the constraint force terms through … r† … v† approximate values ® k or ® k which guarantee respectively that the constraints and their time derivatives are satis ed respectively by the advanced positions (7) for the SH AK E step and the advanced velocities (8) for the R ATTLE step. We give some details on these SH AK E and R ATTLE procedures in the next sections where we also introduce the concept of a uni ed convergence criterion for both routines.

158


3.2. T he S H A KE procedure for an arbitrary holonomic constraint The SH AK E iterative procedure consists of many cycles over all L constraints treated individually and sequentially in order to improve progressively the estimate of the atomic positions at time t ‡ D t given their values at time t. We number the cycles c ˆ 1;. . . ;I max . Satisfaction of a constraint k at cycle c, i.e. at step s ˆ … c ¡ 1†L ‡ k of the iteration, implies an update of † the nk running (‘old’) atomic positions frg…ns¡1 involved k … s† in that constraint to the (‘new’) values frgnk . Convergence is obtained despite the fact that the updating of coordinates for constraint k slightly degrades the results already achieved for those other constraints in which the same coordinates are involved. The key step is thus the update of the set of the nk atomic positions implied by a particular constraint k at a particular cycle number. The result of iterative step s may be written [18, 19] as … s† … s¡1† rj ˆ rj ¡

D

… s† t2 ®k

2 mj

Ak ; j ;

j 2 nk ;

… 9†

… s†

where ® k is given by … s† ®k ˆ

… s¡1†

2 ¼k : 2X 1 t … s¡1† D Ak ; j ¢ Ak ;j mj j2n

¶ k… t † ’

… 11†

where N k is an appropriate factor forcing the left-hand side of (11) to be dimensionless and of an order of magnitude of D r… s† =d, where d is a typical bond length and D r… s† is the discrepancy between the running atomic positions and their asymptotic values.

… r† ¶k

ˆ

I X

… s†

®k ;

cˆ0

… 12†

where the sum extends over all SH AK E steps involving constraint k over successive iteration cycles c, i.e. s ˆ cL ‡ k. The Lagrange multipliers given by equation (12) are useful to evaluate the stress and virial contributions according to equation (5) and equation (6), respectively. F or the stress tensor one obtains Sc ˆ

L X 1X … r† rj Ak ; j : ¶k V kˆ1 j2n k

… 13†

3.3. T he RA T T L E algorithm The procedure is similar to that used for the SH AK E algorithm. All time derivatives of the constraints are treated in succession within a cycle which is then repeated until convergence is satis ed. The updated velocity at iterative step s of atom j 2 nk appearing in the constraint k is … s† … s¡1† r_ j ˆ r_ j ¡ … sˆ0†

In (10), the constraint term in the numerator and the constraint gradient Ak ; j with superscript … s ¡ 1† are … s¡1† functions of coordinates rj relating to time … t ‡ D t † while constraint gradient terms without superscript are functions of coordinates evaluated at time t. We emphasize that at iterative step s, any atom i which is not involved by the constraint k is una ected by the … s† … s¡1† change, i.e. ri ˆ ri . The procedure is initiated with the set of values … sˆ0† ˆ ruj , as given in equation (7). The convergence rj of the procedure is quite robust, requiring typically º 20 cycles over all constraints before the latter are simultaneously satis ed within a prescribed relative tolerance ½ º 10¡6 or 10¡7 . Explicitly, the iterative procedure is terminated at the end of the iteration cycle c ˆ I max , i.e. at iterative step s ˆ I max L , if 7 1 7 7 … s† 7 7¼k 7< ½; kˆ1;...;L N k

max

… 10†

k

max

N ote nally that an estimate of the Lagrange multiplier ¶ k… t † appearing in the force of constraint Gj… t † in equation (7) is

where r_ j … s†

D t …s† A ¯ ; 2mj k k ; j

j 2 nk ;

… 14†

ˆ r_ uj and the incremental Lagrange multi-

plier ¯k is given by … s† ¯k

ˆ

2

X

… s¡1†

r_ j

j 2 nk

¢ Ak ;j

7 : D t X 17 7Ak ;j 72 j 2 nk

… 15†

mj

All quantities involved in the R ATTLE step (14) and … s† (15), including ¯k , correspond to time (t ‡ D t). N ote that, for a particular constraint k, the associated Lagrange multiplier is now given by ma x

… v† ¶k … t

‡ D t† ˆ

I X cˆ0

… s†

¯k ;

… 16†

where again s ˆ cL ‡ k. The value of the Lagrange multiplier is numerically close, but not equal, to its … r† counterpart ¶k … t ‡ D t † obtained from the SH AK E procedure at the next integration step. The numerical values are in general di erent due to the discrete nature of the integration algorithm. Convergence is achieved when the time derivatives of all constraints are su ciently close to zero. The same convergence criterion can be used as in the SH AK E routine, namely equation (11) applied to ¼k… t ‡ 2D t †

159

M olecular dynamics of polymers with ex plicit but frozen hydrogens jr212 ¡ d212 j jr12 ¡ d12 j º < ½: d12 2d212

at an iterative step s completing a cycle over all constraints: max

kˆ1 ;... ;L

1 … s† 1 j¼ … t ‡ 2D t†j ˆ max j‰ ¼…ks†… t ‡ D t † Nk k kˆ1;... ;L N k

‡ ¼_ …ks†… t ‡ D t †D t ‡ ¢ ¢ ¢Šj

1 … s† j¼_ k … t ‡ D t †jD t < ½; kˆ1;... ;L N k

º max

… 17†

where we have used the rst-order approximation to the left-hand side and taken into account the fact that the zeroth order term is zero or made su ciently small by the previous SH AK E step. In equation (17), the term … s† ¼_ k is a function of positions and updated velocities r_… s† at time … t ‡ D t †; the last condition can therefore also be written as 7 7 7 1 7 7X … s† 7 max … 18† 7 r_ j ¢ Ak ;j 7D t < ½: 7 kˆ1;...;L N k 7 j2n k

4. Speci c constraints 4.1. T he bond constraint F or a rigid bond of length d12 between two atoms labelled 1 and 2, we write the constraint as ¼b… r12 † ² r212 ¡ d212 ˆ 0;

… s† ®b ˆ

… s†

m1 m2 … r12 †2 ¡ d212 : 2D t 2 m1 ‡ m2 r… s† ¢ r12 1

12

… 20†

The contribution of the bond constraint forces to W c and Sc are, respectively: W c ˆ ¡2¶ b d212 ; Sc ˆ

2 ¶ b r12 r12 ; V

… 21† … 22†

… s†

where ¶b is the sum of the ® b over all cycles (see (12)). F or the R ATTLE algorithm, the corresponding estimate … s† of the parameter ¯b for the bond constraint (see (11)), is … s† ¯b ˆ

… s†

1 m1 m2 r_ 12 ¢ r12 : D t m1 ‡ m2 d212

The corresponding criterion for the R ATTLE step, as given by (18), is 7 1 7 7 … s† 7 r ¢ r _ … 25† 7 12 12 7D t < ½: d212 4.2. T he conic constraint pair We now consider three successive sp 3 carbons of the chain skeleton C i¡1 , C i , C i‡1 , as illustrated in gure 1. The carbons are connected by two rigid bonds of length dCC and characterized by a bending degree of freedom, the angle being controlled by some bending potential. Consider one hydrogen atom linked covalently to the central carbon, C i . The vibrational motion associated with this hydrogen is eliminated by taking it as xed in a local reference frame Ouvw centred on C i and de ned as follows: the u axis is chosen along the bisector of the C–C–C bending angle, the v axis is chosen parallel to the axis going through C i¡1 and C i‡1 nuclei, and the w axis completes a right-handed system of coordinates. If u, v, w are the corresponding unit vectors, the Cartesian coordinates of the four atoms in the laboratory frame are rH ˆ rC i ¡ dCH cos ¬ u ‡ dCH sin ¬w;

… 19†

where r12 ˆ … r1 ¡ r2 †. The constraint index (here ‘b’) will from now on denote the type of constraint (here bond constraint) because all subsequent developments apply to an individual constraint. U se of equation (9) and the appropriate constraint gradient, namely Ñ 1 ¼b ˆ ¡Ñ 2 ¼b ˆ 2r12 , gives an expression for the SH AK E par… s† ameter ® b :

… 23†

F or bond constraints, the SH AK E iteration procedure is assumed to have converged when (see (11))

… 24†

rC i§1 ˆ rC i ‡ dCC cos

…2†u § d

CC

sin

… 26†

…2†v; … 27†

where dCH is the C–H distance and ¬ is the xed (constant) angle between the relevant (C–H ) bond and the bisector of the variable angle . In writing the above expressions we suppose the hydrogen atom lies on the positive side of the w axis. The scalar product of the vector rH ¡ rC i with any vector rCi§1 ¡ rC i satis es the relation

… rH ¡ rCi † ¢ … rC i§1 ¡ rC i † ‡ dCH dCC cos ¬ cos

…2† ˆ 0:

… 28†

If we now replace the term cos … =2† in (28) by the expression cos

…2† ˆ f ‰1 ‡ d 1 2

¡2 CC… rC i‡1

¡ rC i † ¢ … rCi¡1 ¡ rC i †Šg1=2 ;

… 29†

we obtain two relations between the coordinates of the hydrogen atom and of the three carbon atoms. These can be written in compact form as ¡2 r21 ¢ r41 †1=2 ˆ 0; ¼c ² r31 ¢ r21 ‡ T… 1 ‡ dCC

… 30†

160


where T ˆ … dCH dCC cos ¬† =21=2 and the atoms, now relabelled from 1 to 4, correspond to the two ordered sequences seq 1 ² fC i ;C i‡1 ; H ;C i¡1 g and seq 2 ² fC i ;C i¡1 ;H ;C i‡1g relative to two distinct conic constraints. The new constraint scheme used to x the hydrogen atom consists of specifying one bond constraint (the C– H bond) plus two conic constraints relative to the atom sequences seq 1 and seq 2 . Together, these constraints de ne three surfaces whose intersection determines the position of the hydrogen atom. The bond constraint de nes a spherical surface while the two conic constraints de ne two conic surfaces referring to cones with apex at the C i atom, symmetry axis along one C -C bond and a cone angle ³ º 109¯ given by a dependent instantaneous value, namely cos ³ ˆ ¡ cos ¬ cos … =2†. These three constraints in fact yield two solutions which are disposed symmetrically with respect to the C–C–C plane. The SH AK E procedure, being initiated with a good estimate of one particular hydrogen atom position (in our example, the one on the positive side of axis w), will converge to the physical solution automatically. The constraint gradients Ñ j¼c , j ˆ 1;. . . ;4, are Ñ

1 ¼c

Ñ

2 ¼c

Ñ

3 ¼c

Ñ

4 ¼c

ˆ … r12 ‡ r13† ‡ R… r12 ‡ r14† ;

… 31†

ˆ r21 ;

… 33†

ˆ r31 ‡ R r41 ; ˆ R r21 ;

… 32† … 34†

where, for compactness, we have de ned R ˆ … dCH cos ¬† =‰ 4dCC cos … =2†Š. When these expressions are substituted in (9) and (14), the SH AK E and R ATTLE steps speci c to the conic constraints are obtained. Each elementary step provides the parameter … s† … s† ® c or ¯c , and ¶ c is given by their sum over all cycles (see equations (12) or (16)). The contribution of conic constraints to the stress tensor is Sc ˆ

1 ¶ c‰… r12 r13 ‡ r13 r12 † ‡ R… r12 r14 ‡ r14 r12 †Š ; … 35† V

while the contribution to the total virial W , takes the simple form W c ˆ ¡2¶ c… r12 ¢ r13 ‡ R r12 ¢ r14†:

… 36†

F or the convergence of the SH AK E and R ATTLE procedures, we normalize the conic constraints given by (30) by the factor N c ˆ dCC dCH before comparing it to ½ (see equations (11) and (17)).

5.

Tests of the conic constraints in calculations for polypropylene 5.1. Full atomic models of polypropylene To assess the e ects of employing conic constraints and the gain in CPU time resulting from their use in molecular dynamics simulations, we consider two models, based on the same force eld, which di er only in the number of frozen degrees of freedom. The force eld chosen is the PP model developed recently by Antoniadis et al. [12] but adapted to a full atomic description in which methyl groups are treated in terms of individual atoms and where all hydrogen atoms are considered as mass particles with coincident centres of force. M odel A treats the hydrocarbon chain as fully exible except that the C–C and C–H stretching vibrations are eliminated by imposing bond constraints. Energy contributions from the bending terms and torsion potentials, and from both local and non-local nonbonded interactions are summarized in the Appendix. In model B, we eliminate in addition all remaining vibrations involving hydrogen atoms by considering two conic constraints per hydrogen atom covalently linked to a skeletal carbon. At the same time, we maintain full rigidity of all side or end methyl groups, i.e. we eliminate the C–C–H and H –C–H bending vibrations by adding extra bond constraints between second neighbour atoms [8]. In the methyl group case, we add to the list of constraints three C(C)H bond constraints and either three (model B) or two out of the three (model B 0 ) H (C)H bond constraints. In model B, one constraint relative to the methyl group (say one of the H (C)H bond constraints) is redundant because it is automatically satis ed when all other constraints are applied. In our simulations, the redundant constraints are included because they accelerate dramatically the convergence of both the SH AK E [8] and R ATTLE algorithms, as we shall shortly illustrate. A full description of these models is given in the Appendix. The geometry and potential parameters proposed by Antoniadis et al. have been kept unchanged wherever possible. We recall here that their model considers three types of atoms: carbons, hydrogen s and united atoms to model the methyl groups. F or nonbonded interactions, a distinction is made between local and non-local (non-bon ded) interactions. The main modi cation made in the present work is to replace the united atoms by explicit carbons and hydrogen s described by the same force eld as the corresponding explicit atoms in the model of Antoniadis et al. 5.2. T he aPP model and the generation of initial con gurations The test on polypropylene melt was carried out on a system of ve atactic aPP chains

M olecular dynamics of polymers with ex plicit but frozen hydrogens H –(–CH 2 –CH (CH 3 )–)25 –CH 3 at 400 K at the experimental density (under 1 atm) of » ˆ 0:793 g cm ¡3 . The con guration was prepared by a two-stage procedure [20]. F irst, an equilibrium melt con guration is generated by M onte Carlo simulation for a simpli ed bead– rod model where the number of beads is equal to the number of carbons in the PP chain skeleton and the rod length is equal to the C–C distance. This step provides equilibrated conformations of athermal chains in a melt. N ext, the chemical details, including a correction of the bending angles of the skeleton, are incorporated in the model, taking care to satisfy all the constraints of model B. A subset of the desired number of chains is then brought to the experimental density » by enclosing them in a simulation cell of appropriate volume V . This last step leads inevitably to close atom–atom distances of intra- or inter-molecular nature; the chains are in equilibrium on the global scale but not at the local scale. Energy minimization techniques combined with M D runs and thermalization is therefore used to drive the system to a fully equilibrium state, the last part of the equilibration (º 50 ps) being performed on a N V T trajectory using a N ose´–H oover thermostat. Once a con guration of aPP equilibrated with model B is obtained, it is used as a starting con guration for a M D experiment with model A. Because the additional degrees of freedom allowed by model A are not activated in the starting con guration, an additional equilibration run is carried out at 400 K with coupling to a N ose´–H oover thermostat (typically 20 ps are long enough because of the nature of the starting con guration). In summary, the rst stage provides samples of equilibrated chains on a large length scale while the second part of the procedure drives the system to a local equilibrium at the scale of one or two monomers. 5.3. T hermodynamic results In table 1, we list the thermodynamic results obtained in the simulations of aPP with models A and B. Values are listed of the cohesive energy and the pressure, inclusive of the long-range corrections to the van der Waals (Lennard-Jones) interactions beyond the cut-o distance rc ˆ 0:95 nm, calculated by assuming that all atom– atom pair correlation functions have reached their Table 1.

» =kg m ¡3 T =K p =M Pa … E coh =V †1=2 =… J m ¡3†1=2

161

asymptotic limit (explicit expressions are given at the end of the Appendix). Because we want to compare our results with those relevant to the parent force eld, it should be noted that Antoniadis et al. truncate the Lennard-Jones (LJ) interactions at 1:45¼ and prolong them by a quintic spline going smoothly to zero at 2:3¼ for the dynamics. F or the calculation of thermodynamic quantities, they retain LJ contributions if the distance is below 1:45¼ while the long range contributions are calculated as usual with 1:45¼ as the lower limit in the integral. Some di erences between their results and ours must therefore come from the di erences in the treatment of the LJ interactions. We have quoted in the table the values of the H ildebrand’s solubility parameter ¯ ˆ … E coh =V †1=2 , where E coh is the total cohesive energy of intermolecular origin. As the results of table 1 show, there is very good consistency between the results obtained with models A and B. This is a clear indication that the freezing of the hydrogen vibrations by the new constraints does not a ect the thermodynamics of the system, all other aspects of the model being equivalent. We are presently testing whether the dynamics in the pico-to-nanosecond time scale, as probed by time correlation functions pertinent to N M R or neutron scattering experiments, is di erent for models A and B, some extra D ebye– Waller term being expected for model A. R easonable agreement also exists between the earlier results of Antoniadis et al. and the present ones. Our value of H ildebrand’s solubility parameter appears somewhat too large, as it could be expected to increase if the temperature were lowered to the experimental value of 300 K . The model, if tested under 1 atm and 400 K conditions, would lead to a somewhat larger molar volume than the experimental one, but this would decrease the solubility parameter. On the whole, the full atomic model considered here seems to be more repulsive than its united methyl group counterpart. Our aim in this work is not to design a new aPP force eld on the basis of M D simulations. We could, however, move in that direction by combining the aPP results with those for polyisobutylene (PIB), a system which we are currently treating by two models analogous to the aPP models A and B.

Thermodynamic data of the simulated systems.

M odel A

M odel B

Experiment

793 407 42 ¡15:6 £ 103

793 403 42 ¡15:6 £ 103

793.1 400.1 0.1 ¡15:3 £ 103 (at p ˆ 10 M Pa and T ˆ 300 K)

162

J.-P. R yckaert et al. 7

250

6 200

D E / D K (%)

150

4

3

100

2

CPU time (s) / ps trajectory

5

50 1

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0

time step / fs F igure 2. R atio of uctuations R ˆ D E =D K (in % ) (model A (*) and model B (*)) and CPU time per picosecond of trajectory (model A (&) and model B (&)) as a function of the time-step size. Continuous lines are quadratic ts for R (considering D t up to 3 fs only) and the dashed ones are guides to the eye for the CPU time. The horizontal line at R ˆ 1:5% is used to nd the optimal time step size and therefore the e ciency on the ordinate axis on the right.

5.4. T ests on e ciency In this section, we test the M D program based on the use of bond and conic constraints in combination with the velocity Verlet (VV) algorithm. We consider N VE trajectories of our full atomic aPP sample at the state point T ˆ 400 K , » ˆ 0:793 g cm ¡3 , using either model A or model B. The criterion we use to judge the quality of the numerical integration is the ratio R of uctuations of the total energy (due to algorithm and round-o errors) over the (natural) uctuations of the total kinetic energy of the system. A value R º1–2% is generally considered as a good compromise between e ciency and accuracy for performing M D experiments on the nanosecond time scale. To avoid uctuations in the conserved total energy due to the discontinuity of the pair potential at the cut-o rc , we adopt, for this test calculation, LJ potential functions shifted to zero at r ˆ rc . F or each model A or B, the ratio R has been measured for a series of 12 ps runs starting from a unique initial

con guration, but di ering in the size of time step. In these runs, we choose a xed (probably unnecessarily small) tolerance value of 10 ¡8 for the constraint satisfaction. We note that the use of a single value of the tolerance ½ leads in all cases to an equal number of cycles in SH AK E and R ATTLE routines. F igure 2 shows how the ratio R and the CPU time per picosecond of trajectory (in scalar mode on the Compaq AlphaServer G S140 workstation composed of 12 processors of 699 M H z of the Computing Centre of U LB) vary with the time step D t for models A and B. We rst observe for each model the expected R / … D t †2 dependence of the global error associated with the VV algorithm [17]. The prefactor of this power law is larger for model A, which contains those fast hydrogen motions suppressed by the additional constraints used in model B. F or a given time step, model B leads thus, as expected, to a better total energy conservation.

M olecular dynamics of polymers with ex plicit but frozen hydrogens H owever, for any particular time step, the CPU time required is slightly larger for model B than for model A. This is shown in gure 2 which reports, for both models, how the computational cost per picosecond of trajectory decreases when increasing size of the time step. The larger CPU time needed when using model B comes from the larger number of constraints in this model (L ˆ 510 or L ˆ 229 for models A and B, respectively) which must be handled by the iterative SH AK E and R ATTLE methods. To combine fairly these opposite trends in the present e ciency comparison, we adopt R ˆ 1:5% as an acceptable choice for the numerical integration quality. On this basis, gure 2 suggests that model A and model B allow, respectively, a time step of D t ˆ 1:1 and 2:6 fs. In turn, the CPU time per picosecond of trajectory corresponding to these ‘optimal’ time steps are respectively 94 and 54 s for model A and model B, yielding a global gain in e ciency in favour of model B by a factor 1.75. This gain in e ciency is slightly lower than the ratio of the time steps themselves as a result of the constraint treatment. In addition to the larger number of constraints in model B already mentioned, we need for the latter a larger number of SH AK E and R ATTLE cycles to achieve constraint satisfaction (º 33 versus º 15 for model A), because of the larger size of the time step. As the constraint routines were not fully optimized, we believe that the gain factor in e ciency of model B should be close to 2. 0 We have also considered an alternative model B in which redundant constraints (used in the rigidi cation of methyl groups) are not added to the list in the SH AK E and R ATTLE routines. This scheme suits better most M D packages which link the average kinetic energy and the temperature through the number of degrees of freedom in the system, a number calculated as the di erence between three times the number of atoms in the system minus the number of scalar constraints. (It is however easy to modify the package by removing the number of non-redundant constraints instead.) M odel B 0 is fully equivalent to model B and should thus give identical trajectories up to e ects due to di erent errors in the constraint satisfaction and round0 o errors. At rst sight, model B is expected to be more e cient, since fewer constraints must be treated per cycle. H owever, we observe that typically we need ve more cycles to achieve the desired tolerance on the constraint satisfaction, and the CPU time is doubled, 0 making model B and A equivalent in e ciency! We therefore stress again that redundant constraints [8] used in the rigidi cation of methyl end groups are a major advantage.

163

6. Conclusions In this paper, we have introduced a new functional form of a holonomic constraint which is appropriate to freeze the methine and methylene hydrogen atoms in a typical polymeric molecule. The mathematical treatment of this constraint is particularly simple, which allows a rather trivial extension of the SH AK E and R ATTLE routines combined with the velocity Verlet integration algorithm. With such constraints, time steps in the range of 2 to 3 fs can be used in ordinary molecular dynamics simulations of hydrogen ated polymers with explicit consideration of hydrogen atoms with associated masses and centres of force. We have also stressed in the present paper that redundant constraints in methyl groups are imperative to minimize the CPU time spent in the constraint routines because they speed up the convergence. The methodology introduced in this paper appears as a new valuable tool for the study of a wide variety of polymers, including polypropylene, polyisobutylene or polyvinylchloride melts. In particular, for aPP and PIB, the coupling between methyl side groups and main chain motions is of considerable interest in N M R [21, 22] and neutron scattering studies [23, 24]. Support from atomistic simulations on these melts in the nanosecond range should be decisive in the interpretation in molecular terms of the observed relaxations. The method proposed here has been implemented in the D LPR OTEIN package [25], version 2.1, and is available for academic purposes. The idea of the conic constraints came during enlightening discussions with G . D estre´e who is also acknowledged for extensive technical help. We thank I. R . M cD onald for a critical reading of the manuscript. We are grateful to K . K aratasos, F . Saija and D . N . Theodorou, partners within the TM R network ‘N ew routes to understanding polymer materials using experiments and realistic modeling’ under contract ER B-F M R XCT98-0176, for useful comments and help with some technical aspects. G . Arialdi acknowledges the F onds pour la formation a` la R echerche dans l’Industrie et dans l’Agriculture (F RIA) de la Communaute´ F ranc¸aise de Belgique for a PhD grant. S. M elchionna acknowledges support from the Leverhulme Trust. Appendix: The all-atoms polypropylene models We report here some additional details about models A and B of polypropylene used in this work. Both models, which are minimal variants of the Antoniadis et al. model [12], are all-atoms models treating each nucleus as a mass particle with an associated centre of force. All covalent bonds are taken to be rigid with lengths identical to those of the parent model [12]. The potential function used in our H amiltonian is

164 U… frg† ˆ

J.-P. R yckaert et al. X

angles

‡

1 2k ³… ³

X i