Army Research Laboratory 

Aberdeen Proving Ground, MD 2 10055066 

ARL-m-2389 February 2001 

Transitioning Model Potentials 

to Real Systems 

Scott D. Bembenek and Betsy M. Rice 

Weapons and Materials Research Directorate, ARL 

Approved for public release; distribution is unlimited. 
Abstract 

The parameters of two pairs of potentials that describe argon over its entire liquid phase at a 

fixed pressure were optimized through a novel application of constant temperature and pressure 

molecular dynamics (NW-MD) and Monte Carlo (NPT-MC) computer simulations. The forms 

of these potentials were those of a modified Lennard-Jones potential and a Lennard-Jones 

potential (Lennard-Jones, J. E. Physicn, Vol. 4. p. 941, 1937). The optimized potential 

determined using NPT-MD simulations reproduced experimental densities, internal energies, and 

enthalpies with an error less than 1% over most of the liquid range and yielded self-diffusion 

coefficients that were in excellent agreement with the experiment. The results using the potential 

determined by WY-MC simulations were in almost as good agreement with deviations from 

experiments of no more than 5.89% for temperatures up to vaporization. Additionally, molar 

volumes predicted using this potential at pressures in the 100-600 atm range and over 

temperatures in the loo-140 R range were within 0.83% of experimental values. These results 

show that when properly parameterized, Lennard-Jones-like potentials could describe a system 

well over a large temperature range. Further, the method introduced was easy to implement and 

was independent of the form of the interaction potential used. 

ii 
Acknowledgments 

This research was conducted while Scott D. Bernbenek held a National Research Council 

Resident Research Associateship. 
iv 
Table of Contents 

Acknowledgments ._._ .___. f.. ..f.......f............_._ ..______.f..f_.............f...__ .___f.._..........f...--. . 111 . . 

List of Figures . . . ..f.........._______________._____..__ + .f_~...............f._._.._._...__...f.~.~ff~.~........._.__.__._ vii 

List of Tables . . . ..f.f..................................-..-.-.--....-.-...-.......................--.---..-....*..-...... ix 

1. Introduction .~.. .f........_.f_______._._____f.___............~~f.......f......,,.....______ __f_._..f...._.............. 1 

2. Simulation Method . ..f..........f~.~....f.........._..I..._..._.........................f....._._...__._f_.____.....f 2 

2.1 Molecular Dynamics (MD) .................................................................................. 2 

2.2 Monte Carlo (MC) ............................................................................................... 6 

3. Results ._.__ ._. ____ . ... . _ ...~.ff.f....*..f*.**..*...*...~........~~.......................~~~.~.~....*..*.*..**..~.....~.~.* 8 

3.1 Molecular Dynamics (MD) .................................................................................. 9 

3.2 Monte Carlo (MC) ............................................................................................... 13 

4. Conclusion -..... . . . . . . . . . . .._.f__ff.....f...._..____.______. . _____ _____. ._._............f..........,.ff.____.._ ,._____ 17 

5. References _____+.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..f_f......___________.__-.--......................_.f_ 19 

Distribution List . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..-.-...---.--..-.-.--..............................._... 23 

Report Documentation Page __.____f____._.....................................,,.........--....--.---.---..... 25 

V 
INTENTIONALLYLEFTBLANK. 

. 

vi 
List of Figures 

Figure & 

1. Value of < E: >/k and Corresponding 0 for Liquid Argon at a Pressure 

of 40 bar for Several Temperatures _______._f.. _ .. . . . . .._._________. + ...ff...._._______ff......f.._.__.._..ff.f 10 

2. Value of <&B/k for o = 3.40 A Plotted as a Function of the Cutoff Distance rcut 

(in 0) for Liquid Argon . . ..-_. f..**.........._._.__...._.*....f..........__.....~...ff.f.~~~___......f......~....~f~~ 11 

3. Number Density of Liquid Argon (in Atoms A-‘) as a Function of 

Temperature . . . _ .. __ _ _ _*. . . . . . . ._ .__ __ __* . . . . . . . . . . . . . _*. ._ _ _ ._ .. . . . . . . . . . _ _ _ __ . +. . . . . . . . . . . . . _ . _ _ _ . . . . . . . . . . . ._ .. . 13 

4. Internal Energy (in eV Atom-‘) as a Function of Temperature .._.___._._........_......_....... 14 

5. Enthalpy (in eV Atom-‘) as a Function of Temperature . . . ..____.. .f........._____._...........__ 15 

6. Value of < E >/k and Corresponding G for Liquid Argon at a Pressure 

of 40 bar for Several Temperatures.. ..-..._..... _- + .. . . . . . ..___ ..__.-....... ..,.. . . . . .._.....+......- 16 

vii 
. . . 

Vlll 
Table 

1. 

2. 

3. 

4. 

List of Tables 

The Value of <~>/k for the Common Set as a Function of Cutoff Distance 

Using Equation 1 ff..........f_______..............,..______._f......f.f_._______.._._.f.f....._._f______________f_

10 

Density, Internal Energy, and Enthalpy of Liquid Argon . . . . . . . . . . . . . . .._f.__________............. 12 

Self-Diffusion Coefficients (in cm2-s-‘) for Liquid Argon .____.___.._._ff.....fff...._..__...____ 15 

Molar Volume of Liquid Argon (in cm3 mol-‘) . . . . . . . . . ....f_.___.______ ...__._f..f..f.....f_. . . . . ..___ 17 

ix 
X 
1. Introduction 

One of the major goals in any computer simulation is to accurately predict the properties of 

the system. This can only be accomplished with the proper choice of interaction potential. The 

interaction potential is a fundamental feature of a system. All thermodynamic and transport 

quantities will be determined by its form. Unfortunately, analytical functions, such as the .I 

