Computer simulation package for liquids and solids with polar interactions

Computer Physics Communications 42 (1986) 271—300 North-Holland, Amsterdam

271

COMPUTER SIMULATION PACKAGE FOR LIQUIDS AND SOLIDS WITH POLAR INTERACTIONS I. McMOLDYN / H20: Aqueous systems Aatto LAAKSONEN

*

IBM Corporation, Data Systems Division, Department 48B MS 428, Neighborhood Road~Kingston, NY 12401, USA Received 2 February 1986; in final form 20 June 1986

Title of program: McMOLDYN/H20

Keywords: computer simulations, Monte Carlo, molecular dynamics, aqueous systems

Catalogue number: AALR Program obtainable from: CPC Program Library, Queen’s University of Belfast, N. Ireland (see application form in this issue) Computer: IBM 4341/VAX 750/CRAY iS Operating system: VM or MVS/VMS/COS Programming language: FORTRAN 77 High speed storage needed: dynamical memory allocation by means of work array techniques Peripheral used: line printer No. of lines in combined programs and test deck: 4600

*

Nature of physical problem Study of static and dynamical properties in aqueous systems. Method of solution The water system is mimicked by filling a cubic box with 100 to more than 1000 molecules with desired density, using penodic boundary conditions and minimum mage convention. Water—water interactions are described by ab initio site—site potentials. The Coulombic interactions due to the fractional charges are calculated using the direct Ewald sum technique. Choice can be made between the Metropolis-MC and the conventional MD. In the MD mode, the equations of motion are solved optionally by either using the Verlet leap-frog scheme (rotational by Fincham), the conventional Gear predictor—corrector method or the second order quaternion method by Sonnenschein. The constrained dynamics by Memon et a!. is also included within the leap-frog scheme.

Present address: Dept. of Physical Chemistry, Arrhenius Laboratory, University of Stockholm, S-106 91 Stockholm, Sweden.

OO1O-4655/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

272

A. Laaksonen

/

Computer simulation package for liquids and solids. I

LONG WRITE-UP 1. Introduction McMOLDYN is a family of computer simulation programs for applications ranging from mono-/polyatomic/ionic liquids/ solids to rather arbitrary. mixtures thereof, including biological macromolecules. The main characteristics of the standard versions of these programs is that the interaction potentials are obtained from ab initio calculations, carried out by Clementi’s group [1]. The main goal in the construction of this package has been flexibility from different points of view. First, we wanted to make it specially simple to change the size of the simulated system in terms of particles. We therefore decided to build in a dynamic storage allocation scheme by means of work arrays for the particle vectors. Different simulation schemes in the program (Metropolis-MC and four different algorithms to solve the equation of motion in the molecular dynamics simulation) were included for flexibility, One reason to include Monte Carlo (The “Mc” option in McMOLDYN) is to generate a start configuration for a subsequent molecular dynamics simulation. In the biomolecular McMOLDYN programs it is possible to use a selective MC in specified regions of the system in order to obtain a faster energetic relaxation. The Monte Carlo mode in McMOLDYN can be also used in a simulation of a canonical ensemble in order to obtain static properties. We recommend the leap-frog algorithm or the constrained dynamics to obtain long time step in order to cover the longest possible time scale for a conventional MD method. The leap-frog scheme from Verlet [2] is successful for monoatomic partides and even for small molecules after it has been extended to the rotational motion of rigid polyatomic molecules quite recently by Fincham [3]. The Fincham algorithm is incorporated in our program. In the simulations of chain molecules and macromolecular systems the constrained dynamics has been frequently used in literature. The method developed by Ryckaert et al. [4] is useful for smaller molecules while the procedure “SHAKE” by Van Gunsteren and Berendsen [4,5]

is mainly used for chain- and macromolecules. In this program for the water dynamics we use the constrained dynamics scheme by Memon et a!. [6] which is a further development of the algorithm by Ryckaert et al. [4]. In the applications where an accurate integrator is needed we provide two higher order algorithm to solve the equations of motion: the predictor—corrector method by Nordsieck and Gear [7] or the second order quaternion method (for the rotational motion) by Sonnenschein [8]. For sufficiently small time step the Gear predictor corrector scheme gives a very good energy conservation while its stability decreases rapidly when the magnitude of the time step is increased [9]. The second order quaternion method by Sonnenschein is an elegant way to solve the rotational equation of motion. In this algorithm the first order equation for the quaternions and the first order equations for the angular velocities (as used in the standard predictor—corrector scheme) are replaced by a single second order equation for the quaternions. Compared to the Gear method. Sonnenschein’s discovery improves considerably the numerical stability in the solution of the rotational equation of motion. The method of Sonnenschein [8] is also general, fully capable of handling both linear and non-linear molecules in an efficient way. The predictor—corrector schemes in McMOLDYN, can optionally be chosen from third to seventh order. Using McMOLDYN programs the user is allowed to switch between the MC and the MD modes as well as to change integrators during an equilibration phase in order to effectively minimize the required computer time for a complete simulation. Another aspect concerning the flexibility of the program is found in its structure. In order to create a conceptually clear program for statistical mechanical computer simulations, we have isolated all the basic operations in the Metropolis MC and the conventional MD techniques in separate subroutines. The program is therefore very easy both to learn and to modify by the user. This first paper describes a rather general simu-

A. Laaksonen

/

Computer simulation package for liquids and solids. 1

273

lator for water systems, which is built from more than 50 subroutines. Although this program is specially designed for computer simulations of water, it is easily modified for other homogeneous systems. Only few routines have to be changed as will be explained later in the text. One advantage of the relatively high number of subroutines is that the total size of the program, where both MC and MD with four different integrators are included, is effectively minimized, because the same fundamental routines can be used in different schemes. It is difficult to construct a general computer simulation program in a “black box” manner, but it is possible to minimize the number of user supplied routines, when going from one molecular system to another, by maximizing the number of “utility-type” subroutines,

to a same analytical expression but based on different quantum-mechanical calculations. The MCY potential is obtained using contracted Gaussian (11, 7, 1/6, 1)<4, 3, 1/2, 1> basis set at the CI level from 66 different water dimer configurations. The potential of Carravetta and Clementi [13] is based on quantum-mechanical calculations where (11, 7, 2/6, 1)(4, 3, 2/2, 1) Gaussian basis set was employed. The number of water—water orientations in the Carravetta—Clementi potential is 105. The electron correlation effects are calculated using the approximative functional technique by Colic and Saivetti [14]. The potential is a three-point-charge model, where two positive charges of magnitude + q are placed at the sites of the hydrogens and one negative charge, 2q is lying on the C2,, axis at a

general enough to serve different physical models, There is a long tradition at our laboratory to use ab initio potential functions [1]. We are convinced that it is more natural to use interaction potentials, accurately calculated at the molecular level when used in simulations based on molecular models. The major part of the available ab initio potentials are calculated for the dimers of molecules, thus neglecting the simultaneous interactions between three, four, particles; For examplc, in the case of water, as much as 10—15% has been shown to come from the many-body interactions to the configurational energy of the system [1]. Our aim is initially to calculate quantitatively good pair-potentials from reliable quantum mechanical calculations and later include higher order terms in the potential energy expansion. For water there are now theoretical three-body [10], and even four-body [11] intermolecular potentials available.

distance R from the oxygen towards the hydrogens. The magnitudes of q and R for the both potentials are given in table 1. The analytical form of the pair-potential is written using the original notation from the paper by Matsuoka et al. [12]

—

u

=

qq

-k-- +

2.1. The potential model This standard version of McMOLDYN is supplied with two optional ab initio water pair-potentials, the MCY configuration interaction potential [12] and a recent potential by Carravetta and Ciementi [13]. The interaction potentials are fitted

2qq

+

r14

1

...

2. Method

—~---

T13

—

18

+

r23

1 +

—

+

r24

1 +

28

—

37

r78

1 +

B

+A1 e

—

1r56

47

Table 1 Specification of the parameters ~ for MCY and CaC pair potentials (all values in atomic units)

_____________________________________________ q R d) A1 A2 A3 A4 B1 b2 B3 B4

MCY b) 0.7175 0.5059 1734.1960 1.061887 2.319395 0.426006 2.726696 1.460975 1.567367 1.181792

CaC ° 0.6579 0.47231 724.1307 5.7051 3.3703 0.7297 2.5165 2.0345 1.6808 1.3127

~ The analytic expression for the potential is given in eq. (1). c) Ref d) The distance between the oxygen nucleus and the point charge.

274

A. Laaksonen

/

Computer simulation package for liquids and solids. I

+A2 e82n13 + A2 eB2r,4 + A2 e82r23 +A 2 e_~2r24+ A3 e~3’~16+ A3 e_B~26 +A 3 eBir3s + A3 e~3r45+ A4 eB4n6

where N is the number of particles and NA ~ Avogadro’s3)number, or from the specified density p (in g/cm

3,

+A4 eB4r26 +A4 eB4r35 +A4 eB4r45. (i) The water model is shown in fig. 1. The hydrogens in the first water i are denoted by 1 and 2 while they in the second water j are labelled as 3 and 4 Similarly 5 and 7 are the oxygen and the massless point charge in the first water molecule whereas they are labelled as 6 and 8 in the second water. The values of the fitting parameters A and B are quoted in table 1 for both potential functions,

(3)

L = (NmHO/10Oop)~ where m~ 20 is the mass of the water molecule. In a standard manner the side length L is scaled to 2 in the internal dimensionless program units. Now, placing the simulation box in the Cartesian coordinates so that the corners are at (+ / 1, + / 1, + / I) the periodic boundary conditions are given in a convenient way through the coordinate transformation: p=p—2INT(p), (4) —

—

—

2.2. Computer simulation method A cubic cell is filled with water molecules (normally arranged in a crystal, like fcc or simple cubic lattice). The side length L (l~= I, = L) is determined from the specified molar volume Vm: L

=

( VmN/NA )1/3

(2)

5

2 7

~0

14

where p = x, y, z and INT(p) is the truncated integer part of the centre of mass vector p. By fixing the cut-off radius for the pairwise interactions at L/2 the same transformation can be used to produce the nearest image interactions, where p is now a component of the centre of mass separation between two particles in the box. How justified is the use of L/2 as a cut-off radius for the non-Coulombic interactions depends, of course, on the number of particles in the box as well as on the density of the system. An estimate of the long-range correction for the non-Coulombic interactions (i.e. the exponential terms in eq. (1)) can be made by assuming a uniform number density beyond the cut-off radius. These long-range corrections to both the potential energy U and the molecular virial e are given as: 2 f~ ULRC j r,~U(r,~) dr,~ (5) non-Coul 2~iN L3 L/2 — —

1

and

-2o 8

3

enon~CouI = LRC

6

Fig. 1. Distribution of the fractional charges in the potential model used in McMOLDYN/H20.

2nN2

L3 j~d’

d

dr,~.

(6)

We have noticed from our calculations that these long-range corrections to the non-Coulombic potential and virial are always vanishingly small. The Coulombic terms, however, are long ranged. We have chosen to handle these interactions using the direct Ewald sum method [15].As known, the

A. Laaksonen

/

Computer simulation package for liquids and solids. I

Ewald summation was originally developed to calculate the lattice energy in ionic crystals [16]. Now, because the introduction of the periodic boundary conditions, gives a “crystal-like” structure for the simulated system, the Ewald summation can be applied in the system with interacting (fractional) charges. There is, of course, no such a periodicity in real liquids which has been the main source for criticism against the use of Ewald techniques in the computer simulations of ionic liquids. To date however, the Ewald sum techniques are probably among the most reliable alternatives to treat the Coulombic interactions in the simulations of ionic or partly ionic systems. Other proposed methods are the so called reaction field [17] and the multipole techniques [18] to solve the problem caused by the long-ranged Coulombic interactions. The Ewald summation techniques are well established and in the case polyatomic molecules, routinely divided to three parts eqs. (7), (8) and (9). The constant self-interaction term due to the intramolecular correction [19] (in SI units): —

1

all sites / s~i

=

+ __L

q

1

—

~

)

~

~

s,=1 s,—1

2a

I sites on sj~= i

~s,s

0

q~q5 LERFC(a,.), ~

~,

}

=1

(9)

where k is the reciprocal vector. The Coulombic energy is then: UCoul

=

U~oul+

~

+

uc0~ III

(10)

The total configurational potential energy is found by summing up the exponential terms of eq. (1) with eqs. (5) and (10) U=

Unon-C0w + 1~LRC r,non-Coul + UC0~~.

(11)

