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, Python is used to generate a text file, script.txt, which is then input to the fst executable using the command fst < script.txt. We recommend using Python to generate a text file, instead of manually editting a text file, for the following reasons:

• Variables required by FEASST arguments may need to be converted from variables that the user prefers to input, and these FEASST arguments can then be generated by Python to machine precision. For example, if the user prefers to input temperature and state a temperature of \(T^*=0.9\) in their manuscript, but FEASST requires \(\beta=\frac{1}{k_B T}\), then \(T^*=0.9\) is easier to input than \(\beta=1.111111111111111\).

• Python generation of input files may create multiple related simulations from a single source. For example, the next tutorial shows how to simulate a range of densities.

• Python provides formatted strings and argument parsing.

In this example, we begin a simulation by defining a Python string script, where comments can be input for future reference, and 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.

# This and all tutorial code blocks are in Python.
# Copy/paste or use the URL to find the code.
# For example, https://pages.nist.gov/feasst/tutorial/tutorial.html is $HOME/feasst/tutorial/tutorial.ipynb . Usage: jupyter notebook tutorial.ipynb
# And https://pages.nist.gov/feasst/tutorial/launch.html is $HOME/feasst/tutorial/launch.py . Usage: python launch.py

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')
    if syscode != 0:
        print('syscode', syscode)
        print('If you see "/bin/bash: line 1: ../build/bin/fst: No such file or directory" then either fst was not compiled or you need to change the path to fst in subprocess.call(...) two lines above.')
    else:
        with open('script.log', 'r') as file:
            print(file.read())
run(script)

# FEASST version: v0.25.1
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164

The output of fst in script.log is displayed, and may be used to reproduce the simulation. It includes the FEASST version as well as a confirmation of the seed. 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.23.1-18-gb1bcd7c84e-dirty-user/234
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/user/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, are also initialized.

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.23.1-18-gb1bcd7c84e-dirty-user/234
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/user/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 decimal_places 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.23.1-18-gb1bcd7c84e-dirty-user/234
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/user/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 decimal_places 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 output_file lj.txt
Energy trials_per_write trials_per output_file lj_en.txt
Movie trials_per_write trials_per output_file lj.xyz
Run num_trials 1e6
"""
run(script)

# FEASST version: v0.23.1-18-gb1bcd7c84e-dirty-user/234
MonteCarlo
RandomMT19937 seed 1572362164
# initializing random number generator with seed: 1572362164
Configuration cubic_side_length 38.15714141844439 particle_type0 /home/user/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 decimal_places 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 output_file lj.txt trials_per_write 1e4
Energy output_file lj_en.txt trials_per_write 1e4
Movie output_file 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!