LUSim-A Macintosh Program for Simulation of Molecular and Polymer Systems Olle Teleman
..
VIT ,
Biotechnical Laboratory POB 202 SF-02151 Espoo, Finland Bo ~onsson' Physical Chemistry 2, Chemical Center POB 124 5-221 00 Lund, Sweden Sven Engstrom Food Technology, Chemical Center POB 124 S-221 00 Lund, Sweden
If you ask the students in an undergraduate class about the size of a water molecule and how far i t translates and rotates in a second a t room temperature, you will probably get a very broad range of numbers, especially concerning the motional characteristics. Everything from 100 to billions of revolutions per second will be given as the answer. Although the molecular picture has been accepted for almost a hundred years, the thinking in molecular terms bas just begun. Simulations of Molecular Systems The progress made during the last decade is, to a considerable extent, due to the appearance offast computers that have made it possible to imitate the structure and dynamics of molecular svstems. takine into account the intermowork was made by lecular interactions. he Metropolis et al. ( I ) and Alder and Wainwright (2) during the 1950's. The simulations of liquid argon by Rahman (3) and of lionid water bv Rahman and Stillinzer - ( 4 ) also represent mijor breakt&oughs in the field. There are three methods available todav for simulations of molecular systems:
Monte Carlo (MC) Molecular Dynamics (MD) Brownian (or Stochastic)Dynamics (BD or SD) With the MC method one generates, in a physically sound manner (e.g., via the Metropolis algorithm), a large number of configurations of the system under study, which are used for evaluation of static properties. Monte Carlo programs for polymer simulations (5)and simulations of nematic liquid crystals (6) on a personal computer were described earlier in this Journal. In the MD simulation, Newton's equation of motion (i.e., F = m . a), is solved for a system of interacting atoms. If molecules are to be treated, their rotational motion must also be taken into account, for example, through Euler's equation of motion. In contrast to the MC method, time is explicitly considered here, and an MD simulation will give both static and dynamic properties of the system. (The MC and MD techniques are briefly described in Figure 1. See ref 7 for more details.) In the BD simulation, the effect of the solvent on a solute molecule is considered by means of a stochastic approach. A Statistical Mechanics Simulation Program: LUSim
In this article we introduce a statistical mechanics simulation program-LUSim, written for the Macintosh. I t can be used to simulate three different systems, a monatomic liquid water
an ion in water by means of the MD technique, and a fourth system, polymer melt with the MC method (see Fig. 2). LUSim calculates a number of relevant static and dynamic properties, such as potential energy pressure temperature as well as radial distribution functions for the evaluation of molecular translational diffusion time correlation functions for the evaluation of molecular rotational diffusion We show how LUSim can be used as a tool for teaching about liquid properties and modern simulation techniques using computers within the Macintosh family. Throughout the text we will give examples of suitable exercises to be carried out as numerical experiments with LUSim. Theory and Methods Interaction Potential
I n order to imitate a real molecular system by means of computer simulation, one must describe the intermolecular interactions. I n practice one assumes that the system's total energy can be expressed as a sum of interactions bee . intertween airs of molecules. However. in ~ r i n c i ~ lthe action between any two molecules wiil dep'end upon the positions and orientations of all molecules in the system, due to polarization effects. The polarization effect is too complicated to take into account and thus is either neglected or included only in an approximate way (See the water-water ~otentialbelow.) Fieure 3 shows the input menu for oneof LUSim's system; one ion in water, &at may serve as a n illustration for the other systems as well. LennarMones Particles
The simplest molecular system, from a n interaction point of view, is a monatomic liquid, for example, liquid argon. The interaction between two argon atoms can be accurately described by the so-called LennardJones potential uw.
where r is the distance between two atoms. The meaning of the two interatomic parameters & and o is illustrated in Fimre 4, which shows that & is the depth of the potential wen, and o is the size of the atom. The first term in eq 1represents the repulsion (electron shell overlap) between the two atoms, and the second term represents the attraction (the London dispersion interaction). The r dependence of the second LC& can be derived from quantum mechanics and is of the induced dipolcinduced d i ~ o lM e e . The distance de~endenceof the re~ulsiveterm isbf sh%er range than the'attractive term, but the exponent (traditionallv - 12) . has onlv been chosen for comoutational convenience. Figure 4 afso shows how a variation of the E and o oarameters affects the interaction wtential. This variatidn can serve as the basis for a s t u d e k simulation exercise. 'The executable code, a manual, and limited support are available for US$100 per copy. Please contact Bo Jonsson for further information. Phone: +46 46 10 82 39; Fax: +46 4610 45 43. Volume 70 Number 8 August 1993
641
In order to perfrom MC and MD simulations, an explicit energy expression is needed in terms of atomic positions: Nri, rz, ..., IN). The Monte Carlo Method move accepted In the MC method confiaurations of the ensemble are uenerated according to their ~olt;mann weights, exp(-u/kT). ~ h i impli ies that low energy configurations are generated more frequently than high energy configurations. Calculation of ensemble averages becomes simply a summation, p = Zpi, where the sum runs over all generated configurations.
-
configuration n 1
..........
calculate desired property: p-1 move a molecule a random distance configuration n calculate configuration energy after move: vo IF
v, 5 u ~ move r OK and un = vn ELSE pick a random number a E [0,1] and compare with exp(-(vo- upi)lkT = Au IF
a s Au move OK and un = vn ELSE move not OK, move back and un = u n-i ENDlF ENDlF calculate desired prperty: pn
move another molecule....... configuration n + 1 calculate configuration energy after move: vmi
I t is not sufficient to describe the interaction between two water molecules by means of Lennard-Jones potentials alone because water is a polar molecule with a large permanent dipole moment. I n LUSim t h e dipole m o m e n t i s acc o u n t e d for u s i n g t h e S P C (Simple Point Charge) model (81, i n which t h e oxygen i s given a charge of -0.82e (e i s the elementary charge) and the h y d r o g e n s , consequently, a charge of +0.41e each. Because the OhH bond in the SPC model is 1A and the H a H angle is 109.5', the dipole moment of an SPC water molecule becomes 2.27 D (Debyes). This number i s significantly larger than the gas phase value of 1.85 D. In liquid water, however, the molecules polarize each other, thus increasing their dipole moment. The SPC model takes this many-body effect into account in a n approximate way via a n increased permanent dipole moment. The LennardJones potential in the SPC model acts between the oxygens, which leads to the following energy expression.
The Molecular Dvnamics Method In the M D method the motions of the molecules are simulated by solving Newton's equation of motion, F = ma (and Euier's wuation for orientationai motion). The force F between tWO atom; is given by the gradient of thd potential energy, i.e., F = -(au/ax, &day, Juiaz). The forces are assumed to be constant during the time interval At. different properties may be evaluated at each time, and the results obtained are time averages. time step t = t
..........
calculate desired property: p(Q calculate the force on each atom time step t = t + At predict new positions, velocities, etc. at t + At + ... r(t + AQ = r(t)+ v(fj.At + a(fp~?/2 v(t + ~ t =)v(f)+ a(t).At+ ... a(t + At) = a(f)+ ...
...
correct positions, velocities, etc., at t + Atby using the difference between the predicted force F,,&imd(t + Af) = ma(t + A0 and the true force ~ ~ ~ l+, Af) ( t = -(au(t + A Q I ~ XaNt , + A O I ~ Ka ~ ( t A+Q I ~ Z calculate desired property: p(t+A 0 calculate the force on each atom time step t = ZAt predict new positions, velocities, etc. at t + 2At r(t+ 2At) = r(t+ Af) + ...
where the sum is taken over the nine atom-atom pairs ij in the water dimer. Ion in Water
For the LUSim system of a n ion in water eq 2 also holds, but with only three Coulomb terms for the ion-water interactions a n d a n additional LennardJones term for the ion-oxygen interaction. LUSim has default values for Na', K+, F-, and C1-, but LUSim users may define a n i o n of t h e i r choice. T h e LennardJones parameters for the default ions have been determined from quantum chemical interaction potentials (9,lO). Polymer Melt
Figure 1. Schematic description of the Monte Carlo and Molecular Dynamics methods
642
Journal of Chemical Education
A polymer is modelled by a chain of monomers connected by springs. The polymers in LUSim interact via t h e monomers a s shown below, that is, with a har-
Figure 2. The first window to appear while running LUSim. The desired system is chosen by pressing one of the "radiobuttons", in this case a LennarMones system.
Figure 4. Three different LennardJones curves showing the principal behavior as a function of E and a.
Figure 3. The input window for simulation of an ion in water with the defaultvalues forsodium. Other ions can be preselected by pressing any of the Kt, F,or Cf buttons. The parameters can also be manipulated directly in their respective input fields. monic bond potential between neighbors and with a LennardJones potential between nonneighbors. (Neighbor is defined as two wnsecutive monomers in a chain.)
I"" ireq'% 14[ [
between neighbors
-
between nomeighbom
where r,, is the equilibrium distance between the monomers; and k is a force constant.
Figure 5. Atypical display dunng the simulation of liquid argon at high temperature and pressure. that is, for64 mol. Table 1. Computation Time In S e w n d s for a Few Different Systems on Two Macintosh Modelsa System 27 argon
Time Steps/Passes
llfx
LC
200
28
724
1140
Simulation Details
125 argon
1000
LUSim, in its present version, will accept a maximum of 125 molecules, and they are initially placed in a cubic box, whose size defines the density of the system. In order to reduce boundary effects. LUSim makes use of so-called neriodic boundajwnditions: A particle leaving the simAation box a t one side enters the box from the onnosite side. The interaction between one atom (molecule)&d others is evaluated under the minimum image convention. This means that the only atoms (molecules)taken into account are those within a sphere of radius equal to half the box length from the atom (molecule)in question. The MD program makes use of Gear's predictor+orrector algorithm (71 to solve the equations of motion. The time step is dependent on the system, but 1-10 femtoseconds, s, is typical. The MC part of LUSim uses where 1 fs = the Metropolis algorithm (I).
27 HKJb
200
90
125 H20b
5000
28661
LUSim Output LUSim saves the simulation results in a number of files that can be used for calculation of static and dynamic properties for a restart and as input to graphics programs such
4'5 Polymer
100
20
1'8 Polymer
2000
200
4220 398
L ' USim runs on a Powerkmk 140 as well. '~hesetimes are also valid far the ion-water system. as KaleidaGrapha (see, for example, Fig. 6-8). I t is also possible to make animations provided one has access to a Ball&Sticka program. Exercises with LUSim In this section we describe four sample exercises with LUSim, one for each type of system, illustrating how the results are presented to the student. These might also serve as suggestions for further exercises to introduce the student to the most basic concepts of liquid state physics. A word of warning: LUSim is not meant as a program package for students with no knowledge of statistical mechanics. Rather, it is designed to be a laboratory tutorial at Volume 70 Number 8 August 1993
643
Figure 6. a (top,left). The approach towards equilibrium for liquid argon starting from a primitive cubic lattice. b (bottom,left). The radial distribution function for liquid argon at 94.4 K from a simulation of 64 atoms. The hardly visible dashed line was obtained with 27 argon atoms. c (top,right). The velocity autowrrelation function for liquid argon at 94.4 K from a simulation of 64 atoms. d (bonom,right). The mean square displacement of argon atoms as afunctionof time at three different temperatures and constant density, 1.37 gicm3. the end of an introductory course in statistical mechanics. The examples below are constructed so that the student should be able to work essentially without supervision. Computation Time
One important aspect is the computation time needed for a reasonablv sized simulation. A simulation of s~hericallv symmetric rare gas atoms interacting via a IJ pdtential is, of course, much faster than a simulation of liquid water. simulations of long polymers are also quite tiie-consuming. The original simulation programs in the LUSim package have been optimized for vector computers, and we have made no attempts to reoptimize for the Macintosh. All our effortshave gone into making it as user-friendly and easily accessible as possible. The total length needed for the complete set of examples is about 4 h. Table 1 show typical execution times on two different Macintosh models. The simulation exercises were carried out on either a model LC or IIfx. LiquidArgon
Choice of Starting Configuration In this exercise we attempt to model liquid argon by L e ~ a r d J o n e particles. s An instructive start is to run a very short simulation of say 27 argon atoms lasting for 0.05 ps starting from a primitive cubic lattice. It will be immediately clear that this time is much too short to reach 644
Journal of Chemical Education
equilibrium, but the exercise gives the student an opportunity to follow how the internal enerm -. decreases towards its kquilibrium value. The choice of starting configuration is to some extent arbitrary, and the simulation could also be started from a random confiwration. None of these two initial structures are good representatives of liquid argon, and we must allow the system to equilibrate. It is only aRer this equilibration period that we can expect to obt& meanin&l averages in a so-called production simulation. Approach towards Equilibrium Different quantities are displayed on the screen during the simulation. includine the radial distribution function and an animation of theatomic motion (see Fig. 5,. The initial simulation can be followed bv another that is loneer by a factor of 10. Then an approximate equilibriumis reached. (Figure 6a shows how the internal energy changes during these two simulations.) The approach towards equilibrium can also be studied by looking at the radial distribution function displayed in the lower left corner during the simulation. A final 2-ps production simulation may conclude the exercise, which altogether takes less than 40 s on a model IIk. A more realistic simulation, which actually gives numbers very close to experimental results, can be carried out with 64 argon atoms starting with an equilibration run lasting 2 ps followed by a production run of 10 ps. The radial distribution function from the latter part is shown in
Figure 6b, which shows the first coordination shell, typical fi~rall dense liauids. and ulso the initial osrt of the second more diffuse iayer. The curve is actually in very good agreement with the classical simulation of liquid argon by Rahman (3).F i r e 6b also contains the radial distribution function from a 10-ps simulation of 27 argon atoms, which is hardly distinguishable from the larger simulation. Cage Effect The velocity autocorrelation function obtained in the analysis section is shown in Figure 6c. A small minimum occurs afier the initial rapid decay, which takes place in less than a picosecond. This "cage effect"ofthe velocity correlation function is typical of a dense liquid with frequent collisions. The translational diffusion coefficient is related to the integral of the velocity autocorrelation function, but it can also-be calculated f