In Monte Carlo, a set of configurations, Nconf, is generated after weighting with the Boltzmann factor e U/kT The values of the desired quantities are obtained from the averages: 1

Nconf

(12)

1

In molecular dynamics, the classical equations of motion are solved numerically to determine the trajectories for the molecules and the desired quantities are obtained as time averages.

—

(r5)}

~ 411C

—i—— ~

2 U1~°°’ = L ~o k 1 e~/~” /k2{ all sites q e~~”)2

a,.

(7)

where s, and s~denote denote sites on i and j, respectively, and a is the Ewald convergence parameter. ô,~is the Kronecker delta, The slowly convergent lattice sum is transformed to two more complicated bur rapidly converging sums [15]. The first one is calculated in real space.

U~°”l= ~

space:

Nconf

ERF

2q5 ~•5,5~

275

(8)

ERF and ERFC in eqs. (7) and (8) are the error function and the complementary error function, respectively. The second term is calculated in reciprocal

Computationally, the central task in an MD simulation is the evaluation of the forces acting on the molecules due to the interactions of the remaining molecules in the system. The pairwise forces are given as j’.. = ~j. (13) —

or as in our case when the interactions are spatially spherical (i.e. central forces):

~f-,

F~,= U1.,. (14) The potential given in eq. (1) is the so called —

site—site potential and the potential energy between two molecules can be written as Nsites Nsites

U,~=

~ s, = 1

=

U~jr53),

(is)

276

/

A. Laaksonen

Computer simulation package for liquids and solids. I

0, H1, H2 and point charge for our water molecule and r,., is the site—site separation while U is the specifk site—site interaction. ‘the force, acting on the molecule is given now by: s.

=

Nsites

F. =

F =

(16)

i

Nsites Nsites

~

=

i,.

(17)

=1

=

And the pairwise site-site force i,

is

given as

The equations of motion to be solved for each one of the molecules are the Newton—Eulerian equations mv,

=

F,., =

(19) ~

+ ~P)~(P)(i~~

—

i~),

Ns,tes

“~~°“~ = ~

(22)

and transform it to the principal axis system by eq. (21) so it can be used in eq. (20). The conventional way to solve the rotational equation of motion has been to use the principal axis system on the molecule so that the moment of inertia tensor becomes diagonal and to use Eulerian angles to describe the orientations of molecules, Unfortunately, the numerical methods to solve the equation suffer from numerical instabilities for certain orientations and especially for linear molecules. Therefore the so called quaternions are used to describe the molecular reorientations. The quaternions (i,, ii,, ~, x,) span a fourdimensional space and are defined in terms of the Eulerian angles (0,, ~b1,4,) through the following transformation:

(20)

where I~, I~ and I.~are the elements of the diagonalized moment of inertia tensor and is and w are the linear and angular velocities, respectively. F is the torque and (p) denotes principal axis system, (dot) the time derivative and a, $ and y are the x, y or z components. The y- and z-components of eq. (20) are obtained by cyclic permutations of eq. (20). The water molecule is in eq. (20) treated as a simple asymmetric rotor. An

x F5

~,),

(23)

~,),

(24)

~,= cos(0,/2) sin(i.~,+ ~,),

(25) (26)

x,

=

sin(0,/2) sin(4’,

=

sin(0,/2) cos(ip,

=

—

—

cos(0,/2) cos(ip1 +

~.

‘

I

There are several conventions to define the quaternions (21), eqs. (23)—(26) above represent the so called Y-convention. The rotation matrix R 1S given in terms of quaternions as

2(~ 1x1—~~n~)

—~+~—~+~

~2(~

—

2+x~

2(i~4~~+~1x1) 2(~1~+~)

(27)

.

~-n~-~

E 1x1)

—

2(~,q1+ mxa)

arbitrary vector quantity, A ~ is given in the boxoriented coordinate system, is transformed to the principal axis system by an operation with the rotation matrix Ron it.

—~

—

+

+

The first order differential equation for the quaternions is given as

\

(~,

I —~,— x, + m + ~,\I~~\

I~ J ~J I

I.

I

2~

A(P)=RA~0~.

(21)

Now, for example, if we denote the vector directed from the centre mass of the particle i to a site s, as ‘Si we can calculate torque F in the box-oriented axis system as

+xj_~_~+~P) + + + X~+ ~ —n, — ~,— ~, + x

(28)

I

0

1

This equation has a central role in the integration

A. Laaksonen

/

Computer simulation package for liquids and solids. I

of the rotational equation of motion. The following relations hold for the quaternions (orthonormality): + + + +~ +

i.~

=

t, + k,x,

=

(29) (30)

0.

Several quantities (static MC, static and dynamical — MD) can be calculated using this simulation program. The molecular virial 6’ is defined as —

N =

0non-Coul

i(
—

N

~

2.3. Monte Carlo mode Monte Carlo simulations are performed using the Metropolis sampling technique [22]. A flow chart of the MC simulation is given in fig. 2. The R, parameter in fig. 2 denotes the molecular centre-of-mass coordinates together with quaternionThe parameters. MC run is carried out by specifying the -

number of MICRO moves in each MACRO move. After each MACRO move the internal energy of the simulated system is printed out. The number

N

~ ~

+

277

(U( r))

of MACRO moves is called GIANT moves. After each GIANT move a detailed information of the

‘J

N ~r

1~F,1.

~

START

(31)

i(
the virial is then used to calculate the pressure P, of the system from: PV=NkT—~6’,

~oid1d u~

MOVE

N =

i1, 2,..

N

0(5oid)

(32)

where V is the volume and k is the Boltzmann constant. Statistical mechanically the MD method is essentially an ENV ensemble while the MC method an NVT ensemble. temperature in the MD issimulation is obtainedThe from the translational and rotational kinetic energies, Etr and Erot Etr

0

=

=

1,NCONF

N

5’= s~ = u”~

+

R

T

~

YES

mv~v 1,

(33)

j=1

N Erot

~

=

(I~w~w~+

i=1

14)/kT

YES

+I>~w~)

(34)

e _150eW ‘1

u?

=

5new 1

by the relation:

=

1 0)goid) 50014

3NkT=Etr+Erot.

=

ASUM+A)R0~)

(35)

The total energy of the system is given as E

=

(36)

U + Etr + Erot.

Also the pair correlation functions are calculated for the radial distributions of the 0—0, 0—H and H—H pairs. g1(r)

=

n~(r— r

+

dr)/(4rrr dr)(N1/V).

(37)

=

~o.i

=

random number [0.13

ASUM

/ NCONF

STOP

Fig. 2. Flow chart for the Metropolis Monte Carlo used in McMOLDYN/H20.

278

A. Laaksonen

/

Computer simulation packagefor liquids and solids. I

energies is given together with different statistics (i.e. the number of accepted/rejected moves or hard core collisions). The total number of MC configurations will be the number of MICRO x MACRO X GIANT moves.

Main part jn—l/2

2

jn+i/

R~’~’~2 + j;~/2,

=

~n+i/2(P)p.i/2(PyJP

2.4. Molecular dynamics mode

“~-

(48) (45)

+ ~tF,.”,

a=x,

y,

z,

(50)

1

q 2.4.1. Leap-frog integrator Verlet’s algorithm [21can

r, (st)

=

—

r, ( — ~t)

+

1 =q1 ~2 Store jn+i/ momentum.

be given as

2 r,. (0) +

i~t2j.

(0)

+ ~ ( st”),

(38) ~(0)=(r 1(~t)—r1(—L\t))/2

2), (39)

z~t+O(i~t

where i~t is the given time step. Another, and a more convenient notation for eqs. (38) and (39) is

(40)

= 2i” — r”1 + a” L~t2, v7 = (r~”~’ — rn~)/2 itt.

and qfl+i/2,

(51) J is the angular

2.4.2. The Gear predictor—corrector integrator The higher order predictor—corrector schemes are based on the Taylor expansion of the basic equations. Following Nordsieck and Gear (7) we predict i) the centre-of-mass coordinates: —

PRED ~k

—

(41)

f=7 j!/(k!(j

~

—

k)!)rk,

(52)

j=k

where In eqs. (40) and (41) the higher order terms are neglected. The omission of the second and higher order terms in the velocity expression may sometimes be critical. Eqs. (40) and (41) may be written in the leap-frog form: =

=

v7~”2+

,“

(42)

a~’/~t,

(43)

+~

The mid-velocities in eqs. (42) and (43) are more accurate that the corresponding velocities in the

= (z~t”/k!)(a”r/at”) (53) and k = 0, 1, 2, 3, 4, 5, 6, 7 in the scheme used in this program. The coordinates are then corrected by

rk

CORR

rk

—

PRED

—rk

+

(54)

~

where s, = ~

1120) —

PRED

r

2

Verlet algorithm. The leap-frog algorithm is recently extended to the rotational motion of rigid polyatomic moleëules by Fincham [3], giving the following algorithm:

(55)

and c~ is a specific corrector constant. In the next predictor step, the how corrected coordinates ~ are used as rk in eq. (52). ii) similarly the angular velocities:

Auxiliary part

7

(44)

PRED

— ~j!/(k!(j—k)!)wk,

(45) =

n+1/2

Use

J/JP), —

a

=

q7 + ~z~tQ”

qfl+I/2

x, y, z,

(46)

~,~n+l/2(p)

(47)

to calculate R~’~2and

Qn±l/2

j=

(56)

J=k

where =

(~tk/k!)(a~/atk)

(57)

and k = 0, 1, 2, 3, 4, 5, 6, 7, the angular velocities

A. Laaksonen

/

Computer simulation package for liquids and solids. I

are corrected by (58)

w~0RR = W~RED + ~

279

forces, torques and the 0 matrix (given in eq. (28)) are needed to correct the coordinates, angular velocities and quaternions, respectively.

where 2.4.3. Second order quaternton method =

and iii)

2/2!)(TtD/I)

(~t c~

— W~IRED

(59)

is a specific corrector constant,

also the quaternions:

q~RED=

f=7

In the second order quaternion method by Sonnenschein [9] the translational motion is treated in a normal way through the eqs. (52)—(55). In the treatment of the rotational motion, the angular velocity can be found from eq. (28) as:

(60)

~j!/(k!(j—k)!)q~

j=k

(p)

and k = 0, 1, 2, 3, 4, 5, 6, 7 in the scheme used in this program. The quaternions are then corrected by CORR

qk

PRED

(1)

+ ck

= qk

i

—

‘i”

~oz,

0

\

,~61)

q,~,

(63

_2AT

~

.

.

and its time denvative

t~,

where — q,—

2’~”~O

(~i\

PRED

q

1 (1)

t~ui) .

and ck is a specific corrector constant. The corrector constants c~ and cj~ [8] are given in table 2. Notice also that the predicted

Xl

.(p) ~ .

(i,)

=2(AT4,+A~.). ‘

(64)

Wz

0

Table 2 Corrector constants C, used in the predictor—corrector algorithm First order equation1 C~’ 0k 3/8 1 2 3 4

C~’~ 251/720 1 11/12 1/3 1/24 0

C~’1 95/288 1 25/24 35/72 5/48 1/120

C~11 19087/60480 1 137/120 5/8 17/96 1/40

5

1 3/4 1/6 0 0

6

0

0

0

1/720

7

0

0

0

0

C~21

C~2)

Second order equation k C~21 0 1 2 3 4 5 6 7

1/6

5/6 1 1/3 0 0 0 0

19/120 3/4 1 1/2 1/12 0 0 0

3/16 251/360 1 11/18 1/6 1/60 0 0

C~21 863/6048 665/1008 1 25/36 35/144 1/24 1/360 0

C~’~ 5257/17280 1 49/40 203/270 49/192 7/177 7/1440 1/5040

C~21 275/2016 19087/30240 1 137/180 5/16 17/240 1/120 1/2520

280

A. Laaksonen

/

Computer simulation package for liquids and solids. I

Now, inserting eq. (64) in eq. (20) gives after some

