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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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!