Lennard Jones in the canonical ensemble
In this example, we will reproduce the average energy reported in the NIST SRSW at T*=0.9 and \(\rho*=0.009\) using a Monte Carlo simulation of a bulk Lennard Jones fluid. The simulation will be initialized with a desired number of particles in a cubic box, equilibrated while tuning the maximum trial displacement before obtaining the ensemble average energy.
To begin, a Python string, script
, is used to write the file, script.txt
, which is then run by the fst
executable. Within script
, a MonteCarlo object is created. Then, a random number generator is given a seed. The available arguments of RandomMT19937 can be found in the documentation of that object, or its base-class, Random. Arguments of objects are separated by
a space and are provided first as the name of the argument and then the value. In this case, seed
is the argument name and 1572362164
is the argument value.
[1]:
import subprocess
script="""
# comments begin with the '#' symbol
MonteCarlo
RandomMT19937 seed 1572362164
"""
def run(script):
with open('script.txt', 'w') as file: file.write(script)
syscode = subprocess.call("../build/bin/fst < script.txt > script.log", shell=True, executable='/bin/bash')
with open('script.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)
run(script)
# FEASST version: v0.22.2-11-g5377b2e66a-dirty-hwh/txttrlwhite
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
exit: 0
The output of fst
in script.log
is displayed. It includes the FEASST version as well as a confirmation of the seed. Note that the version shown may not match your version.
If there is a typo in the arguments, this may result in an exception which will print to the terminal. You may test this now by changing the seed
argument to sd
. Because the argument sd
is not recognized by RandomMT19937, a nonzero system code should be returned and a verbose error should appear in your terminal.
The next step is to add a Configuration. In this example, a simple cubic periodic box is defined based on the number of particles and density. Because the cubic length needs to be input to high precision, a Python format string is used to input a variable enclosed in curly brackets.
In addition, the first particle type with an index of 0
is defined by the file lj.fstprt
. Any argument value beginning with “/feasst” will have that beginning replaced with feasst::install_dir(), and thus the full path to the file lj.fstprt
is displayed in script.log
shown below. See Particle for more information about the format of the data file, which is a LAMMPS-inspired file with some major differences.
[2]:
num_particles=500
density=0.009
script+="""
Configuration cubic_side_length {length} particle_type0 /feasst/particle/lj.fstprt
""".format(length=(num_particles/density)**(1./3.))
run(script)
# FEASST version: v0.22.2-11-g5377b2e66a-dirty-hwh/txttrlwhite
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/hwh/feasst/particle/lj.fstprt
exit: 0
Two Potentials are introduced. First, the pair-wise Lennard-Jones (LJ) model. Because the domain is large compared to the potential cutoff, a cell list is used. Second, long-range corrections, which approximately account for the cut off of the LJ potential by assuming a pair-wise radial distribution function of unity.
[3]:
script+="""
Potential Model LennardJones VisitModel VisitModelCell min_length max_cutoff
Potential VisitModel LongRangeCorrections
"""
Next, set ThermoParams, such as temperature and chemical potential of each particle type. The initial configuration will be generated with grand canonical particle additions, so a chemical potential is included, but will not contribute to canonical ensemble production simulations.
In addition, Metropolis acceptance criteria is utilized.
[4]:
script+="""
ThermoParams beta {beta} chemical_potential0 -1
Metropolis
""".format(beta=1./0.9)
A TrialTranslate is then introduced which attempts to translate a random particle by a random distance which is bound in each dimension by a tunable_param
. This parameter may be adjusted to obtain a desired acceptance ratio, tunable_target_acceptance
, with the help of Tune.
[5]:
script+="""
TrialTranslate weight 1 tunable_param 2 tunable_target_acceptance 0.2
Tune
"""
run(script)
# FEASST version: v0.22.2-11-g5377b2e66a-dirty-hwh/txttrlwhite
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/hwh/feasst/particle/lj.fstprt
Potential Model LennardJones VisitModel VisitModelCell min_length max_cutoff
Potential VisitModel LongRangeCorrections
ThermoParams beta 1.1111111111111112 chemical_potential0 -1
Metropolis
TrialTranslate tunable_param 2 tunable_target_acceptance 0.2 weight 1
Tune
exit: 0
CheckEnergy asserts that the optimized energy calculations match the unoptimized calculations within a tolerance. Because the unoptimized energy calculation can be expensive, this check is only performed every trials_per
. The line set_variable name value
will replace any use of name
in subsequent values, as shown in script.log
.
[6]:
script+="""
set_variable trials_per 1e4
CheckEnergy trials_per_update trials_per tolerance 1e-8
"""
With the help of TrialAdd, we can now generate an initial configuration with a Run until the desired number of particles are reached. RemoveTrial is then used to revert back to the canonical ensemble.
[7]:
script+="""
TrialAdd weight 2 particle_type 0
Run until_num_particles {num_particles}
RemoveTrial name TrialAdd
""".format(num_particles=num_particles)
run(script)
# FEASST version: v0.22.2-11-g5377b2e66a-dirty-hwh/txttrlwhite
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/hwh/feasst/particle/lj.fstprt
Potential Model LennardJones VisitModel VisitModelCell min_length max_cutoff
Potential VisitModel LongRangeCorrections
ThermoParams beta 1.1111111111111112 chemical_potential0 -1
Metropolis
TrialTranslate tunable_param 2 tunable_target_acceptance 0.2 weight 1
Tune
CheckEnergy tolerance 1e-8 trials_per_update 1e4
TrialAdd particle_type 0 weight 2
Run until_num_particles 500
RemoveTrial name TrialAdd
exit: 0
The simulation is then Run for a number of equilibration trials before RemoveModify disables tuning of the maximum displacement.
[8]:
script+="""
Run num_trials 1e5
RemoveModify name Tune
"""
Finally, it is time to obtain the ensemble average energy from a production simulation. First, a Log periodically outputs the instantaneous status of the trials, but it is more useful for watching the progress of a simulation than actually calculating quantities. Instead of analyzing the Log file, Energy may be used to more precisely compute the average energy by accumulating its value after every trial. A trajectory in the XYZ format may also be produced using Movie.
Additional Analyze or Modify may be added to perform some task contingent upon the number of attempted trials.
[9]:
script+="""
Log trials_per_write trials_per file_name lj.txt
Energy trials_per_write trials_per file_name lj_en.txt
Movie trials_per_write trials_per file_name lj.xyz
Run num_trials 1e6
"""
run(script)
# FEASST version: v0.22.2-11-g5377b2e66a-dirty-hwh/txttrlwhite
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/hwh/feasst/particle/lj.fstprt
Potential Model LennardJones VisitModel VisitModelCell min_length max_cutoff
Potential VisitModel LongRangeCorrections
ThermoParams beta 1.1111111111111112 chemical_potential0 -1
Metropolis
TrialTranslate tunable_param 2 tunable_target_acceptance 0.2 weight 1
Tune
CheckEnergy tolerance 1e-8 trials_per_update 1e4
TrialAdd particle_type 0 weight 2
Run until_num_particles 500
RemoveTrial name TrialAdd
Run num_trials 1e5
RemoveModify name Tune
Log file_name lj.txt trials_per_write 1e4
Energy file_name lj_en.txt trials_per_write 1e4
Movie file_name lj.xyz trials_per_write 1e4
Run num_trials 1e6
exit: 0
Compare the average energy to the NIST SRSW
[10]:
import pandas as pd
df=pd.read_csv('lj_en.txt')
print('U_FEASST/N=', df['average'][0]/num_particles, '+/-', df['block_stdev'][0]/num_particles)
print('U_SRSW/N=', -8.9936E-02, '+/-', 2.44E-05)
import math
assert abs(-8.9936E-02-df['average'][0]/num_particles) < 3*math.sqrt(2.44E-05**2+(df['block_stdev'][0]/num_particles)**2)
U_FEASST/N= -0.09030193615567321 +/- 0.00048147353358693576
U_SRSW/N= -0.089936 +/- 2.44e-05
You should also find the lj.xyz
trajectory file with an automatically-generated lj.xyz.vmd
file for use with VMD (e.g., vmd -e lj.xyz.vmd
).
Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!