(

rearrangement. ((FX + w~a~1~— =

(l)~W~ ( ‘zz — ‘XX ‘XX

\

)

bk,=

1~ )/I~~

(~ +

(r~+ c’xw~(

[6] are expressed as:

~3’~

=

0.

(67)

Memon et al. [6] showed that taking the square root of eq. (66) gives a significant improvement to the convergence rate in the iterative part of the algorithm. They also developed a leap-frog formalism for their constrained dynamics scheme: vfl+l/2 = ~ a7/z~t — r,”~1/.~t, (68)

(65)

~,

~I~rk(t)— r,(t))2

— ivV))/Izz

q / The angular velocities are obtained from eq. (63). The matrix A is defined by eq. (28). —

(69) 2.4.4. Constrained dynamics Using the constrained water dynamics the simulation cell is filled with hydrogen and oxygen atoms, rather than rigid water molecules. In the simulation of pure water the ratio between the number of oxygens and the number of hydrogens should be, of course: 1: 2. The oxygen—hydrogen and hydrogen—hydrogen bond lengths are given as input. After the force calculation between the particles the bond length constraints are used to form rigid triangles corresponding water molecules. Ryckaert et al. [4] presented a scheme for constrained dynamics and applied the method to n-alkanes. The constraints of Ryckaert et al. can beexpressed b~,=(rk(t)—rl(t))2=0,

Following Memon et a!. and writing the constraints C in eq. (67) as C r~”~1 +~ — rf”~1— — bk,— 0, (70)

I

=

where the prime denotes an unconstrained position, In the case of water, the constraints are needed only for the bond lengths. In large and flexible systems, constraints have to be used even for bond angles, dihedral bond angles, etc. Memon et al. Taylor-expanded eq. (70) and expressed it in matrix form: S=Ag+O(g2).

(66)

Forwaterwehave:

iJ r~

where bk, is the bond length between atom k and atom 1 and rk ( I) and r,( I) are the corresponding position vectors between the atom k and the atom

n+I

~

=

r~

—

I

rH

,,+i

1\

+

0

1

jPoHcP~H~*1

I I

i~

mrOHiroH~+I

mH

1 m0 m~ 1 roHrHH —

n*1

mH

1 m0

1

Hf0113’

~~‘H,Hi~’OH5*I

mH

mH 1 —rowrn

~+

,H2

(72)

boH — 0HH~

n+~

—rH

m (1

boH

— rH

2

1. The constraints in the method of Memon et al.

—

(71)

(

mH 1 + mH)H1H2H1H2~ 1

I

,

(73)

A. Laaksonen

where vector, Tkl

I rk —

g

is a column vector

g

=

I

I

Computer simulation package for liquids and solids. I

denotes a normalized (to unity) position

‘k/

(rk — ‘~)/

=

/

1’/

(74)

1’

goH

(75)

2

gH,H2

and has to be solved iteratively. calculate the corrections ~rk ~ gk,~

—

g is

needed to

(76)

m

~

{

I I

tions for the atoms. These new positions are then treated by interdependent constraints, eq. (71), a procedure where the molecules are “refound” at each cycle. Calculation of quatermons and torques is therefore not needed to reorient the molecules. There is, however, a complication in using the four-site MCY potential model. Because one of the sites is massless, namely the point charge on the C2,, axis, it will cause singularities in the constrained dynamics equations. In order to prevent this problem, the forces acting on the point charge have to be distributed on the remaining three sites (i.e. on the oxygen and the hydrogen). This problem is discussed by Impey et al. [23] in their MD work using the MCY potential inside a

3. Description of the program

I

3.1. General features

CALCULATE

g

The program McMOLDYN/H2O is build on a network of three work arrays: One for the centre-

FROM LINEARIZED

eq.

)7o)

YES

S TO P

NO

Fig.

=

replace the forces acting on the point charge to act on the atoms.

CALCULATE

s-

=

constrained dynamics scheme. should stressed, however, that there is no It unique way be to

I

CALCULATE

=

Note that g,,~, g,~ and ~kk 0. A flow chart for the iterative calculation of g is given in fig. 3. The procedure by Memon et al. [6] is called “ADJUST”. Note that the rigid molecules do not exist in the force calculation when the constrained dynamics is used. The forces are assumed to be on the “separated” atoms giving a new set of posi-

goH,

=

281

+

3. Flow chart for the procedure “ADJUST” by Memon et al. [6] used in the constrained dynamics scheme in Mc MOLDYN/H20.

of-mass vectors, one for the site vectors and one for the look-up tables. Division to the three work arrays was rather natural because the dimensions of these large arrays are multiples of the dimensions of the vectors stored in the arrays. The reason for placing the look-up tables in a work array is that this provides the ability to change easily the number of table entries. The user can, by changing the table size, optimize the numerical accuracy of the stored values. Because a linear interpolation is used in reading the energies, virials and forces from the tables, there is always a risk for a poor numerical accuracy if the tables are too small. The accuracy of the look-up tables can

282

A. Laaksonen

/

Computer simulation package for liquids and solids. I

be improved by using a quadratic (or even a cubic) interpolation or using a non-linear tabulation: for example, creating the tables with respect to 1/r2 instead of r2. However, the use of look-up tables is made optional. For example, when runfling the program on a Cray computer, it is very inefficient to read from the tables because it is not vectorizable, On a scalar machine it is faster to use look-up tables than to compute the exponential terms and square roots for all the energy and force terms between all the sites of all the molecules in the simulation cell, 3.1.1. Program level 1

lattice. The last routine to be called depends on the mode in which the program is running on. If we work on the MC mode, MONTEC (COM,SITE,TABLE) is called, where CYCLES (COM,SITE,TABLE) is called on the MD mode, At the top of the MAIN program, a list is given containing all the COMMON blocks with their current dimensions. Also the total number of local variables used in the program is given. Adding up these figures with the three work array dimensions will give the total number of words needed to allocate the program. All the real variables in the present version of the program are in double precision.

The MAIN routine calls three routines, as shown in fig. 4. The work arrays are named COM, SITE and TABLE. The first routine SYSTEM (TABLE) sets up the system. The routine START (COM) picks up the entreof-mass arrays either from a restart file or from a

3.1.2. Program level 2 The organization of the routine SYSTEM (TABLE) is shown in fig. 5. The first routine to be called is INPUT, where the input parameters are read in and analyzed. The next routine INDEX

I

MAIN LEVEL

BLOCK DATA

SYSTEM(TABLE)

1) LEVEL

___________

2

i~uj

SYSTEM(TABLE)

I)IDEX

~coM~

MONIEC

(CIM,SITETABLE)

CYCLES

(COM,SITE.TABLE)]

EN~”II Fig. 4. Program structure in routine MAIN.

T A B L ES (TABLE)

C I ER N S

[RETURN

Fig. 5. Program structure in routine SYSTEM.

A. Laaksonen

/

Computer simulation packagefor liquids and solids. I

calculates the dimensions of the work arrays COM, SITE and TABLE and constructs the address arrays. When using the program, the MAIN routine and BLOCK DATA should be stored separately from the rest of the program. The first run is done only to check the input data and to obtain the work array dimensions. These dimensions are updated in the MAIN program and the BLOCK DATA. The MAIN program is then re-compiled together with BLOCK DATA and linked together with the load module for the rest of the program. By doing the initial set up of the work array dimensions, described above, no memory space will be wasted. The look-up tables for the Ewald complementary error function energies and forces and nonCoulombic energies vinals and forces are created in TABLES (TABLE). Different constant terms are calculated in CTERMS. The routine START COM) (see fig. 6) calls RESTAR (COM) if the run is a restart case. If the time step (MD) is different in the restart file (i.e. it has been changed from the previous run) all the time derivatives read from the restart file are scaled. Also if temperatures are changed, all the temperature dependent quantities coming from the previous run are scaled. For a completely new run, the positions are taken from a simple cubic lattice, and the quaternions, linear and angular velocities are initiated by randomizing because the integrator algorithms are not self-starting. The randomized quaternions are then normalized to unity (eq. (26)) and the velocities are scaled in order to remove the excess heat caused by the randomization. Initiation of the quaternions and velocities are done in QRAND(COM), NAIL(COM) and SCREW(COM), respectively. It should be stressed that the initiation of the quatermons and the velocities as described above is not normally sufficient to start the predictor—corrector integrator, The procedure used in this program to start a new predictor—corrector run is described later in this section. Working on the MC mode we do not need the velocities (or other time derivatives), because it is a “timeless” procedure where the particles are translated and rotated in a stochastic manner. The

S IA R I ( LEVEL

__________

283

CON)

2

___________

LRES

________________

R E S I A B (coN) ________________

.NOT.LRES

.NOT.LRES

L A I I I C

(coM)

U R A N D (cos) _______________

I

.NOT.LRES

N A I L (coM)

,NOT,LRES

S CR E W (caM)

R E IU R N

Fig.

6.

Program structure in routine START.

routine INDEX (called from SYSTEM (TABLE)) calculates automatically the required work array space for all the methods available in the program (i.e. the Monte Carlo and the four different schemes to solve the equations of motion for rigid molecules). The work arrays for an MC run are much smaller than, for example, for an MD run with the leap-frog algorithm. The predictor—corrector method requires most core space of all the five methods, because of all the higher order time derivatives of the positions, angular velocities and the quaternions. 3.1.2.1. MC mode. Working on the MC mode, the program has the routine MONTEC(COM,SITE, TABLE) as the organizing block. The structure of MONTEC is shown in fig. 7. The first routine to be called from MONTEC is UPDATE

284

I

A. I.aaksonen

MONTEC (

/

Computer simulation package for liquids and solids. I

COM,SITE,TABLE

LEVEL

2) __________

~E(COM.S

TIE)

~(O(COM~SITE,TABLE)

~(COM,SITE)

KIII~Ii ~

molecule, also randomly specified. If the Ewald summation is used the reciprocal part of the lattice sum is calculated in BROKE(COM,SITE). The new energy is calculated for the molecule and returned to MONTEC, where it is compared with the old energy for the same molecule. The move is accepted if the new energy is lower than the old energy. The move is not rejected if the energy is higher unless the Boltzrnann factor for the simulation temperature is less than a random number (0—1). This is the Metropolis [23] sampling technique. 3.1.2.2.

~E(CON.S

ITT)

_~EllBR0KE(COMSTTE) __________________

I

R C I U B N

______________

Fig. 7. Program structure in routine MONTEC.

(COM,SITE). UPDATE places one water molecule on each of the centre-of-mass positions in the simulation cell and rotates the molecules by the quaternions producing a configuration for an energy calculation After updating the positions the routine CASINO (COM,SITE,TABLE) is called. CASINO calculates the energy of one molecule due to the other molecules in the simulated systern. Before starting the MC sampling, CASINO is called once for all the molecules, in order to have the “old” energies for all the particles. The routine CASINO calculates both the non-Coulombic and the Coulombic energy. The direct Ewald summation can be optionally used to calculate the Coulombic energy. Inside the MC loop, CASINO calls the routine ROULET with a randomly chosen molecule and suggests a new position to the

MD

mode.

The

routine

CYCLES

(COM,SITE,TABLE) organizes the MD cycling in all of the four MD methods available in McMOLDYN/H2O. i) Leap-frog scheme (fig. 8) Using the leap-frog integrator, the first routine to be called from CYCLES is UPDATE(COM,SITE). As in the case of MC simulations, UPDATE calculates the orientations of the molecules from the quaternions and places the molecules at their centre-of-mass positions. Having the configuration of the molecules, the routine ACTION(COM,SITE,TABLE) is called to calculate the non-Coulombic energies, virials and forces, and the corresponding complementary error function quantities to the Ewald summation. Optionally, the radial distribution function accumulators are updated in ACTION. The reciprocal part of the Ewald sum is calculated in FURIER(COM,SITE). The routine FURIER is taken from ref. [25] and modified for the work arrays. The Fourier series is computed explicitly at each time step. After FURIER the centre-of-mass forces and torques are calculated in COMFT(COM,SITE). Having obtained the centre-of-mass forces and torques, we are in a position to use the leap-frog algorithm to calculate the linear velocities and the new positions in ALONG(COM,SITE). The angular velocities and quaternions are computed in the routine ABOUT(COM,SITE). The translational and rotational kinetic energies are also calculated in ALONG and ABOUT. In the next routine, VALUES, all the therrnodynarnical quantities are

A. Laaksonen

/

Computer simulation package for liquids and solids. I

~ CYCLES (COM,SITE,TABLE) (LEVEL2)

~S~cOM.STTE,TARLE~ LEVEL 2

Mi

-

285

I I

~LRES

U PDATE

CaM SITE

~LRES

ACT ION

COM 5 lIE, TAB LE)

LOOP

~N(cOM,SITE,TABLE~

~LRESPCINIT(COMSITE)

~R)cOM,SlTE)

~(coM,sITE)

~E(coM,sITE)

~G(COM,SITE)

~N(COM,SITE,TABLE)

~(COM,SITE)

~(coM.sITE)

(cSM, SITE)

~C(coM)

~ET(COM)

_ I H

CALAVS

C A LA VS RETURN

LSAVE

RESTAR

I

RESTAR (coM)

RETURN _____________________ Fig. 9. Program structure in routine CYCLES — predictor—cor-

(caM)