hard-sphere potential, square-well potential [ 1 J and the Lermard-Jones potential [2 3, provide only 

an approximate description of the underlying interaction potential. The Lennard-Jones potential 

has played an important role in the theory of classical fluids. Its simple form makes it desirable 

for computer simulations and theoretical calculations alike. It has been extensively used in 

computer simulations of liquids [3-S], glasses [(i-12], and phase coexistence [13-261, to name 

only a few. It has also been widely used as a reference fluid in perturbation treatments for more 

complex fluids [27]. Further, several attempts have been made to obtain an equation of state; the 

most notable being by Nicolas et al. [28] and another, and more recently by Johnson et al. [29]. 

The Lennard-Jones potential requires only two parameters in its description: E, which is the 

depth of the potential well, and CT, which is the effective particle diameter. Often, there are many 

parameter pairs that will predict with equal accuracy a thermodynamic property at a given phase 

point. However, the potentials using these different parameter sets might produce widely 

varying predictions at other phase points. In principle, there should be a single parameter set that 

will best describe the system at every phase point within the limitations of the model. In this 

study, this set is referred to as the %ornmon set.” This report will describe the method of fmding 

the common set for a modified Lennard-Jones potential and a Lennard-Jones potential using 

either molecular dynamics or Monte Carlo simulations. It will be shown that a potential using 

the common set reproduces features of the system over the entire range studied. The LennardJones form of the interaction potential can be regarded as one member of any number of 

interaction potentials that could describe the system properly over the entire liquid phase and, 

therefore, is not unique for parameterization using this method. In principle, this approach can 
be used to optimize any functional form of the interaction potential. The method will be 

illustrated using liquid argon. 

2. Simulation Method 

The method presented can be tested by performing either Monte Carlo (MC) or molecular 

dynamics (MD) simulations under constant temperature and pressure. 

2.1 Molecular Dynamics (MD). In MD simulations [30, 311, driving a system into a 

nonequilibrium steady state by coupling it to the appropriate external fields is well described 

[32]. These nonequilibrium molecular dynamics (NEMD) simulations allow for faster and, in 

some cases, more accurate determination of certain transport quantities than equilibrium 

simulations. A currently unexplored concept is presented here, where a nonequilibrated liquid 

system at constant temperature, pressure, and particle number is driven to a steady state under 

the constraint of achieving a particular bulk property. For this study, the bulk density was 

chosen to be the constraining “external field.” Here, the external field will continuously “push” 

the system towards the desired bulk density. For these purposes, the desired density is the 

experimental density of liquid argon. 

Consider a nonequilibrium system of N particles at constant temperature, T, and pressure, P, 

interacting through a modified Lennard-Jones potential, 

v(r) = 4~ [(g -(3”]+M +B, 

where A and B are constant and chosen such that both the potential and the force are continuous 

at the potential cutoff, rcUf. Modifying the potential in this way ensures that the energy is 

conserved. 
The bulk density constraint is imposed by starting the system with initial values of E: and 6 

and allowing a to vary in time, t, according to 

f(t)= E. +&A&), 

t’=l 

(2) 

where h is the coupling strength, EO is the initial value (seed) for E and the difference between the 

desired and instantaneous bulk densities is Ap(t’) = o,jesjd - p(t). The effect is that E will be 

made larger, thereby increasing the attraction between particles, if the system is ?oo dilute.” 

Similarly, E will be made smaller (thus decreasing the attraction) if the system is Yoo dense.” A 

value for h that is too large will cause the system to dramatically overshoot the desired density, 

thus requiring a longer amount of time for the system to relax under the “external field.” The 

system momentum will be conserved, as in a usual equilibrium simulation, since the external 

field is actually coupled into the pair potential itself. Failure to conserve momentum may be 

further indication that the value of h is too large. If h is too small, the external field will be 

unable to compete with the fluctuations inherent in the system and, therefore, the system will 

never approach the desired bulk density. 

f 

Eventually, the instantaneous densities of the system will oscillate around the desired bulk 

density. It is necessary that the system oscillate several times around the desired bulk density 

before E is chosen. This ensures that the E that is chosen will yield the desired bulk density. In 

the case of liquid argon, it was found that the integration of 25,000 time steps (1 time step 

= 5.1 fs) was needed for the system to approach an E that gave the desired density. After this 

time, an E that reproduced the desired bulk density (for a given time step) to within a certain 

amount was chosen. This amount can be defined with any precision and was chosen here to be 

5 x lo4 particles A-‘. The ability of the system to come within this established amount depends 

on the natural fluctuations of the system. 

3 
Repeating this procedure over a range of o’s for a given temperature and pressure will result 

in a set of several combinations of E and 0’ that will reproduce the desired density for this phase 

point. The cornmon set, which refers to the parameter set that reproduces the density over all 

temperatures in the liquid phase at a given pressure, requires generating ensemble averages of s’s 

for the various o’s at several temperatures over the range to be studied. With the exception of 

simulations for which rcut = 6.00~, each ensemble consisted of ten trajectories; the remaining 

ensembles consisted of five trajectories. Once these parameter sets ‘have been generated, the 

average value of E, < E > , for each G is calculated and curves of < E > vs. 0 for each temperature 

are drawn. The 0 of the common set corresponds to the point where all the curves intersect 

(within statistical error). The E of the common set is the average value of < E > determined for 

all the temperatures. 

The common set obtained for the modified Lennard-Jones function in equation 1 is 

dependent on the value of rcut. As rcut becomes large, however, the changes in the common set 

become less pronounced, as will be shown. In order to examine this dependence on the value of 

rCUt and satisfy the minimum image convention imposed in this simulation it was necessary to 

increase the system size as reut increased. The systems consisted of 500; 864; 1,024; and 

2,000 particles for r Cuf’~ of 2.5G3.500 (0.02% intervals) 4.00~, 5.000, and 6.000, respectively. 

Particle masses were set to 39.95 an-m, and all simulations were run at a constant pressure of 

40 bar. The pressure was held constant using Andersen’s method [33], making it necessary to 

choose a piston mass. A larger mass will result in better momentum conservation than a smaller 

one, but the system will explore volume space more slowly. In an effort to satisfactorily 

conserve momentum and sample volume space in a reasonable amount of time, the mass was 

chosen to be 2 x 10S3 amu A+ for the systems consisting of 500 and 864 particles, and 5 x 10” 

amu AA for the larger systems Finally, the temperature was held constant through a simple 

scaling of velocities, and the equations of motion were integrated using ,a modified 

velocity-Verlet algorithm [30, 34, 351, which eliminated the need for scaled coordinates. The 

coupling strength, h, was chosen to be 1 x lOa eV A3 particle-‘. 

4 
Average values for E were determined using the procedure previously described and 

equation 1 with the different values of rCUt. Each ensemble, except for the series with 

r cut - - 6.OOo, consisted of 10 trajectories; the ensemble for rcut = 6.00~ contained 5 trajectories. 

The initial conditions were selected as follows: for each value of rCuf, o is set to the minimum 

value in the range of o’s to be investigated and s is assigned an arbitrary value. The particles 

were placed in lattice sites of either a bee or fee crystal to generate the initial coordinates. Initial 

velocities were set to zero, and the coordinates slightly displaced” from the equilibrium lattice 

sites. The equations of motion were integrated for 50,000 time steps at a temperature of 135 K, 

and a pressure of 40 bar. This was found to be sufficient for equilibration. Once the system 

equilibrated, the equations of motion were then integrated for an additional 75,000 time steps, 

but E was allowed to vary according to equation 2. The trajectory was stopped, and the final 

value of E (denoted ~0) and the coordinates and velocities were used as initial conditions for each 

trajectory with that value of r EUt over all the values of o and temperature studied. It is important 

that a suitable EO is obtained so that the system will not prematurely “vaporize” (the volume of 

the simulation box will expand without bound during an APT simulation). The problem is 

circumvented by performing the equilibration of the system near the boiling point (135 K) and 

using the smallest value of CT (in this case, 3.35 A) in the range. This results in the shortest 

attractive tail for the interaction potential (since r CUt is chosen to be a function of 0) over the 

range of o’s. 

Some care must be taken not to put the system too close to the boiling point, otherwise, the 

system might begin to make the phase transition, which would result in a poor value for ~0. Each 

trajectory to be used for < E > was integrated for 25,000 time steps beginning with the initial 

conditions previously described and allowing E to vary according to equation 2. The 25,000 time 

steps were sufficient for the system to adjust to each new value of o. By this time, the density of 

the system was oscillating about the desired value. After this period, the instantaneous density 

was monitored. If the instantaneous density was within 5 x 10s6 of the experimental value at that 

temperature and pressure, the trajectory was stopped, and the value of E at that point was selected 

as the correct value for that o value. Once this entire procedure has been performed for 

5 
the series of o’s at a single temperature, the process is repeated for that temperature, except the 

initial conditions for the subsequent run (coordinates and velocities) are those corresponding to 

the immediately preceding trajectory. This is repeated until the desired number of values of E’S 

that will be used in averaging are calculated for a given temperature. The G’S that were sampled 

ranged from 3.35-3.45 A (in increments of 0.025 A) and the temperatures sampled ranged from 

85-145 R (in intervals of 5 K). Once the set of s’s and G’S have been determined for all of the , 

temperatures, the cormnon set can be obtained by superimposing pl&of the curves showing 

< E > vs. G for all temperatures sampled in the study. 

In order to evaluate the performance bf the modified Lennard-Jones potential using the 

common set obtained for rcUt = 6.00~, time averages of bulk properties of the liquid over the 

temperature range 85-145 K were calculated. For each temperature, a trajectory consisting of 

250,000 time steps was integrated to generate the ensemble averages. The initial coordinates and 

velocities of a system of 2,000 particles were generated by first setting up the particles in an bee 

lattice, and performing a trajectory integration for 50,000 time steps at the desired temperature 

and pressure, using equation 1 and rcut = 6.000. The ensemble averages were obtained by 

averaging over the calculated values for the remaining 200,000 time steps. 

2.2 Monte Carlo (MC). This technique can be implemented using MC simulations as well. 

MC does not require continuous forces at the cutoff distance; therefore, this method will be 

demonstrated using the unmodified Lennard-Jones potential. Long-range corrections &e 

included to account for all interactions of pairs whose distances exceed rcut [30,31]. The density 

constraint can be implemented in the same way as was done for MD using equation 2. However, 

the time dependence in equation 2 is not applicable to MC, and therefore, it is replaced with a 

“step” dependence, 

step 

Estep=Eg+XzA&, 

i=I 

6 

(3) 
where h is the coupling strength as previously described, sStq is the value for E for that step of the 

Markov sequence, EO is the initial value for E and Api is the difference between the desired and 

instantaneous bulk densities. 

Parameter sets for various values of o over the temperature range 85-145 R were determined 

using NIT-MC simulations and equation 3. As in the MD simulations, these sets were 

determined for o values ranging from 3.35-3.45 A in increments of 0.025 .A,. The simulation 

box contained 150 atoms, placed randomly within the box, and the size of the box was initially 

set such that the initial density was equal to the experimental value at 135 K, 40 bar. As in the 

MD simulations, 0 is set to the minimum value in the range to be studied, and E is given an 

arbitrary value. The NPT-MC simulation proceeded for 100,000 moves, with E remaining fixed. 

This was sufficient to equilibrate the system to the desired thermodynamic state. Imposing the 

same temperature and pressure constraints, the NIT-MC simulation proceeded for an additional 

100,000 Markov moves but the E was now allowed to vary according to equation 3. The 

sequence continued with a trial move of the particle positions, which generated an instantaneous 

density. A new value of E was generated according to equation 3. The energy of the system 

using the new value of E was calculated, and the changes accepted or rejected according to the 

probability min [ 1, exp (- W/kT) J [30,3 l] where 

w = p(Lv - L ) + (vnew 

V -U,)+NkTlnL, 

V (4) 

new 

and P, V, and U denote the pressure, volume, and potential energy of the system, respectively. 

The magnitude of the displacements of the particles during the Markov sequence were chosen 

such that 50% of all attempted moves were accepted. If a trial move was rejected, the properties 

associated with the immediately preceding configuration were included for averaging. 

The final value of E at the end of the 100,000 moves became EO and the resulting coordinates 

were used for the initial conditions for NE”-MC calculations to determine E as a function of o 

over the entire liquid phase at 40 bar. As in the MD simulations, this process is repeated for the 

7 
series of o’s at each temperature, except that the initial conditions for the subsequent run 

(coordinates) are those corresponding to the final configuration of the immediately preceding 

Markov sequence. This is repeated five times in this study. The common set is determined, as 

before, by superimposing the plots of < E > vs. o for each temperature. 

In order to assess the performance of the Lennard-Jones potential using this set of parameters 

and those derived from other studies [36, 371, NPT-MC simulations of liquid argon at 40 bar 

over the temperature range 85-145 K were performed. The initial system consisted of 108 atoms 

in a cubic box where each edge length was 17.46 A. Periodic boundary conditions were 

imposed, and a cutoff distance equal to one-half of the edge of the box was used. Long-range 

corrections to the calculations were included. An initial equilibration simulation using the 

Lennard-Jones parameters of Wijker et al. [37] for T = 97 K, P = 17 atm was performed in order 

to reproduce one of the values generated by McDonald and Singer [36]. An excess of 500,000 

Markov moves were used for equilibration, and 100,000 moves were used in averaging the 

results. The results using this parameter set were reproduced. This coordinate set was then used 

as the initial set for the calculations of properties at X5 K using parameters obtained in this work, 

from McDonald and Singer [36], and from Wijker et al. [37]. An equilibration Markov sequence 

of 10,000 steps was performed at T = 85 K, P = 40 bar and thermodynamic properties were 

calculated and averaged over the next 100,000 steps. A new temperature was then selected, and 

the process repeated until the entire temperature range had been sampled. Longer Markov chains 

were generated, and different sequences were initiated, but the results did not vary. 

3. Results 

Although argon has been well described using more complex interaction potentials [38, 391 

the purposes of this study are (1) to illustrate the newly described method for obtaining 

parameter sets using either MD or MC, and (2) to show that the Lennard-Jones potential can 

obtain accurate results over the entire liquid phase upon proper parameterization. 

8 
. 

3.1 Molecular Dynamics (MD). Parameters for the modified Lennard-Jones potential 

(equation 1) were determined using the method described in the preceding section for 

temperatures of 85-145 H. The value of c&>/k and the corresponding o for liquid argon at a 

pressure of 40 bar at several temperatures is shown in Figure 1. This is presented for cutoffs 

rcut’s of 2.500 and 6.00~. Figure l(a) shows c&>/k vs. o using rcUt = 2.500, and Figure l(b) 

shows the values using rcut = 6.000. In the figure, the common set corresponds where all of the 

curves intersect (within statistical error). This changes considerably from’ 2.50~ until about 

4.000, where it slowly converges to its value at 6.00~. The common set at 6.00s was used to 

calculate all properties of interest for the MD simulations. The error bars for all points are on the 

order of the size of the symbols except for those points at 145 K, which show the error bars 

explicitly. The deviations at 145 K are due to the proximity of the system to the boiling point; 

the system is beginning to make the expected phase transition. Since the values generated at 

145 K have such greater uncertainty than the remaining points, these have not been included in 

the determination of the common set. For all values of rCUf, the G in the common set is 3.40 A but 

the corresponding E values are very different. For example, for rcut = 2.500, the value of <&k is 

155.876 K, but the value of <E>/k resulting from the simulations in which rcUt = 6.000 is 

119.808 K. The changes in the value of the E of the common sets with increasing values of rcut 

are given in Table 1, and illustrated in Figure 2. The values of E for the common set appear to be 

converging as rcut increases (i.e., <aHk converges for 0 = 3.40 8, to = 120 K as rcut is extended). 

The error bars have been eliminated since they are comparable with the size of the symbols, 

The performance of the modified Lennard-Jones potential using the common set obtained for 

rcut = 6.00~ was tested, in order to keep the system size small enough that molecular simulation 

would not be computationally prohibitive. NPT-MD simulations were used to calculate 

ensemble averages for density, internal energy and enthalpy over the temperature range 

85-145 K at 40 bar. Table 2 gives the calculated ensemble averages for densities, internal 

energies and enthalpies along with the experimental values [40], and the comparisons are 

illustrated in Figures 3-5, each of which are based on MC calculations using the parameters of 

Wijker et al. [37], McDonald and Singer [36], and those obtained in this work along with the 

9 
0 

I 

Y 

T x 

0 A I 

<E> 

0 120 160 

1 

P 

B 

- q T=95K 

x T=105K 

AT=115K - 110 

+ T=125K 

VT=135K 

* T=140K 

OT=145K 

140 3.34 ' 3.38 ’ 

3.42 

100 -. 3.34 3.38 3.42 3.46 

(4 @I 

Figure 1. Value of < E >/k and Corresponding [T for Liquid Argon at a Pressure of 40 bar 

for Several Temperatures. 

Table 1. Value of <&B/k for the Common Set as a Function of Cutoff Distance Using 

Equation 1 

155.876 -- 

117.653 144.707 

Deviation 

0.161 

0.126 n 11< 

I 

, .--- 

I33.112 

“.A I4 

3.25 0.171 

7 tin 13OQ7< n IO" 

5.00 

6.00 

aThe value of o for the common set is 3.40 A. 

*--. .I I v.1 I” 

121.516 0.069 

119PQ ,“” I 

n no0 

“.“OO 

experimental values [40]. Also shown in Figures 3-5 are the values obtained from MD 

simulations using the parameters found in this work. The error bars have been eliminated since 

they are comparable to the size of the symbols. The density (Figure 3), the energy (Figure 4), 

10 
165 T 

155 

, 145 

<E> 

k 

125 

ll5-. 2.0 3.0 4.0 5.0 6.0 

rcut 0 ( > 

Figure 2. Value of c&>/k for (5 = 3.40 A Plotted as a Function of the Cutoff Distance rrut 

(in a) for Liquid Argon. 

0 

0 

] 

0 I 

0 1 

0 

0 

1 

and the enthalpy (Figure 5) have been multiplied by lo2 for convenience. The calculated results 

are within 1% of the experimental data for most of the results over the temperature range 

investigated. (Once again, the deviations at 145 K result from the system being at the boiling 

point.) These results are comparable with results that were obtained with a more complicated 

interaction potential [39]. 

The ability of the modified Lennard-Jones potential using the common set obtained at rcut 

= 6.000 to describe the self-diffusion coefficient has also been evaluated. In their paper [41], 

Naghizadeh and Rice present experimental data of self-diffusion for liquid argon over a 

temperature range of 9C120 K for pressures of 12.9, 57.5, 103.0, and 135.0 atm. The behavior 

of the self-diffusion coefficients with temperature at each pressure was satisfactorily described 

using an Arrhenius-like expression. Naghizadeh and Rice fitted their data to this expression, and 

provided an analytical description of the temperature dependence of the self-diffusion coeffcient 

for liquid argon could be estimated at a pressure of 40 bar for each temperature, to allow for a 

comparison with the MD simulations results. 

11 
Table 2. Density, Internal Energy, and Enthalpy of Liquid Argon 

3/l@ atoms A-’ -U/l@ eV atom-’ -H/lP eV atom-* 

This work This work This work 

r/K Expt. McDonald Wijker MC MD Expt. McDonald Wijker MC MD Expt McDonald Wijker MC MD 

and et aLb and et al.' and et alb 

Singer’ Singer’ Sing& 

85 2.1388 2.1495 2.1667 2.1513 2.1334 5.024 5.094 5,296 5.022 5.048 4.907 4.978 5.181 4.906 4.931 

90 2.0933 2.0976 2.1169 2.1039 2.0877 4.807 4.844 5.048 4.789 4.821 4.688 4.725 4.930 4.670 4.701 

95 2.0467 2.0475 2.0681 2.0541 2.0413 4.584 4.602 4.804 4.551 4.593 4.462 4.480 4.684 4.429 4.471 

loo 1.9989 I.9987 2.0257 2.0041 i.9939 4.357 4.368 4,587 4.316 4.364 4.232 4.242 4.464 4.191 4.239 

105 1.9488 1.9477 1.9680 1.9501 1.9443 4.124 4.129 4.323 4.069 4.132 3.996 4.001 4.196 3.940 4.003 

110 1.8966 1.8930 1.9232 1.8972 1.8933 3.887 3.881 4.101 3.829 3.897 3.754 3.749 3.971 3.697 3.765 

115 1.8410 1.8368 1.8640 1.8336 1.8382 3.643 3.632 3.838 3.562 3.652 3.507 3.4% 3.704 3.425 3.516 

120 1.7816 1.7786 1.8151 1.7750 1.7806 3.392 3.380 3.609 3.309 3.403 3.252 3.240 3.471 3.168 3.263 

125 1.7167 1.7158 1.7491 1.7126 1.7199 3.130 3.120 3.329 3.052 3.149 2.985 2.974 3.186 2.906 3.004 

i30 1.6442 1.6414 1.6840 1.6416 1.6514 2.855 2.828 3.060 2.772 2.875 2.703 2.675 2.911 2.619 2.724 

135 1.5601 1.5708 1.6146 1.5539 1.5719 2.551 2.550 2.779 2.453 2.575 2.391 i.390 2.624 2.292 2.416 

140 1.4567 1 A683 I.5358 1.4490 1.4781 2.208 2.190 2.474 2.090 2.243 2.037 2.020 2.311 1.917 2.074 

145 1.3079 1.1768 1.4237 0.8719 1.3520 1.770 I.356 2.086 0.542 1.837 1,579 1.091 1.909 0.114 1.652 

1361 :. 

b[371 

r 
P . 

1.6 - 0 Experiment 

0 McDonald and Singer @. 

V Michels et al. 

1.4 0 This Work-MC - 

t This Work-MD 

gT7_ 

* 

q . 

1.2 

60 90 100 110 120 130 140 160 

T(K) 

Figure.3. Number Density of Liquid Argon (in Atoms Am3) as a Function of Temperature. 

The self-diffusion coefficient from MD is calculated from the mean-square displacement and 

is an average over ten trajectories each consisting of 20,000 integration steps. Table 3 compares 

the experimental and calculated self-diffusion coefficients as a function of temperature. The 

calculated values have a larger percentage deviation (with a maximum error of 8%) from 

experimental values than those of the thermodynamic properties, but are still in reasonably good 

agreement. However, Naghizadeh and Rice indicate that their data contains uncertainties of less 

than 5%, which could account for some of the discrepancies. 

c 

3.2 Monte Carlo (MC). Lennard-Jones parameters for liquid argon were obtained using the 

method described in section 2 and APT-MC calculations over the temperature range 

85-145 K, at 40 bar. Figure 6 shows < E b/k over the o range .3.35-3.45 8, in increments of 

0.025 8, obtained from NPT-MC calculations using this method. In this figure, the common set 

corresponds to o = 3.40 8, as it did in the application of the method using NPT-MD. The value 

of E, averaged over all temperatures is < E >/k = 116.359 + 0.128 K. The value for E is 3% 

smaller than that obtained by Wijker et al. [37] and is within the error of the value given by 
6 

5 

4 

-u 

3 

2 

1 

V 1 

I3 v 

@ ,v 

8 v 

b v 

8 v 

d v 

80 

80 

a Experiment 

0 McDonald and Singer 

V Michels et al. 

0 This Work-MC 

+k This Work-MD 8vti 

0 I 

60 90 100 110 120 130 140 150 

T(K) 

Figure 4. Internal Energy (in eV Atom-‘) as a Function of Temperature. 

McDonald and Singer (117.2 K f 1.4 K) [36]. The value for o that was calculated is close to 

that of Wijker et al. [37] and McDonald and Singer (3.40 A) [36]; however, when generating the 

curves of E as a function of d at each temperature, those for a fine grid of 0 values were not. 

Rather, these points were generated at intervals of 0.025 A. 

The parameters of Wijker et al. [37] were obtained using experimental gas-phase information 

for argon, rather than liquid state information. SmaJl discrepancies between experiment and 

results from constant volume and temperature (MT)-MC -simulations using these parameters to 

describe thermodynamic properties in the liquid state led McDonald and Singer to propose a new 

set of parameters that would better reproduce values in the liquid state [36]. In order to compare 

the results &scribed herein with those of experiment .(Wijker et al. [37] and McDonald and 

Singer [36]), NPT-MC simulations were performed to predict densities, internal energies and 

enthalpies for liquid argon at 40 bar, over the temperature range 85-145 K using Lennard-Jones 

potentials with the parameters of McDonald and Singer [36], Wijker et al. [37], and the set that 

14 
6 

5 

4 

-H 

3 

2 

1 

I I 

t 

$ v 

6 v @I v I 

I B 7 

b v I 

8 v 877 1 

0 Experiment 

McDonald and Singer 

80 

0 V Michels et al. 80 

1 

0 This Work-MC $90 j 

N This Work-MD doi 

Figure 5. Entbalpy (in eV Atom-‘) as a Function of Temperature. 

Table 3. Self-Diffusion Coefficients (in cm’+-‘) for Liquid Argon 

T/K Dexp H 1 O5 DMDH lo5 

85 1.681 1.783 

90 2.122 2.281 

95 2.614 2.709 

100 2.954 3.191 

105 3.737 3.764 

110 4.361 4.533 

115 5.022 5.374 

120 5.71 s 5.950 

15 
130 

<E> 

k 

r 

0 0 T=85K 

+ 0 T=BSK 

a 0 x T=105K 

X AT=llSK 

E + T=l25K 

0 

VT=135K 

m T=14OK 

0 T=145K 

3.35 3.42 3.46 

Figure 6. The Value of < E >/k and Corresponding 0 for Liquid Argon at a Pressure of 

40 bar for Several Temperatures. 

was derived using the method described in this work. The results of all calculations are given in 

Table 2 for comparison with experiment and among the predicted values; Figures 3-5 provide 

visual comparisons. With the exception of values at 145 K, the calculations using the Wijker et 

al. [37] set are in poorer agreement with experiment than potentials using the set derived herein 

or the McDonald and Singer set [36]. Properties calculated at 145 K using both the McDonald 

and Singer set and the parameters determined in this work have a significant deviation from 

experimental values. The values, however, do not represent an average over a single phase since 

the system is actually fluctuating between the liquid and vapor states [30], which was evident 

upon inspection of the instantaneous densities. At temperatures away from this phase transition, 

the properties predicted using either McDonald and Singer or our set are comparable, thus 

affirming McDonald and Singer’s conclusion that the parameters generated by Wijker et al. 

produce an attractive well-too deep for liquid argon [36]. F 

16 
The parameters generated in this work were determined using experimental information of 

the density of liquid argon at 40 bar; in order to assess the performance of this potential at other 

pressures, APT-MC calculations were performed at pressures ranging from lo&600 atm, at 100, 

120, and 140 K for comparison with experiment. The resulting molar volumes are given in 

Table 4 and compared with the experimental results of Streett and Staveley [42]. The predicted 

values over the entire pressure and temperature range are in agreement with experiment to 

0.83%. 

Table 4. Molar Volume of Liquid Argon (in cm3 mol-‘) 

4. Conclusion 

It has been shown that it is possible to optimize a pair-additive interaction potential using 

molecular simulations in which the system is constrained to reproduce the experimental bulk 

property over a range of temperatures at a constant pressure. Specifically, this was shown for 

liquid argon using either a modified or unmodified Lennard-Jones potential, and using the bulk 

density of the liquid as the constraining factor. The modification to the Lennard-Jones potential 

used in the APT-MD simulations was made to conserve the energy during the trajectory 

integration and to ensure the force was continuous at the cutoff distance. Since APT-MC 

17 
simulations do not require this, application of the method using NPT-MC was performed using a 

regular Lennard-Jones interaction potential with cutoff. 

NPT-MD simulations using this method produced values of <E >/k = 119.808 K and 

o = 3.40 A with rcut of 6.00~ were found to be the best set of parameters over the temperature 

range of 85-145 K. Using these parameters, excellent agreement with experiment for density, 

internal energy, enthalpy, and self-diffusion coefficient was obtained. This suggests Lthat density 

is a good constraint. NPT-MC simulations using this method produced values of <E >/k 

= 116.237 K and (J = 3.40 8, for the unmodified Lennard-Jones potential, which also gave results 

that were in very good agreement with experiment over the temperature range. 

Results from the simulations using either NPT-MD and A?PT-MC were compared against 

available experimental information and NPT-MC calculations using parameters for the 

Lennard-Jones potential suggested by McDonald and Singer [36] and Wijker et al. [37]. The 

results using the parameters generated in this study were in better agreement with experiment 

than those using the Wijker et al. set, and were comparable to those generated using the 

McDonald and Singer results. Additionally, NPT-MC calculations were performed using the 

common set of parameters determined for the Lennard Jones potential for temperatures ranging 

lo&140 K and pressures of 100-600 atm for comparison with experimental results. The molar 

volumes calculated using NPT-MC were within 0.83% of the experimental values. The results 

presented herein seemed to indicate that E and o are independent of temperature and pressure. 

Although this method should be applied to other systems to determine if it will show the 

same degree of success in parameterizing simple potential energy functions, the results herein 

suggest that this might be a very useful and powerful tool for easily and effectively developing 

model interaction potentials to describe real systems. Such models can then be used to obtain 

thermodynamic and some transport properties from computer simulations. 

18 
5. References 

1. Hkschfelder, J. O., C. E CL&S, and R. B. Bird. Molecular Theory of Gases and Liquids. 

New York, NY: Wiley, 19@. 

2. Lennard-Jones, J. E. Physica. Vol. 4, p. 941, 1937. 

3. Verlet, L. Physical Review. Vol. 159, p. 98, 1967. 

4. Buchner, M., B. M. IAanyi, and R. M. Stratt. Journal of CkmicaZ Physics. Vol. 97, 

pm 8522, 1992. 

5. Larsen, R. E., G. Goodyear, and R. M. Stratt. Journal of Chemical Physics. VoL 104, 

p. 2987,1995. 

6. Madsen, J. U., and R. M. J. Cotterill. Physics Letters, Vol. 83A, p. 219, 1981. 

7. Keyes, T. Journal of Chemical Physics. Vol. 101, p. 5081, 1994. 

8. Vollmay, K, W. Kob, and K Binder. Europhysics Letters. Vol. 32, p. 715, 1995. 

9. Vollmayr, K., W. Kob, and K Binder. Physical Review B. Vol. 54, p. 15808, 1996. 

10. Bemhek, S. D., and B. B. bird. Physical Review Letters. VoL 74, p. 936, 1995. 

11. Bedmek, S. D., and B. B. bird. Journal of Chemical Physics. Vol. 104, p. 1, 19$%. 

12. L&d., B. B., and S. D. Bembenek. Journal of Physics: Condensed Matter. VoL 8, 

p. 9569, 1996. 

13. Hansen, J., and L. Verlet. PhysicaZReview. Vol. 184, p. 151, 1969. 

14. Liu, K. S. Journal of Chemical Physics. Vol. 60, p. 4226, 1974. 

15. Iadd, A. J. C., and L. V. Woodcock. Chemical Physics Letters. VoL 5 1, p. 155, 1977. 

16. Toxvaerd, S., and E. F’raestgaard. Journal of Chemical Physics. VoL 67, p. 5291, 1977. 

17. Hiwatari, Y., E. Stall, and T. Schneider. Journal of Chemical Physics. VoL 68, p. 3401, 

1978. 

. 

18. Iadd, A. J. C.; and L. V. Woodcock. Journal of Physics C. Vol. 11, p. 3565,1978. 

19 
19. Bonissent, A., E. Gauthier, and J. L. Finney. Philosophical Magazine. VoL 39, p. 49,1979. 

20. Bushnell-Wye, G. M., J. L. Finney, and A. Bon&sent. Philosophical Magazine. VoL 44, 

p. 1053, 1981. 

21. Broughton, J. Q., and G. H. Gilmer. Journal of Chemicd Physics. VoL 84, p. 5741, 1986. 

22. Broughton, J. Q., and G. H. Gilmer. Journal of Chemical Physics. Vol. 84, p- 5749, 1986. 

23. Broughton, J. Q., and G. H. Gilmer. Journal of Chemical Physics. VoL 84, p. 5759, 1986. 

24, Nijmeije, M. J. P., A. F. Bakker, C. Bruin, and J. H. Sikkenk. Journal of Chemical Physics. 

VoL 89, pb 3789,1988. 

25. Holcomb, C. D., P. Clancy, S. M. Thompson, and J. A. Zollweg. FZuid Phase Equilibria. 

VoL 75, p. 185, 1992. 

26. Lot& A., J. Vrabec, and J. Fischer. Molecular Physics. Vol. 76, p. 1319, 1992. 

27. Barker, J. A., and D. Henderson. Reviews of Modern Physics. VoL 48, p. 587, 1976. 

28. Nicolas, J. J., K E. Gubbins, W. B. Street, and D. J. Tilde&y. MoZecuZar Physics. VoL 37, 

p* 1429, 1979. 

29. Johnson, K, J. A. Zollweg, and K E. Gubbins. Molecular Physics. VoL 78, p. 591,1993. 

30. AJlen, M. A., and D. J. Tildesley. Computer Simulation of Liquids. Oxford, England: 

Clarendon Press, 1987; Collaborative Computational Projects-5 (CCP5) Newsletter, no. 31, 

October 1989 (see http://~.dLac.uk/CCP/CCPS/newsletter-index.html). 

31. Frenkel, D., and B. Smit. Understanding Molecular Simulations. London, England: 

Academic Press, 1996. 

32. Evans, D. J., and G. P. Morriss. Computational Physics Reports. Vol. 1, p. 297, 1984. 

33. Andersen, H. C. Journal of Chemical Physics. VoL 72, p. 2384, 1980. 

34. Brown, D., and J. H. R. Clarke. Molecular Physics. Vol. 51, p. 1243, 1984. 

35. Bembenek, S. D. ‘Ph.D. Dissertation, University of Kansas, Lawrence, KS, August 1997. 

36. McDonald, I. R., and K. Singer. Journal of Chemical Physics. VoL 50, p. 2308, 1969. 

37. Wijker, H., A. Michels, and H. K. Wijker. Physica. Vol. 15, p. 627,1949. 

20 
38. Barker, J. A., R. A. Fisher, and R. 0. Watts. Molecular Physics. VoL 21, p. 657, 1971. 

39. Miyano, Y. Fluid Phase Equilibria. Vol. 95, p. 31, 1994. 

40. Sychev, V. V., A. A. Vassemm, A. D. Kozlov, G. A. Spiridonov, and V. A. Tsymamy. 

Themwdynamic Properties of Neon, Argon, Krypton, and Xenon. New York, NY: 

Hemisphere, 1987. 

41. Naghizadeh, J., and S. A. Rice. Journal of Chemical Physics. Vol. 36, p. 2710, 1962. 

42. Strmtt, W. B., and L. A. K Staveley. Journal of Chemical Physics. Vol. 50, p. 2302, 1969. 

21 
22 
NO. OF 

COPIES ORGANIZATION 

2 DEFENSE TECHNICAL 

INFORMATION CENTER 

DTIC DDA 

8725 JOHN J KINGMAN RD 

STE 0944 

FT BELVOJR VA 22060-62 1X 

1 HQDA 

DAM0 FDT 

400 ARMY PENTAGON 

WASHINGTON DC 203 1 O-0460 

1 OSD 

OUSD(A&T)/ODDDR&E(R) 

RJTREW 

THE PENTAGON 

WASHINGTON DC 20301-7100 

1 DPTY CG FOR RDA 

US ARMY MATERIEL CMD 

AMCRDA 

500 1 EISENHOWER AVE 

ALEXANDRLA VA 22333-0001 

1 INST FOR ADVNCD TCHNLGY 

THE UNIV OF TEXAS AT AUSTIN 

PO BOX 202797 

AUSTIN TX 78720-2797 

1 DARPA 

B KASPAR 

3701 N FAIRFAX DR 

ARLINGTON VA 22203-1714 

1 US MILITARY ACADEMY 

MATH SC1 CTR OF EXCELLENCE 

MADNMATH 

MAJ HUBER 

THAYER HALL 

WEST POINT NY 10996-1786 

1 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL D 

D R SMITH 

2800 POWDER MILL RD 

ADELPHI MD 20783-I 197 

23 

NO. OF 

COPIES ORGANIZATION 

1 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL DD 

2800 POWDER MILL RD 

ADELPHI MD 20783-1197 

1 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL CI AI R (RECORDS MGMT) 

2800 ,POWDER MILL RD 

ADELPHI MD 20783-l 145 

3 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL CI LL 

2800 POWDER MILL RD 

ADELPHI MD 20783-l 145 

1 DIRECTOR 

US ARMY RESEARCH LAB 

AMSRL CI AP 

2800 POWDER MILL RD 

ADELPHI MD 20783-I 197 

ABERDEEN PROVING GROUND 

4 DIRUSARL 

AMSRL CI LP (BLDG 305) 
NO. OF 

COPIES ORGANIZATION 

ABERDEEN PROVING GROUND 

21 DIR USARL 

AMSRL WM 

B RINGERS 

AMSRL WM BD 

B E FORCH 

W R ANDERSON 

S W BUNTE 

C F CHABALOWSKI 

A COHEN 

R DANIEL 

D DEVYNCK 

R A FIFER 

B E HOMAN 

A J KOTLAR 

K L MCNESBY 

M MCQUAID 

M S MILLER 

A W MIZIOLEK 

J B MORRIS 

R A PESCE-RODFkIGUEZ 

BMRICE 

R C SAUSA 

M A SCHROEDER 

J A VANDERHOFF 

24 
REPORT DOCUMENTATION PAGE *Approved 

OMB No. Oi’ObOl88 

wblit mponlng burden forms fOlktim of infomtlon 15 edmtedo avenge i now per msponxs. hcludhw me time foravlawhg hmclion5. 98ptching axkting dam SOYIOIJ. 

lanwing L”d maInPInIng the dam needd. and completing aid mvhwlng the coIIectlon or Information. smd conwmln$ regarding this tmrdan ~ahlut~ or any MnerafHct of this 

of inform~ttoe. including wgge%tlons for mduclng this burden, M WashIngton HaadqYallwS Servt.~s. Dlmctome for Info,mation OmBr;*lom and Rspons. 1215 Jotfels~l 

rransitioning Model Potentials to Real Systems 

RGANIZATION 

LJ.S Army Research Laboratory REPORT NUMBER 

ATTN: AMSRL-WM-BD 

Aberdeen Proving Ground, MD 2 10055066 ARL-l-R-2389 

11. SUPPLEMENTARY NOTES 

12a. DISTRIBUTIONIAVA~AB~L~N sTm3AENT 12b. DISTRIBUTION CODE 

Approved for public release; distribution is unlimited. 

13. ABSfRACT(Maximum 2W words) 

The parameters of two pairs of potentials that describe argon over its entire liquid phase at a fixed pressure wen 

optimized through a novel application of constant temperature and pressure molecular dpmics (NPT-MD) and Mont1 

Carlo (NPT-MC) computer simulations. The forms of these potentials were those of a modified Lennard-Jones potentia 

and a Lennard-Jones potential (Lennard-Jones, J. E. Physica, Vol. 4, p. 941, 1937). The optimized potential determine1 

using NPT-MD simulations reproduced experimental densities, internal energies, and enthalpies with an error less tbar 

1% over most of the liquid range and yielded self-diffusion coeffxcients that were in excellent agreement with thr 

Experiment. The results using the potential determined by NPT-MC simulations were in almost as good agreement witi 

deviations from experiments of no more than 5.89% for temperatures up to vaporization. Additionally, molar volume: 

predicted using this potential at pressures in the 100-600 abn range and over temperatures in the 1 O&140 K range werr 

within 0.83% of experimental values. These results showed that when properly parameterized, Lermard-Jones-lb 

potentials could describe a system well over a large temperature range. Further, the method inmoduced was easy tc 

implement and was independent of the form of the interaction potential used. 

14. SUBJECT TERMS 15. NUMBER OF PAGES 

molecular dynamics, Monte Carlo, oxygen, intermolecular potential 30 

16. PRICE CODE 

17. SECURtTY CLASSIFICATION 16. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LlMlTATlON OF ABSTRACT 

OF REPORT OF THIS PAGE OF ABSTRACT 

UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UL 

N$N754041-28Od5QO Standard Form 298 (Rev. Z-89) 25 Prescthd by ANSI Std. 239-i 8 298102 
. 

26 
USER EVALUATION SHEET/CHANGE OF ADDRESS 

This Laboratory undertakes a continuing effort to improve the quality of the reports it publishes. Your comments/answers to 

the items/questions below will aid us in our efforts. 

1. ARL Report Number/Author ARL-TR-2389 (Bembenek) Date of Report February 2001 

2. Date Report Received 

3. Does this report satisfy a need? 

. 

(Comment on purpose, related project, or other area of interest for which the report will be 

used.) 

4. Specifically, how is the report being used? (Information source, design data, procedure, source of ideas, etc.) 

5. Has the information in this report led to any quantitative savings as far as man-hours or dollars saved, operating costs 

avoided, or efficiencies achieved, etc? If so, please elaborate. 

6. General Comments. What do you think should be changed to improve future reports? (Indicate changes to organization, 

technical content, format, etc.) 

CURRENT 

ADDRESS 

Organization 

Name E-mail Name 

Street or P-0. Box No. 

City, State, Zip Code 

7. If indicating a Change of Address or Address Correction, please provide the Current or Correct address above and the Old or 

Incorrect address below. 

, Organization 

9 OLD Name 

ADDRESS 

Street or P.O. Box No. 

City, State, Zip Code 

(Remove this sheet, fold as indicated, tape closed, and mail.) 

(DO NOT STAPLE) 
DEPARTMENTOFTHE ARMY 

NECESSARY 

OFFICIAL BUSIHESS 

UNITED STATES 

POSTAGE WILL BE PAID BY ADDRESSEE 

I 

DIRECTOR 

US ARMY RESEARCH LABORATORY 

ATTN AMSRL WM BD 

ABERDEEN PROVING GROUND MD 21005-6066