Fig. 8. Program structure in routine CYCLES scheme.

~LLSAVE

I

—

leap-frog

computed and added to their accumulators (calculation of the average values is done at the end of the run). The last routine SCALET(COM) is called if temperature scaling of the velocities is desired. This is normally the case during an equilibration procedure. ii) The Gearpredictor—corrector scheme (fig. 9) Initiation of a completely new run using the predictor—corrector algorithm may be difficult without a special operation, because in most of the cases it is not sufficient to initiate only the linear and angular velocities and the quaternions only. Normally it is, however, enough to try to determine the first derivatives of the velocities and quaternions. To initiate a new predictor—correctur run we first call UPDATE(COM,SITE) to obtain

rector scheme.

the configuration of the molecular positions and then ACTION(COM,SITE,TABLE) and COMFT(COM,SITE) to obtain the centre-of-mass forces and torques. The next step is to call PCINIT(COM) where the second time derivatives of the positions are obtained from the forces. Using the randomized quaternions and eq. (28), also the first time derivatives to the quaternions are obtained. Finally, the time derivatives of the angular velocities are calculated from the torques. No time derivatives of the angular velocities are needed in the second order quatermon method. The first routine to be called in a predictor—corrector cycle is PREDIC(COM). In PREDIC all the zeroth and higher order centre-of-mass quantities are predicted. With the predicted coordinates, the routines

286

A. Laaksonen

/

Computer simulation package for liquids and solids. I

ACTION(COM,SITE,TABLE), FURIER(COM, SITE) and COMFT(COM,SITE) are called. Haying the predicted centre-of-mass forces and torques, all these quantities are corrected in CORREC(COM). For the quaternions, the first order ordinary differential equation (eq. (28)) has to be solved in order to correct the higher time derivatives. Again, the values of the thermodynamical quantities are added to their accumulators in VALUES and if temperature scaling of the velocities is required SCALET(COM) is called. This completes the predictor—corrector cycle. If the time-step counter exceeds the desired number of time steps set for the run, the routine CALAVS is called to calculate the average values from the accumulators. This is done routinely in all of the available MD schemes. Also, before closing the run, the restart file is dumped. The first information written on the restart file is a label specifying the method used in the run. This information is needed when the restart file is picked up next time either to continue the same calculation or to serve as an initial configuration for a new run using for the restart file according to a special hierarchy. some other method. The label also gives a position At the top of this hierarchy is a restart file dumped from a predictor—corrector (both the Gear- and the second order quaternion scheme) run. It can be used in the other three methods without any special treatment to the file. The predictor—corrector restart file contains all the information needed to restart a run either using MC, leap-frog algorithm or constrained dynamics. It also contains some information not needed in these lower order methods and those quantities will be skipped. The next method in the hierarchy is the leap-frog scheme, followed by the constrained dynamics and at the lowest level is the restart file generated after a Monte Carlo simulation. It is therefore always straightforward to read a restart file of a higher rank. If the restart file, however, belongs to a lower level than the method aimed to be used in the current run, all the quantities not found in the restart file will be generated by randomization. iii) Constrained dynamics (fig. 10)

The constrained dynamics scheme is rather different from the rigid molecule leap-frog or predic-

CYCLES (COM,SITE,TASLE)

II I

LEVEL 2 )

I I—I.NOI~RE5i—.- —-J

____________________

UPDATE (COM.SITE,TABLE)

1

___________

~

_____

~~(COM.SIJE.TABLE) _______________[R(COM~T~~BLE7~I

ALONE (coN SITE

~T~coM(

_______

RETURN

________

_______

i

I

[C A L A V S

LSAV~1—.~ RESTAR (ToM) ___________________

Fig. 10. Program structure in routine CYCLES dynamics scheme.

—

constrained

tor—corrector schemes. From the molecular point of view, the COM work array will lose its original function, because this time we do not have any centre-of-mass quantities, only the interacting atoms. From the atomic point of view, the SITE work array is not needed for its original purpose. Therefore there is no such a strict division between the contents of the COM and SITE work arrays when the constrained dynamics is used. If it is a new run, the first routine to be called from CYCLES is UPDATE(COM,SITE). In other words the COM and SITE arrays are used as before. This is because both ACTION and FURIER can be called, in the same fashion as in the case of leap-frog or predictor—corrector algorithms. If it is a new run (i.e. started either from

A. Laaksonen

/

Computer simulation package for liquids and solids. 1

lattice or from a restart file generated by another method) even the quaternions are used to start the run. The quaternion arrays are laster used as dummy space for other purposes. To start the MD cycles, the first routine to be called is REORDE(COM,SITE). REORDE calculates a position for the centre-of-mass of each of the “molecules” and also the site vectors from the centre-of-mass to oxygen, hydrogens and the point charge. The next routine called is ACTION(COM,SITE,TABE), to calculate the energies, virials and the forces. After FURIER an additional routine, ATOMF(COM,SITE) is called to distribute the forces, acting on the massless point-charge to act on the oxygen and on the hydrogens. We also call ALONG(COM) to calculate the unconstrained velocities and positions using the leap-frog scheme. The molecules are “refound” after calling the routine ADJUST(COM) to fix the oxygen—hydrogen and hydrogen—hydrogen distances. The routine ADJUST is described by Memon et al. [6]. The kinetic energy is calculated in ADJUST. VALUES is called in a normal manner to calculate and store the thermodynamical quantities and the cycle is completed. 3.1.3. Program levels 3 and 4 At levels 3 and 4 are all the pure utility routines, as well a the user-supplied routines. We now return to the routine INPUT, The processing of the input data (called from SYSTEM) is shown in fig. 11. The first routine called from INPUT is TYPE, which prints out the headers. INPUT controls three separate NAMELIST inputs: SETUP to specify the main features and the thermodynamic state of the run. If an MC run is planned, the input parameters are given through the NAMELIST CARLOS. NAMELIST MOLDYN is read only if working on the MD mode. A detailed description of SETUP, CARLOS and MOLDYN is given in section 3.1.6. The second routine to be called from INPUT is COORDS. COORDS is a pure user-supplied routine to provide the site coordinates and fractional charges according to the potential model. The next routine, FIXCOM, is used to fix the centre-of-mass of t1e molecule at the origin in the molecular

287

I N P U ( LEVEL 3

C 0 0R B S

F I X C 0 11

U N I IS

I) P 0 L E S _________

________

R E T U R N

Fig. 11. Program structure in routine INPUT.

coordinates. FIXCOM also calculates the moment of inertia tensor and diagonalizes it. UNITS calculates the internal reduced units for the energy, length, time, mass and moment of inertia, and also the corresponding conversion factors to more convenient units. UNITS calls MPOLES, to calculate the dipole and the quadrupole moments. Several quantities are calculated in the routines called from INPUT, which are not actually needed in this standard version of McMOLDYN/H20, but which may be useful when the program is modified for other potential models. The construction of the look-up tables is divided to several routines. TABLES(TABLE) (called from SYSTEM) first calls ERFTAB(TABLE), where the complementary error function terms for the energy and the forces are calculated. The next routine PARPOT is a detailed description of the potential model. It supplies the potential parameters and charge terms for eq. (1) and an index list of the different types of site—site interactions EVFTAB(TABLE) (ENETAB(TABLE) in the case of Monte Carlo) is called with

288

A. Laaksonen

/

Computer simulation package for liquids and solids. I

.the information from PARPOT. In EVFTAB, the non-Coulombic energies, virials and forces are calculated. EVFTAB also calls LRCEVF where the integration from L/2 to ~ is carried out for the non-Coulombic energies and virial for the calculation of the long-range corrections in CTERMS (called from SYSTEM). Only the energies are tabulated in ENETAB. COORDS, TABLES, PARPOT, EVFTAB, ENETAB, LRCEVF and CTERMS are potential model dependent routines, whereas ERFTAB is a utility routine. Other pure utility routines are, for example ROTATE which constructs the rotation matrix and its transpose from the quaternions. ROTATE is called from UPDATE, ABOUT and CORREC. JACOBI performs a Jacobian diagonalization of symmetry diagonal matrix and RAND is the random number generator. MINV3 inverts a 3 X 3

PARAMS CONSTS TFACTS TABLES

look-up table parameters work array dimensions and interaction index table parameters to ERFC polynomial natural constants conversion factors to potential parameters potential parameters and fractional charges

CONVRT PHYSIC LRCCEV AVERIZ IOUNIT CORRCC MOUNTE

conversion factors and internal units physical specification of the system long-range corrections to energy and vinal accumulators of calculated quantities I/O units correction constants to corrector scheme MC control parameters

matrix.

CHARLO

3.1.4. COMMON areas

SCRTZR SCRTZL RDFRDF LRCVAL

MC system specification scratch area for real variables scratch area for logical variables RDF parameters and accumulators long range correction quantities

A rather large number of COMMON blocks

~5

used in the program. A list of these is given in table 3, with a short description of the content. Because the centre-of-mass and the site quantities are found in the work arrays there will be no need to change any dimensions in the COMMON blocks. Only the parameters NCOM, NSITE and NTABLE in COMMON/INDEXQ/ should be updated to correspond to the current dimensions of the work arrays. This is done in the MAIN program and BLOCK DATA. Several quantities stored in the COMMON areas are initiated in BLOCK DATA. Different accumulators are zeroed, values are assigned to physical constants and the default values for the input parameters are provided, 3.1.5. Work arrays The total length of the work arrays are: i) COM N * NOP ii) SITE N’ * (NOP * NSITES) iii) TABLE N” * NENT where N, N’ and N” are integers, NOP is number of particles in the simulation cell, NSITES is number of sites on each molecule and NENT is the number of look-up table entries. The routine INDEX (called from SYSTEM) returns the multi-

Table 3 COMMON blocks used in the program WINDEX WJNDEX WKNDEX

CONTRL LCONTR TPARAM INDEXQ

RANDNO

a scheme to split COM work array a scheme to split SITE work array a scheme to split TABLE work array integer control parameters parameters logical control

quantities to control the random number generator

plicities N, N’ and N”. These are different for each of the methods. INDEX also calculates the addresses where the work arrays are split in different routines. The contents of the work arrays in the different methods are given in table 4. In the work array COM, where all the centre-of-mass quantities are stored X, Y and Z are the Cartesian components of the positions of the mass points. QA, QB, QC and QD are the quaternions, and AX, AY and AZ are the x, y and z components of the angular velocities. 0, 1, 2 7 denote the zeroth, first and so on time derivatives of these quantities. FX, FY, FZ, TX, TY and TZ are the components of the centre-of-mass forces and torques, respectively. In the work array SITE for the site quantities SX, SY and SZ are the x, y and z components of the site vectors i~,as defined before eq. (22). FSX, FSY and FSZ are the components of the forces acting on the sites. The sites in the potential mode are illustrated in fig. 1. The array charge contains the fractional charges on the sites. The remaining part of the SITE array is used

289 Table

5

NAMELIST input parameters Table 4

Parameter

DESCRIPTION of the work arrays

value

X0,Y0,Z0,QAO,QBO,QCO,QDO SX, SY, SZ, CHARGE,(DUMMY)

TABLE: Leap-frog

TABE, (TABPE)3

NAMELIST/SETUP/ MC TRUE. for a MC simulation

Monte Carlo COM: ... SITE: ... ...

Default

Explanation

30

COM:

X0,Y0,Z0,QAO,QBO,QCO,QDO,VXO,VYO,VZO, X0,AYO,AZO,FX,FY,FZ,TX,TY,TZ

SITE:

SX,SY,SZ,FSX,FSY,FSZ,CHARGE,

simulation

.FALSE.

MD NUMBER VOLM RHO TRTEMP

.TRUE. for a MD

ROTEMP LMCY

rotational temperature (K) 300 TRUE. if the MCY potential is used .FALSE.

LCAC LTB NENT

.TRUE. if the CaC potential is used TRUE. if look-up tables are used number of look-up table entries

number of molecules 3)

molar volume (10-6 m3/mol) density (g/cmtemperature (K) translational

FALSE. 111 1 18.0 300

DUMMY) TABLE:

42 TABE,TABF, (TABPE)3,(TABVIR)3,(TABFOR)3

Predictor— corrector COM: X0,Y0,Z0,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3, {X4,Y4,Z4,X5,Y5,Z5,X6,Y6,Z6,X7,Y7,Z7) QAO,QBO,QCO,QDO,QA1,QB1,QC1,QDI, QA2,QB2,QC2,QD2,QA3,QB3,QC3,QD3,

(QA4,QB4,QC4,QD4,QA5,QB5,QC5,QD5, QA6,QB6,QC6,QD6,QA7,QB7,QC7,QD7) AXO,AYO,AZO,
FX,FY,FZ,TX,TY,TZ SX,SY,SZ,FSX,FSY,FSZ,CHARGE, (DUMMY)42

TABE,TABF, (TABPE)3,

(TABVIR)3(TABFOR)3 Constrained dynamics COM: (X0)3, (Y0)3, (Z0)3, QAO,QBO,QCO,QDO, (X1)3, (Y1 )~‘(Z1)3 (FX)3, (FY)3, (FZ)3 SITE: SX,SY,SZ,FSX,FSY,FSZ,CHARGE, (DUMMY)42 TABLE: TABE,TABF, (TABPE)3, (TABVIR)3 (TABFOR)3 where: X0,Y0,Z0,X1,Y1,Z1,X2,... Cartesian COM positions and their time derivatives QAO,QBO,QCO,QDO,QA1,... quaternions and their time derivatives AXO,AYO,AZO,AX1,AYI,... angular velocities and their time derivatives FX,FY,FZ,TX,TY,TZ,... COM forces and torques SX,SY,SZ,FSX,FSY,FSZ,CHARGE site positions, forces and charges

(DUMMY) dummy space to calculate the Ewald reciprocal space terms TABE, TABF ERFC energies and forces (TABPE),(TABVIR),(TABFOR) nonCoulombic energies virials and forces no. of sub-arrays of basic size. optionally chosen. <~~> not used in the second order quaternion method. For a more detailed description, see the text.

(.~

(...

),, =

}

ALPHA Ewald convergence parameter KMAX K max for Ewald Fourier sum KSQMAX K max squared LSAVE .TRUE. to dump the restart file LRES .TRUE. to pick up the restart file LGR .TRUE. to calculate g(r) ITIN unit for formatted input ITOUT unit for line-printer output ITFROM unit to pick up the restart input ITON unit to dump the restart file

FALSE. .FALSE. 2000 5.6 5 27 .FALSE. .FALSE. .FALSE. 5 6 7 8

NAMELIST/MOLDYN/ LLF .TRUE. leap-frog integrator FALSE. LPC TRUE. predictor—corrector scheme .FALSE. LQ2PC .TRUE. second order quaternion method .FALSE. LCNSTR TRUE. constrained dynamics .FALSE. NSTEPS number time steps 1 DT time step (s) 1. E—15 NSBPO number of steps between print out 1 LTS .TRUE. for temperature scaling .FALSE. NSBTS number of steps between temp. scaling 1 NSSTS total number of steps for temp. scaling 0 NTRAPC order of translational predictorcorrector equations 5 NROTPC order of rotational predictor— corrector equations 6 NITERQ number of micro-iterations in the second order quatermon scheme 1 EPSILO * “temperature” parameter to scale the reduced time step in the predictor— corrector scheme 300. LSCFT * .TRUE. for scaling the forces and torques (“catastrophe parameter”) FALSE. NAMELIST/CARLOS/ NGIANT number of giant moves 1 NMACRO number of macro moves 5 NMICRO number of micro moves 200 TRANP translation parameter (A) 1 ROTP rotation parameter (radians) 1 LMCEQ TRUE. for pre-equilibration FALSE. LMCEW TRUE. for MC with Ewald sums .FALSE. ____________________________________________________________ * Should be used only in the equilibration phase.

290

A. Laaksonen

/

Computer simulation package for liquids and solids. I

to calculate the Ewald sums. The TABLE array contains the look-up tables for complementary error function energies, virials and forces, respectively, followed by non-Coulombic energies, virials and forces. As mentioned before, these quantities can be calculated explicitly. In that case the TABLE array is not used. Although the quaternions are not used in the constrained dynamics, these are needed to initiate a new run and are used later as dummy arrays. 3.1.6. NAMELIST input

The input parameters are given through NAMELIST decks. A) The specification of the molecular system

The first NAMELIST input-deck is SETUP. SETUP contains all the control parameters cornmon to MC or MD. The input parameters are listed in table 5, together with the default values stored in the program. All the potential-model specific parameters are stored in the user-supplied routines COORDS and PARPOT and are not given as input. The box size is calculated either from the density RHO or the molar volume VOLM. If both parameters are specified it is calculated from the density. Default values for the Ewald sum parameters ALPHA, KMAX and KSQMAX are taken from CPC program MDIONS. The use of these parameters is explained in ref. [25]. B) Molecular dynamics Giving input parameters for a MD run, the type of the integrator has to be specified. This is done in NAMELIST MOLDYN. The magnitude of the time step DT has to be given, as well as, various parameters to control the I/O and the temperature scaling of the velocities, C) Monte Carlo The input for an MC run comes through CARLOS. The default values for NGIANT, NMACRO and NMICRO gives 1000 MC moves, Computationally it corresponds to 3—4 MD time steps depending on the method, when default values in table 5 are used. As a rule of thumb in specifying the displacement parameters TRANP

and ROTP should approximately give as many accepted as rejected moves. There are also two logical control parameters LMCEQ and LMCEW in the NAMELIST CARLOS. If LMCEQ is .TRUE. a special “dealyed” MC scheme is started to perform a faster equilibration when the sirnulation is started from a lattice. In this scheme the position and energy vectors are not updated after each move but after several moves. This procedure is meant to be used only in an equilibration phase. If LMCEW is .TRUE. the Ewald sum technique is used for the Coulombic interactions. The reciprocal space term is updated at every macro move.

4. Using the program 4.1. Use of the standard program The program is written in FORTRAN 77. It should be easily implemented on systems with FORTRAN 77 compilers. The program does not call any system routines, but if desired the function RAND can be replaced with a corresponding system random number generator. Some routines appear in the program with the same names as the corresponding system routines in Cray. For example SSUM, SDOT and CVMGP. This is done in order to have a “pre-vectorized” version of the program. Originally the program was written for an IBM computer in double precision. When starting a new run, the positions of the centre-of-mass points are taken from a simple cubic lattice. There are no special restrictions to the number of partides, because the lattice generator minimizes the number of vacancies in the crystal. The orientations of the molecules are randomized, as well as the velocities if it is an MD run. Temperature scaling should be used in the beginning of a MD run, because a large amount of heat is produced if the system is far from the equilibrium state, thus causing heavy fluctuations in the calculated quantities. A suggested procedure to start an MD run is to start with the MC mode and to equilibrate (or nearly equilibrate) the configuration. The Monte Carlo restart file is then picked up as a starting configuration for an MD simulation with randomized velocities. The restart file is an unfor-

A. Laaksonen

/

Computer simulation package for liquids and solids. I

matted file and contains a label from the run by which it is generated (i.e. MC, leap-frog-MD, predictor—corrector-MD or constrained dynamics). After the label follow temperatures (MC,MD), the time step (MD), the number of particles and some information concerning the content and the size of the restart file. If the number of particles does not match, the execution is aborted. The next quantity dumped on the restart file is the COM work array. At the end are the accumulators of calculated quantities. There are four I/O units in the program. The three NAMELIST decks are read from unit 5. The line printer output is produced on unit 6, the restart file is read from unit 7 and dumped on unit 8. The program does not dump any trajectories simply because there are so many different quantities available to calculate correlation functions or mean square properties. The user is expected to modify the program to collect trajectories of interest. All the line printer output begins with three characters: * * *, ??? or ! !!. * * * is a normal output from the calculation. ??? contains a warning: the user should check if it may cause any errors. The !!! case, finally, is always followed by an abortion of the execution because an error is traced. There is a large number of checks for possible errors and inconsistencies in both the control parameters and the calculated quantities. 4.2. Modification of the program This program was originally written for water simulations, using two different ab initio water pair-potentials. For water there is a large number of potential energy functions available in literature, both theoretical and empirical. If the user wishes to incorporate some other site—site potentials with fractional charges in the program (ST2 for example [25]), the routines COORDS, PARPOT, EVFTAB, ENETAB, LRCEVF and CTERMS must be modified: possibly also routines REORDE and ATOMF for the constrained dynamics, Also, if the program is desired to be modified for simulations of other systems than water, ammonia or hydrogen fluoride as examples, the routines where the modifications are needed are the same as for the new potential. Finally, it is simple to remove one or several computational

291

schemes from the program because of its modular structure.

5. Test runs Three test runs are included in this communication: A) MD with Verlet—Fincham leap-frog scheme, B) MD wUth Gear predictor—corrector scheme and C) MD with second order quaternion method by S nnenshen .

-

°

These are probably the most interesting cases of all the available simulation techniques in McMOLDYN/H2O. The three test runs were performed on a CRAY 1 computer. All the runs are started from a lattice configuration. Because of the limited text space, output listings from the runs with different equilibration procedures (i.e. mixing MC and MD or different MD methods as discussed in the text) are not included here. The input data to these three runs are given in table 6. The number of water molecules in the test jobs is 216 and the density is equal to 1.0 g/cm3. The simulation temperature is constrained to 300 K for both translation and rotation, In run A) the time step is specified to 1 fs while it is 0.5 fs in runs B) and C). The MCY potential is used in runs A) and B) while the CaC (Carravetta Clementi) potential is employed in run C). Radial distribution functions are calculated in run C). The CPU times on a CRAY 1 were A) 1970s (400 steps) B) 2960s (600 steps) and C) 2970s (600 steps + RDF). This standard version of the code is not very heavily vectorized nd the four-site potential model is rather complex to use in the MD simulations. The energy conservation is very good in all of the methods used. Concerning the water models the MCY potential gives a dipole moment of 2.19 D while it is 2.12 D for the CaC potential (the experimental dipole moment is 1.85 D). The quadrupole moments are (in Cm2) qxx 0.0, qyy = 1.32E—3 and qzz 5.28E—40 for the MCY model and qxx 0.0, qyy 1.21E—39 and qzz 5.OE—40 for the CaC model. The components of the diagonalized moment of inertia tensor are lxx 2.94, Iyy = 1.02 and Izz 1.92 (in the units of 1E—47 kg m2). =

=

=

=

=

=

=

292

A. Laaksonen

/

Computer simulation package for liquids and solids. I

Table 6 Input parameters used in the three test jobs A) Verlet — Fincham leap-frog ‘$SETUP MD = TRUE., LSAVE = TRUE., LRES = FALSE., LMCY = .TRUE., LCAC = FALSE., LTAB = .FALSE., LGR = .FALSE., NUMBER = 216, TRETEMP = 300., ROTEMP = 300.,

B) Gearpredictor — corrector scheme $SETUP MD = .TRUE., LSAVE = TRUE, LRES = .FALSE., LMCY = TRUE., LCAC = FALSE., LTAB = .FALSE., LGR = FALSE., NUMBER = 216, TRETEMP = 300., ROTEMP = 300.,

C) Second order quaternion scheme $SETUP MD = TRUE., LSAVE = TRUE. LRES = FALSE.. LMCY = .FALSE., LCAC = .TRUE., LTAB = FALSE., LGR = TRUE., NUMBER = 216, TRTEMP = 300., ROTEMP = 300..

$END $MOLDYN NSTEPS =400, NSBPO = 50, LTS = .TRUE., NSSTS = 300, LLF = .TRUE., LPC = .FALSE., DT = 1.OE-15, $END DT = 0.5E-15,

$END $MOLDYN NSTEPS = 600, NSBPO = 50, LTS = .TRUE., NSSTS = 450, LLF = .FALSE., LPC = TRUE., NTRAPC =3, NROTPC =3, $END $END

$END $MOLDYN NSTEPS = 600, NSBPO = 50, LTS = TRUE., NSSTS = 450, LLF = FALSE., LQ2PC = .TRUE., NTRAPC = 6. NROTPC =6, DT = 0.5E-15,

Acknowledgements

This program was entirely written during my stay at the IBM research laboratories in Poughkeepsie and Kinston (N.Y.). 1983—1985. The author wants to thank IBM Fellow, Dr. E. Clementi for his support and encouragement and also for the hospitality during this period. I would like to thank all the colleagues for sharing their ideas and giving constructive criti-

cism. Especially I would like to mention Drs. P.

G. Corongiu, J. Detrich, K. Hermansson, 5, Kandadai, G. Lie, H.L. Nguyen, D. Rappaport, R. Sonnenschein and M. Wojcik. I also would like to thank Ms. M. Van Duyne for introducing me into the fascinating world of the IBM Script Mathematical Formula Formatter, Ms. P.E. Bowen-Jenkins for correcting the language and Mr.O. Louhivuori for producing the figures.

Bopp,

References [1] E. Clementi, Computational Aspects for Large Chemical Systems in Lecture Notes in Chemistry, vol. 19 (SpringerVerlag, Berlin, Heidelberg, New York, 1980). [2] L. Verlet, Phys. Rev. 159 (1967) 98. [3] D. Fincham, Information Quarterly CCP5 (2) (September 1981) 6. [4] J.P. Ryckaert, G. Ciccotti and H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [5] W.F. van Gunsteren and H.J.C. Berendsen, Mol. Phys. 34 (1977) 1311. [6] M.K. Memon, R.W. Hackney and S.K. Mitra, J. Comput. Phys. 43 (1981) 345. [7] A. Nordsieck, Math. Comput. 16 (1962) 22. [8] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hall, Princeton, NJ, 1971). [9] See for example Integration Algorithms in MD by D. Fincham and D.M. Heyes, Information Quarterly CCP5 (6) (September 1982) 4. [10] R. Sonnensehein, J. Comput. Phys. 59 (1985) 347. [11] E. Clementi and G. Corongiu, Intern J. Quantum Chem. 10 (1983) 31.

A. Laaksonen

/

Computer simulation package for liquids and solids. I

[12] J. Detrich, G. Corongiu and E. Clementi, Intern J. Quantum Chem. Symp. (1984). [13] 0. Matsuoka, E. Clementi and M. Yoshimine, J. Chem. Phys. 64 (1976) 1351. [14] V. Carravetta and E. Clementi, J. Chem. Phys. 81(1984) 2646. [15] R. Colic and 0. Salvetti, Theoret. Chim. Acta (Berlin) 37 (1975) 329. [16] LV. Woodcock and K. Singer, Trans. Faraday Soc. 67 (1971) 12. [17] PP. Ewald, Ann. Phys. 21(1921)1087. [18] W.F. van Gunsteren, H.J.C. Berendsen and J.A.C. Rullmann, Faraday Discuss. Chem. Soc. 66 (1978) 58. R.O. Watts, Mol. Phys. 28 (1974) 1069. [19] W. Smith, Information Quarterly CCP5 (4) (March 1982) 13.

293

[20] D.M. Heyes, Information Quarterly CCP5 (8) (March 1983) 29. A.J.C. Ladd, Mol. Phys. 33 (1976) 1039. [21] D.J. Evans and S. Murad, Mol. Phys. 34 (1977) 327. [22] H. Goldstein, Classical Mechanics, 2nd ed. (Addison-WesIcy, New York, Reading, 1980). [231N. Metropolis, A.W. Metropolis, M.N. Rosenbluth, A.H. Teller and E. Teller, J. Chem. Phys. 21(1953)1087. [24] R.W. Impey, M.L. Klein and I.R. McDonald, J. Chem. Phys. 74 (1981) 647. [25] CPC Library Program MDIONS (catalogue no. AARM); N. Anastasiou and D. Finchan, Comput. Phys. Commun. 25 (1981) 159, erratum 28 (1983) 323. [261 F.H. Stillinger and A. Rahman, J. Chem. Phys. 60 (1974) 1545.

294

A. Laaksonen

/

Computer simulation packagefor liquids and solids. I

TEST RUN OUTPUT Test run I Ph~.slCAL / THERMOD~NA4clCAL SPEC1FICPTBQN CF THE 5151FF: CONSiSTS OF ~

216

PARTiCLES

/ CELL

AND THERE ARE

—

1.000000 KG/~DM**3 CORRESPONDING IC MOLAR VOLUME OF O.1TO2F—G~ ~*E3/MoL

THE DENSITY FOR THE SYSTEM IS GiVEN AS WHICH GIVES A FOX (CUBIC) LENGTh CF

1.8626 NM

TOTAL MASS/MOLECULE:

2,95151E—26

DESiRED TRANSLATIONAL TEMPERATURE:

(K)

300,00000

:

(K)

300,00000

DESIRED ROTATIONAL TEMPERATURE

MOMENT OF INERTIA: ~

DIPOLE MOMENT

:

QUADRUPOLE MOMENT:

~‘*

SiTE

(KG*NM**2)

2.93828E—29

1.O2ONOE—29

1.91788E—29

(C’M)

O,00000E+OO

O,00000E+0O

—7.314148E—3O

(DEFIE)

O.00000E+OO

0,00000E+OO

-2. 192C3E+G0

(C~44M**2)

O,00000E+OO

0,00000E+OO

0.00000E+00

0.00000E+OO

1 .31717E—21

O.00000E+OO

0.00000E+OO

0.00000E+O0

5. 28447E—22

MASS (KG)

0 H H

CHARGE (E)

2.65679E—26 1.673614E—27 1.67364E—27 0.00000E+0O

I,,

X (NM)

O.00000E+OO C.71750E+00 0.71750E+00 —O.14350E+O1

Y (NM)

0.000000 0,000000 0.000000 0.000000

~

CONVERGENCE PARAMETER ALPHA FOR EWALD SUM

‘~

NR OF STEPS PLANNED FOR THIS RUN IS

400

“

TEMPERATURE SCALING WILL BE USED UNTIL

300TH

~

N P01N~ MASSES / CHARGES (SIlES) I1~EACH PARIICLE

a (NM)

0.000000 0.075095 —0.075695 0.000000

0.006556 —0.052032 —0.052032 —0.020214

5.60000 AND THE TIME STEP TO RE USED IS

1,0000 FS

TIME STEP

MOLECULAR REORIENTATIONS ARE DESCRIBED USING QUATERNION ALGEBRA

BOTH TRANSLATIONAL AND ROTATIONAL EQUATIONS OF ROTATIONAL SCHEME BY FINCHAM

“

~

STEP ‘~‘

U

“*

B—TB

*‘~

E—ROT

~

E—TOi

MOTION ARE SOLVED USING LEAPFROG ALGORITHMS

~e*

TB—I

*‘~

ROT-T

TEMP (K) ‘~‘ VIR (KJ/MOL)

P (MPA)

50

—28.4266

3,7757

3.8202

—20.8307

302.75

306,31

304.53

—114.1208

2251.2506

100

—32.8847

3.7672

3.7347

—25.3829

302.06

299,146

300.76

—165.1373

31914.8722

150

—33.0661

3,7474

3.71409

—25.5778

300.148

299.96

300.22

_155.96146

30214.14216

200

~314.4693

3,75143

3.7216

—26,99314

301.014

296,141

299.72

_11l6.9070

2857.0898

250

—314.2203

3.7392

3.7315

—26.7497

299.82

299.20

299.51

—142.5886

2776.6255

300

—314.6252

3.7377

3.7593

—27.1282

299.70

301.143

300.57

—137.6008

26814.14323

~143.52114

2797.6509

—150.8401

2932.6540

~

TEMP.

SCALING SIsITHED OFF AT TIME STEP

300

350

—344,7188

3.83914

3.7561

—27.1233

307.86

301.18

304.52

400

—34.7765

3.8298

3.82CC

—27.1179

307.09

307.01

307.05

—

/

A. Laaksonen AVERAGE VALUES AFTER

***

Computer simulation package for liquids and solids. I

295

100 STEPS:

P01—ENERGY

(KJ/MOL)

—34.6~O4 +/—

0.1894

ThAN KIN—E

(K.J/MOL)

3.8004

+1—

0.1138

ROT

KIN—E

(KJ/MOL)

3.7146~( +1—

TOT—ENERGY

(KJ/NOL)

—27.1232

0,1~464

÷/—

0.0021

THAN TEMPERATURE (K)

304.73 .1—

9.13

ROT

TEMPERATURE (K)

300,143 ~/—

11,114

TEMPERATURE (K)

302.53

+1—

7.63

V1R1AL

(KJ/MOL)

—1143.2324

*1—

3.3530

PRESSURE

(MPA)

2190.1052

+1—

‘12.0707

Test run 2

*‘*

NB OF STEPS PLANNED FOR THIS RUN IS

**‘

REDUCED TIME STEP

~‘

TEMPERATURE SCALING WILL BE USED UNTIL

600

AND THE TIME STEP TO FE USED IS

0.5000 FS

0.1998E—O3 450TH

TIME STEP

***

MOLECULAR REOBIENTATIONS ARE DESCRIBED USING QUATERNION ALGEBRA

**~

TRANSLATIONAL EQUATIONS OF MOTiON 41LL FE

SOLVED USING

3Th ORDER PF.EDICTOR—CGRRECTOR ALGORITHM

*‘*

ROTATiONAL

SOLVED USING

3111 ORDER PREDICTCR—CQRBECTOR ALGORITHM

STEP *~

U

EQUATIONS CF MOTION WILL BE ~

F—IN

E—RO1

*‘~

E—TcJI

~

TR—T

***

~

ROT—I

1BHP (K)

‘*~

VIF (KJ/MOL)

P (FPA)

50

—20.1171

3.7779

4.0330

—12.3061

302.93

323.38

313.16

—95.4291

1905.4875

100

—29.1883

3.7500

3.7944

—21.6439

300.65

304.25

302.117

—116.9442

2302.5412

150

—30.7l425

3.7581

3.7533

—23.2310

301.314

300,95

301.15

—142.7410

2780.11483

200

—31.5253

3.7611

3.7389

—24.0259

301.58

299.79

300.69

—157.6247

3055.6442

QO*Q0

O.1000E+O1 +1—

0.1192E—O6

Q0*Q1

—0.6827E—09

+1—

0.3EC1E—0E

“‘i

250

—32.8985

3.7555

3.7290

—25.4115

301.37

299.00

300.18

~1l46.9533

2858.0982

300

-33.3193

3.7361

3,7380

-25.8431

299.74

299.73

299.73

-139.2992

2715.7255

350

—33,2632

3.7408

3.7257

—25.8167

299.95

298.74

299.35

—135.8174

2651.140141

400

—33.3407

3.7455

3.7355

—25.8557

300.33

299.85

300.09

~1145.0319

2622.0673

—26.1544

300.85

301.67

301.26

—154.2271

2992.441457

—153,1947

2961.6984

..~.

—33.6687

1450 *~*

3.7521

3.7622

TEMP. SCALING SF1THED OFF AT TIME STEP

C0’QO

0.1COOE+O1 +1—

0.1460E—06

4450

QO~C1

_0.4l4148E_09 +1—

O.5171E—08 ~

500

_33,93117

3.9778

3.7916

—26~1653

318,96

3011.02

311.149

550

—34.1902

4.1885

3,8380

—26.1637

335.85

307.75

321.80

~1140.0790

2746.8187

600

—33.9732

4.1096

3.6992

—26.16144

329.52

296,61

313.07

—1143.4645

2806.5420

296

A. Laaksonen AVERAGE VALUES AFTER

~‘

/

Computer simulation package for liquids and solids. I

150 STEPS:

—34.0518 +1—

0.2139

4.0575

+1—

0.1393

(KJ/MOL)

3.8296

+/—

0.1056

(KJ/MOL)

—26.1647

+1—

0.0015

IRAN TEMPERATURE (K)

325,34

+1—

11.17

ROT

TEMPERATURE (K)

307,07

+1—

8.47

TEMPERATURE (K)

316.21

+1—

8.56

POT—ENERGY

(Ki/MOL)

IRAN KIN—B

(KJ/MOL)

1101

KIN—E

TOT—ENERGY

VIRIAL

(KJ/MOL)

—1447.3097

+/—

6.5148

PRESSURE

(MPA)

2875.7575

+/—

115.7812

Test run 3 ~

NB OF STEPS PLANNED FOR THIS RUN IS

“~

REDUCED TIME STEP

~

TEMPERATURE SCALING WILL BE USED UNTIL

600

AND THE TIME STEP IC BE USED IS

0.5000 ES

O.1998E—03 I5OTH

TIME STEP

~

MOLECULAR REORIENTAT1ONS ARE DESCRIBED USING QIJATERNION ALGEBRA

~“

TRANSLATIONAL EQUATIONS OF MOTION WILL FE

SOLVED USING

5TH ORDER PREDICTOR—CORRECTOR ALGORITHM

‘~‘

ROTATIGNAL

SOLVED USING

5TH ORDER PREDiCTOR—CORRECTOR ALGORITHM

‘*~

SECOND ORDER QUATERNION SCHEME BY SONNENSCHEIN WILL FE USED

EQUATIONS OF MOTION WILL BE

1* 5*****~****~****5*5J~****~*4

STEP *~‘

U

E—TR

~

*5 *****I**********

B—ROT

‘~

*~

E—TOT

******~***5******5

TR—T

~

5*5 *I**5*****I*************K*****5**5I******5

ROTCT

~

YEMP (K) ~“

VIR (KJ/MOL)

P (HER)

50

—24.8279

3.7714

4.0023

~17.05112

302.40

320.92

311.66

—45.3215

976.1261

100

—32.8013

3.7500

3,7963

—25.2550

300.69

3014,40

302.55

—68.6282

1408.5685

150

—3Ll.2029

3,75711

3.7327

—26.7129

301.28

299.30

300.29

—95.5919

1907.7390

200

—35.6151

3.7588

3.7306

—28.1257

301.140

299,13

300.26

—111.0108

2193.0827

5QO

O.1000E+01

.#/—

0,14&OE—O6

QO~C1

—O.1219E—Ob

+/—

0.2100E—O7

QO 250

—36.5707

3.7582

3.7577

—25.0548

301.35

301.31

301.33

—103.2174

20148.8625

300

—37.1620

3.71400

3.7374

—29.6906

299.89

299.68

299.76

—97.3702

1940.0005

350

—37.2381

3.7472

3,7394

—29.7515

300,116

299.84

300.15

—97.5584

19443.7476

400

—37.2516

3.7385

3,74911

—29.7637

299.77

300.644

300.20

—103.5452

2054.1977

-37.50146

3.7441

3.7636

-29.9971

300.21

301.78

301.00

-109.88544

2171.7130

450 TEMP.

~

Q0500

SCALING SW1THED OFF AT TIME STEP 0.1000E+01

+/—

O,1333E—06

005G1

450 —O.8760E—09

+/—

O.1799E—07 ~~l~

500

—37.3322

3.8840

3,44498

—29.9984

311.43

276.62

294.03

—105.8073

2101.4375

550

—37.6761

3.8938

3.7829

—30.0013

312.22

303.33

307.78

—116.3635

2297.1178

600

—37.61445

~j.7239

3.9182

—30.0024

298.59

314.18

306.39

—102.7596

2150.1357

A. Laaksonen ~‘

AVERAGE VALUES AFTER

/

Computer simulation package for liquids and solids. I

150 STEPS:

POT—ENERGY

(Ki/MOL)

_37,48714

+/—

0.2270

IRAN KIN—F

(KJ/MOL)

3.8299

+/

0.0752

ROT

KIN—F

(KJ/MOL)

3.657e +1—

0.2167

TOT—ENERGY

(KJ/MOL)

—29.9997 +1—

0,0014 6.08

IRAN TEMPERATURE (K)

307.10

+1—

ROT

293.30

+/—

300.20

TEMPERATURE (K)

17.53

+1—

9.07

VIR1AL

(KJ/MOL)

—111.6905 +1—

14.2004

PRESSURE

(MPA)

2208.2894

TEMPERATURE (K)

+/—

78.0236

297

298

/

A. Laaksonen

Computer simulation package for liquids and solids. I

* * as ** Sa ass * *5 * s__sass_sm asss * * 5*555* 5* * msa**55**55555****5555sa55*ssssmsaa.s 5*5 5* •sssss*sss..a * 5~* 555*5***s*sSas5a s__as_s-s ,aasssas*sssa * *555 s*55*5*ss*s**sS* **ssS ass__s ma ssas**s*sa*sa 55*55 *assmassssasssss * ssmsa * as_as * * — * * * * ma * * 55*5 * 5585*aas * * *55* ss*s * 5**s*~s5 * S**** * 5555*55*5 * 5~~5 * *5* * *5* * * *5* * *55*5 * * * *5* a*5~5 * *55*5 * 5* * S•55* * sss~ ssaass * s__s_as_a. ass * * 555*55*555555 as__s_a s~mss * * ssasassss. ass*a * 5*5 * 5*55*58* 5*5* *5* *5 s.sss~.*s.ss*as~ss, a *s*s..ssssssssssasss a_s_s sssassssssassa.sa 55*5 s*s*sss.asss.sas * ssaassa*55 * * ss**ssssss sssasssssass.ss s__s ssssasss ssssssassssss 5 * sss sssSassssasssas. *5 s*s**saas *~s 5*55*5555 sea s_sm_s... ~~sasss sm__sm__s_s_ms sassas ss*sa*sss ~ss~ss ***~555 ~ sS.sss**as. 5a555S5s~s*5as~ 5~5 ~55 _~5 5555555 ~s~sss s*saas*ssssss*~a ~ sss*s.ssass.sasassa **55* 555*5*555 555s5*~*5 sa**ssa*sssssssssss ass.. a a**Sssas*saa.sss * m.~aa~ ~

I-’ C/~ 0 In

0 ‘—I

~ CLI ~ 0 -I 0 I

IZ -0 a)

0 .S-,mS 0

0 0 I_ o a)

,C.J

4a)’-OacSJCCJIcl.ION 0’OInU’SCflIOCflC,JCJJa)OnJflo.-NN-Lfl noo0oo’SornooLna)nnfla)nncCJ.~o.o~cC4nn4Cooa)Lot_C0N-0Cea’Ccna)Lnqs0.

—11410 10100.-ID IN

~~LflO’Sflo~-c,j

a N-04N-,~fla).-10a.rnIOO

4U’CN-cOO”IflLflN.IN.-OdL’CN-CCflSOa).-IOLCCOOC’0a),..fl 1400(0 0(0 4•#LOLOLOLOIOC4IOIO ON-na) IN LI’S (1 N 4’ N-1fl Nfl C) 0(0(5410 a) COCOa) 0,0051(1 ~000a)0OSO,00Os0O500O,O5O’S

0 0 0 0 0 0’—-~0’--’--0’-~’--’--

0 0

0 0 0

—

0 0— 0 0 0 0 0 0 0

1(54 (00)

.~ .~

c

(‘51051010*

010 ‘—

0

n

Cs.

00000000 4-100 000000

-~

1-100000000

IQ

N-N-a)a)COCCCCCC000.-.-’04CC410

‘0410 (‘La) 0410104 ‘0 InOCO t~-1— CO CNN- 0 as (-Lfl05fl CSJ,—O)’0 4 I01000S1OLOC) 100505*

Cs.

0

000ooooanjor—,--

(fl’ —

‘-

*5*5*5

-‘-4 ~

a) . Os

~-fl

Cs. 00 0 000 000000000

(Xl I4

DO

a)000000QQ

0 0 0

0

0

0

0 0

0 0

0

ID 000000000000000000000 —~ oso ‘—N (04’ LflS0N-CL) ‘050010

010

‘-~~

.041010

0000001D00000000000000 0101-N (‘04’ (5(0 N-a) 050’~~N (041(110 N.a) N-a)a)a)C)a)0) 0.0,0.0,0’

N-a)

~ 0 10

00000000 (04’ 1010

.4 ~

.1. ~ N ~90000000

~

~aoooooooo,

(~-C1)

-4

* —a *5. as *55 5*5

0 ‘-4 NJ

as..

s_as 5*5* sass s_as

s_ms 5

a_a~

s_s_s *5**s

0. (SC

— (1)

0

‘(-1

0

CO 0’

o

0 *sa**S5*m****5 0

a

0 0 .L..~saas55a55 0

sa

S sassams

sasaa~sa

5~5*5 s_ss_ a *5*5* * *e~s* a 5 * 5 *55*5 as 5 as a 5* *~* 5 .55*5 ** a as a*as ass * m.ss~saa 5 5*5 *sssa sass 5*5 .sa~ssas a ~ss 5*ssa*5-5***a*5 swss~~s~ *sss*5*s***s**sssaSs a * amass * s*asasss asssa*sssssss 5 a s55*55* * *5*55*5 ~ * .s**s * **smss a * ~a * *5* * *aa5~55*s~5*ss **aamas 5 * sss*ssssa*s ass 55 * 5,55 ~ s*ssas * * s.s*a*ss ma * ss~assaas~s*sss~aa sss*SSs*Sssss** m*ass*assa*sas a sass__sam * * * sasss a * asss.a* ~ vs a *assaSssssa * *ssssSs.**sa **sass *s*S*a.a S * ss~s.ass*sss~.ss*s a sssssssss*ss*ss* a ss.ssssasssass s**_ss*s*55s__sss sssss*sss***aSs* ***a***_*ss_ssa S ~ 555*5*55*55*555* * 5* **Sss*asaa 5 5 5 *5* 5 5 5 ~ a a * a ~ 5 5 a * a a S a S S ~ a as S’s *5 a seams. sassssasas*saaa*s*aas *

a

0 1(1 4.1 LI

a

ID

(S.

a) 0 1-’ ID CO CI) C-I

ID

5, ID

a 000

0 000 00 0 0 0 0000 000 0 00 0 0010 N 4 100) ‘~ “1 IN N- 0 CS N- N-CO N 10 0-4 0000000000000000000000000CONa)NOC1Sa)0N-a)a)fl,-0,N-1(4.na)41flfl100,n,01010a)C)N10 4-,000000000000000000000000000’CINCO0400040)10010NCOCOflN,-0a)N-CON-a)’0,.-CO10

z

(51 0 0 0 0 0 0 0 0 0

N 0 4 4 4’

II’S 10 10 N- N-

CO

CC’ 0’

0

N (0,0

‘~

05 0-0

04 “-1 .0 4’

1(5 10 N .-—,--.-‘-‘--‘-

a)

10 111

05 ._

0

I’d (0

(~.0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a In o 1’) 01(4 CCC 10 1(10 C) N 4’ 1010 111 1010 0 0 a) .0.0 0) 0 11511) Dl I’d , 01 .4’ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 N- .‘ COO 0) 4 ‘-~ N- N- 104 N-IC 0 01 . 0 .0 ‘0 151 0 4’ IlL LI’S 05 .21 5010 05 1-’ CO000000000000000000000000010N-01010a)0000N-.-N-000a).-0500CC000CC0’a)CCa)CC 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 151 N N 0 0 0 0 0 ‘~ ‘~ 0 0 0 ‘~ ‘ 0 0 0 0 0 0 0 0 0 00000000000000000000000000

0 (51(04 Z00000000 00

a

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

a)

IrSIC

N-CO

050

0

C~J (044010N-a)

0~’’~’-’-’-

0 0 0 0 0 0 0 0 0 0 0 0000

0 0000000000000000000000000000

050 I’d 104100(010 0’0’ N (041010 N-CO 050 N 041(110 (5JN(5JNNCNNNNOJ1flIflInIn1flIflInIn(flfl*4~fl4fl~fl**4~LA01fl10InL0I(5

0 0 0 0 0 000000

0 000000000

N-a)

050’—

0 000000000000000000

N

(nfl 1(5(0

A. Laaksonen

/

Computer simulation package for liquids and solids. 1

a 5 * * a 5* scm. *5* 55~5 5 5555* 55* *55*555*5*55*5 s****aS5*555~5*****S *s *a5555555555*aS sa*ssss*Sss*S5*Sas*55 *sssaaassa*sss*s5*ass**SaS*aSSsaSasS S ~ 5 * 5 * 5 * 5 5 5 * * * * * * 5 * 5 * 5 5 * 5 5 5 * 5 5 5 * 5 ~ 5 * 5- 5 * S 5 *s*assss*ss_*s*Sa*sa**sSs*ssa**s__a*s_ ssSsssSamssSSs*aS S saSa***S*ssS*S**a*** a*ss**ssaas**ssaa*s*s*55***sa*a****s* * 5 * * 555555*5 * * * a * * sssassSs * 55 * 55 * ama. * * * 55 * 5S * 555*5555 * 5 *5555555*55 * 55 * 55 * * * mm * aSs. * * 5* * *55555*sssss * as * * ssasSSSas ISs*as*sS*s***S 5*5 as*asasss*saaasass** 55555*5555*55*555 555555a55*5*aS**5*5S5 msssasssas*aasmssssaaaasss*ssaa*S***s* a*sa**.s,S*SaSs*ssas*.saaa*Sa*a*a**SSS * smaS * * * * 555*5 * * 5 * 5*55*5 * assssa*a * SSaS 555555~5S55555**S5S*5S5**5*5S5*5S*5*5S *aasssssS.aSs**Saa*Saa*.Saa*saSSaSs*a* *5*555 *555 *SS*SSSS a*5*a*s*a*sSassmssSS * *5*5 * * * a * * *5 * * a * * 55 * * * * * * 5*5*555 * 5 * 55

299

0

(/) 0.

‘(-C

14 CI)

L’S 0

a) Cs) > 0 14 I 14 0

0 .1-1*5* a

*

*55

*

ms

9 a) 0 C’S

~ooooooo~~~

jLI

(‘40000000000

‘- C00,00LOO’Ca)OLCO 04N-0’N-0544’-NN 05Ifl’5.-0s0,flC)fl10C~LON-NlJCN-N-4N-Iflfl

4 .0

(5,0

0 O

~

10N.~~N0CCCL(O1(lC0IflL0CII10’-CN~,flN50OJlflNL0N’-10’-NL0OIflO10LflInLflN0C”NLfl’ 00500) ‘-050)5401(10 IflGbOINQ,N-C0a)’ILSCON4’

4’

100ILOCOflOa)CO

0044 1flN—1J’’-I’dU’lSOC5~,flU’SCCCO ,((1010.-,05005N10a)14(0O4’N-.-10C)(54SO0( (OIflIOlflIflCV’5fl4 410101010101010 .ON-I~N-N-a)a)a) 0.0105000’-’N-COO 05.01010104’ N 1(SN .-0S(fl,fl.4’CQLCS.-OCOLj’ 1Sj0I0OOI’d~C0 .I.C

((4

~(5J1(j

0000000000 4.19000000000

.0

000000000

00000000000

a)0000900000 0 0 0

‘-

0 0’-’-0’-’-’-’-0’-’-

‘-0

0 0 0 0 0 0 ‘-0

0 0 0 0 0 0 0 0

00000000000000000000000000000000000000 ‘-N 104 U’(SO (00)050’151 0*1010 N. CC) Qs0’-N 104’ 1010 N-a) 050’10,050C0’0S0501010 N-N-N-N-N-N-N-N-(0N-10a)C)C)a)C)a)a)a)a)990’O500’0’O’O’ 00000000000000000

154104

Lfl10 N-a)

000000000000000000000

0

0 0 0 0 0 0 0 0 0

10 ~ .1)

.-.0000000000

~

aoooooooodd

X

‘-N

(04’

4010

N-CC) 050

Z009000000’-

‘C 0

a

*5

5*

*saa* *5*5*5 5*_ass am_ass * a * * S 5s

sa*5*5a * assasasas s*s**as_s **sasS_a_*am*s a*saass**ss***sas Ss*Ss * SS*S * 55555

a *5 * 5* Sa * Ssss*ss *5 **s*s*ss* * 55*5 * *5 555 * 55 * * * * 5*5*S***55555 * m*asas* * * *5 * * *5* *5*5*5 Ss**S~55*a* ssss***a ssS*** ssa 555*555*5* a__S_Sm *ss****Sss*s** * *5 as sass * a * as * * ss,*a ** * s__a_a * * a *5 55*5 as as * mm * * as *5555* * Sa a*S*ass * as * * 5 * 55*5 *55 s*a* aassSsasSS * **s*s55s5 * masS * * *5*5 5* * * * 5 * * 5* * *55 * s*s * *5 * 5*5* * 5*5 * ss a * * sea *55*5*555~5 *55_S 5*555 * Sa* * * * * as * * 5 * 55 * 5 5 *5*555 * a * mss*aS * * * * **5*~ * * * S * 5*5* * * 55 * sssa sass * 55 * 5 * * sass. * * ass * em * *msss * 5* * * * mesas * s 5*5*55 * *5 * * ~s*5S * *am *5* * * s a * *smsas a * * * aSsa a a as * * * 555 * 5*s * * *55 * * sass * * * * * ass * a a a 55*55 * 5 * *55*5* * a * *555*5 * * 5555*5*5*~5 * ama * ass * a * ,** * sass * 5~ * 5*555555*5 * a a * as * * * * sa * a * * * * * a * * * a * s * ass * s * 5*5 * 5 * ss * 555*55 a 5*55*55 a ass * * a * sssmssssas* * * a * as * sasS. *5*5*5 * * a * as * * * a * * 5 * 55555*5*5*5*555* s**msSas * 5**Sss**ssm*S* *

55555555*

00 000001JNN-O’-CC41001’dI’d 0 0 0 0 0 0 0 0 .‘ CC., 0 CL) 4 51) 05 N.’O 000OO00OOONC”C10sON-eJOsO~

00 00000C0(’-’-IlJ IN1p’-C’-N-J(0l 00 00000.-C0104’ N10”1fl,fl4 00000000010(00(01011(10,511 0

0

0

0

0

0

0

0 0

0

‘-0

00000000000000000 os0’-CM .04’ 11110 N-CO 050.-nj 00000000000000000

4N~0

a’-’-’-’-’-

0000000000000000

O

N’SO C.)

0 0 0

0

(0

10 0

0 50 0

‘ON

0

0 (04’ 1050 N ‘dilL 0

0

00

1,0 ((‘C 114 CCI SO ‘-

NC’dInfl10101UN-N-0J

N’ ((CISC 0JN-1flI5(L”..O’CIS.COOCJ,flL..4’ (0 CC 45CC) ~(04 104* (1400 ~

* ~S

0,Os0,050%COCO C

0

0

1

0

0 0

0 000000000000000000000000000000000 (‘-10050.-N (041050 (010010.-N (041(550 N 0

0

N 4’

In aJ (0 04

‘- (00)

N ‘-

CCI 0’ 0

LYCLLSQ’-,(4(’O(OflIOCO(’-1000ICJ,fl4lCN-CYC N N (54 NNJ

000000000001000000000000000000000

0)1010

0.0)

0

0

0

N-C)

0

3\05005000,00

0 0

O50’N

N N

0

-

0’

10.04050

-

0

N-CO 050

300

A. Laaksonen

/

Computer simulation package for liquids and solids. I

* S a ama as 5*5 * *5 5 5* 5* 5 *55555 5*5 SsmS*ssSa*a*a*sss*smS*s *S* * s*sa*a * am_as * ss * sas * saa* * 555 * * a * s***ss* * ssss*sssasas * as Sass S * 5 * *555*5555_a * asasaass**ms*s*aas * *** * sass *55*5 * * S 55* * sss*ss*S * *55 * *SS* * sam * *sssss * * sS*SS * asS. * ss**mm**5 * a * * Sa**SSS*55*5 * 555555 S * 555*55* * *5 * * 5 * *mmsssaa * ma * *sma * SaSa sass * * aSs 5 * * * 55*55 * *5*5*5 * 555*555*55* * * * * 55 55555*5*5 * * Ssa*Sa*5*s * Smsssaa**SSS a*Sssm*a*aa*saa 555 * * 5 * * 5 * 5*5 * S * *5 * * *s*ss * * ssmm * ass 55*S * * * *5 * 55 5*5 * 55 * as * *s * a_ass * sas*sassas * sasaS * 555* 555*5*5 * * * * 5SS*S * * 5 * *S55* * 5 * 5*5 * 55 s*sa * a * * * a * m * *5 * 5*5 * 555 * *5* * sSSS* * * *5* * * *55 * * smssss*s*msas * SS * a * s**S * s*SS as ass ** * ms_s_S_a * *msa**a * * * 55 5 * * * sa * SS*5 * m**asssssa * 5*5 * amasS *5 * * * * a * * * 5 5 * 5’ 5 * * s a S S 5 5 5 5 * * 5 * * 5 * * 5 * 5 _ass*_a*s * * * *5*5 * * S * * * * * ss * * a ** * * *

N0’40110CO00I0100100N-10O\05C)L450~N-C0N0,1010N0N0COC)05 ON N-C(COCCOCONO.CCO ION-II’S 10N04 0(0101010405100) ‘—010(010’10CO.-00.-.-1010C)NC0NN-10C)S04’NOONN4100510C0.-N-flGCIO.-

01(0*

N-OS’- 1010 N-0S.-fl1fla) 0(flIfl(00(010 OCN 1010.-fl (‘-04’ ((-.-.010 COCOflflflflfl40101010S0I0I010(0(0N.t’-a)a)W0,05050001-~’-’—NflJN

‘-4001

N-a)050a)LOOS10N.-.-a)*(OC-NC)NI’0IOOSOC01OCSC)flSOl00010 .-Ost..IØN.-.-0.-QCOOION-a)1000N-4a)N1010flN-N-CO(0N-(0NN(0 0010,0,0 00,000,0 0 0 01010,0005010,0 Os050105010105Os01050’C)

—

0

0

0

‘-

‘-0

‘-‘-0’’-’-

000’-’-

000’-

0510

000000000000

000000000000 0000000000000000000000 10I0(0a)Os0.-N(0fl10C0N.CO050~NCOO10I0N-a)01O~N10fl10CON-10 (0 10 50 CO CO N-N-N-N- N- N- N- N- N-N-a) a) 100)10100)0)0)0) 0,0.050,010501005

1

0

0

0

0

0000000

0000000000000000000000

55

*

* * * * * * 5 *

a *5*5 *5 555*5 55* 5*5*5 as_ssm 5* * * ama sm.mma*sas* as S*as 55 5 ass * * *ssss*a*ss***a * * * * * s*a*ssas * * * * 5* *55* * 55 * * S*S* * 55 * Sass_ms * * * *S55 55 *5 * * * 5 sms *5 * * a * *SS 5*mS*mSsas*s*5*S*aasSsSSaSm*ssa * ass 5555* * 55 * 5 * 5* * * as * msms*m55 ** — 55*5*5*555555 * SS*S*55*55555 * 5555 * * as * as *5 -, 5 555 * 5,5 * Sassams * * mass * s * 5s * as * * 55*5 * * * 5 *5 a s**s* * * 5*555 * 555 * * * sass s 5 * 5 * * * 5 * 55*5*5 * * .555*5*5 * 555*5 * *sssssas*a* sSas * mss* * * * * * * * * * S*am* * * 5 * sssass * a*s~s*sssass SS * msS*ss* a * a * * ssmsm * *55***5a**5*5S * *5*5*5 a *ss * S * * * SS 5 * as*5**S***5S 55 * * S*5**sss asSsaSs*ss*aa ass * as_a * sma *5* * Sa*s*ssss*5555 **m *sSmsSa*a * * 55 * mSs * SSas***aSS5 * sass * 5*55 * 5555*5 ***m * 55*555 s*s * * 5*5555*5 * * 55 * *5* *5* *a**mssaSS*aa * 55*55 ss** * a Ssssss***a * * ass * a * sss * sass a * sass * ss ass * sSSaS*SSS* * * 5 * * * *555*5555*5555*55 * *5 * * SSSSS assaSs aaSaS*SS* * a * ass* * 55555 * *5 * *SSaS aSs s*5s *5*55 * 555*55*55*5*5 * **mS55*a*SssS * * 55555 * 5 * as 555*5 a * a*s* * * * * a * 55 * 55 * * 5 * 55 * * * Sm 5 * SSSmm*s* * a 55555S * * SS*S * * *55 * Sa * ssmSs * **SS * Sass_s. * * sa * *55*

55*5

* * 5 *

55

0000000OflN-NNNfl0’0)N-CO1’dN010N10O10100.-C’dCO10COOs10N-S00CN.10N&.10~~I0CO0COCO1010.-N0000 00 0 0 0’- 1050 N-a) 4’ 0 (“COCOS 04’ (0100 NO (0 (“Ca) 0CC14N-CYS(Ofl N•~ 10 N N Na) N- 4~ N-Nfl LOON 00000000000’(05004’ a)’-fla) ‘-4’ 0)100)40 (‘-4’ ‘-01504’ NO a) N-N-N-N-N-CC) 05010100) NS0 00000000000000.-.-

~~NN

N (0CO(04’~0C01NN(04’10I0Na)Q5.-

(54104’ ‘010(11010

OS

NCOOCON-OSQNJ*U’S N NN N N N N 10 (“C (“C (51

0000000010OSIOflNNCOSOOS(nCOOsa)N-10fla)10.-010a)1OflC000’LOa)S0000flu’SO\CflN-COCOlON-OOsO 00000000NN-flOsflOSa)S0Osa)4’.-CON-N1O.-flOLflN-LflCO.-N-01fl10fl0O’LO1OCOO’-a),’01010N-Na)fl10a) 0000000000N4’ OSN 1(5(54 100a)C)C)N-0)0100’-’-’-000050101010,0010’IN0.01OsOs0’0S0’05Os01000(

0000000’00000

0’-’-’-’-—0~OOO0

0.-.-’-.---’-.---

00000000000000000000000000000000000000000000000000000 ‘-Nfl~C4~LOL0N-a) 0.0~N (nfl 1050 ((-10010’-N (5(4’ 4(110 N-C) 0.0.-N 01 N N N N N N N 01 N 1010101010101010(0(04’ 4’ 4’ 000000

C

00000000

0

(04’ COCCI—a) 050’-N 4’ *4’ 4’ 4’ 4’ 4’ 1010Lfl

00000000000000000000000000000000000000

0 IflO 1010 N-a) 010’N (04’ 101(1 1010101010(010(01010