Lennard Jones in the canonical ensemble

In this example, the average energy reported in the NIST SRSW at \(T^*=0.9\) and \(\rho^*=0.003\) is computed with Monte Carlo simulation of a bulk Lennard Jones fluid (Section IIIB of the manuscript). The simulation initializes a desired number of particles in cubic periodic boundaries, tunes the maximum trial displacement during equilibration, and finally obtains the ensemble average energy.

To begin, Python3 generates a text file, fstin.txt, which input to the fst executable using the command fst < fstin.txt. FEASST does not require the use of Python and users may create the FEASST text input file in other ways. We use Python to generate a text file for the following conveniences:

• Variables required by FEASST arguments may need to be converted from variables that the user prefers to input, and these FEASST arguments can 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 simulates a range of densities with a single Python file.

• Python provides formatted strings and argument parsing.

In this example, we begin a simulation by defining a string variable fstin. Optional comments can be included for future reference. First, a MonteCarlo object is created. Then, the 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 input as pairs of argument names then values separated by an equal sign (e.g., “name=value”). Each argument pair is separated by a space. In this case, seed is the argument name and 1572362164 is the argument value.

# This and all tutorial code blocks are in Python and included in your download of FEASST.
# For example, https://pages.nist.gov/feasst/tutorial/tutorial.html is in $HOME/feasst/tutorial/tutorial.ipynb . Usage: jupyter notebook tutorial.ipynb
# And https://pages.nist.gov/feasst/tutorial/launch.html is in $HOME/feasst/tutorial/launch.py . Usage: python launch.py

import subprocess

fstin="""
# comments begin with the '#' symbol and are not displayed in the output
MonteCarlo
RandomMT19937 seed=1572362164
"""

def run(fstin):
    with open('fstin.txt', 'w') as file: file.write(fstin)
    syscode = subprocess.call("../build/bin/fst < fstin.txt > fstout.txt", 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('fstout.txt', 'r') as file:
            print(file.read())
run(fstin)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=1572362164
# Initializing random number generator with seed: 1572362164

The output of fst in fstout.txt 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, a particle type “lj” is defined by the file lj.txt. Any argument value beginning with “/feasst” will have that beginning replaced with the installation directory of FEASST). See Particle for more information about the format of the particle file.

num_particles=500
density=0.003
fstin+="""
Configuration cubic_side_length={length} particle_type=lj:/feasst/particle/lj.txt
""".format(length=(num_particles/density)**(1./3.))
run(fstin)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=1572362164
# Initializing random number generator with seed: 1572362164
Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj.txt

Two Potentials are introduced starting with the pair-wise Lennard-Jones (LJ) model. Because the domain is large compared to the potential cutoff, a cell list speeds up the energy calculation. Then, long-range corrections approximately account for the cut off of the LJ potential by assuming a pair-wise radial distribution function of unity.

fstin+="""
Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff
Potential VisitModel=LongRangeCorrections
"""

ThermoParams set temperature and the chemical potential of each particle type. A chemical potential is required when the initial configuration is generated with grand canonical particle additions, but the chemical potential will not contribute to canonical ensemble production simulations. Metropolis acceptance criteria is also utilized.

fstin+="""
ThermoParams beta={beta} chemical_potential=-1
Metropolis
""".format(beta=1./0.9)
run(fstin)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=1572362164
# Initializing random number generator with seed: 1572362164
Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj.txt
Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff
Potential VisitModel=LongRangeCorrections
ThermoParams beta=1.1111111111111112 chemical_potential=-1
Metropolis

A TrialTranslate attempts translation of a random particle by a random distance bound in each dimension by a tunable_param. The desired acceptance ratio, tunable_target_acceptance, is obtained by adjusting this parameter with Tune.

fstin+="""
TrialTranslate weight=1 tunable_param=2 tunable_target_acceptance=0.2
Tune
"""

CheckEnergy asserts that the optimized energy calculations match the unoptimized calculations within a tolerance. This check is only performed every trials_per_update because the unoptimized energy calculation is expensive.

fstin+="""CheckEnergy trials_per_update=1e4 decimal_places=8"""

TrialAdd, generates an initial configuration with a Run with the desired number of particles. Afterward, we Remove TrialAdd to simulate the canonical ensemble.

fstin+="""
TrialAdd weight=2 particle_type=lj
Run until_num_particles={num_particles}
Remove name=TrialAdd
""".format(num_particles=num_particles)
run(fstin)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=1572362164
# Initializing random number generator with seed: 1572362164
Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj.txt
Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff
Potential VisitModel=LongRangeCorrections
ThermoParams beta=1.1111111111111112 chemical_potential=-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=lj weight=2
Run until_num_particles=500
Remove name=TrialAdd

The simulation Run for a number of equilibration trials before Remove disables tuning of the maximum displacement (to satisfy detailed balance).

fstin+="""
Run num_trials=1e5
Remove name=Tune
"""

Finally, initialize a production simulation. First, a Log periodically outputs the instantaneous properties of the simulation. Log is useful for monitoring simulation progress, but is not recommended for calculations. Instead, Energy computes the average energy by accumulating its value after every trial (more often than once every line in the Log file). Movie produces a trajectory in the XYZ file format. Additional Analyze or Modify may be added to perform some task contingent upon the number of attempted trials.

Because all three of these Analyze share similar arguments, Let minimizes redundancy by defining a string of characters as the variable [write]. Variables must begin and end with square brackets followed by an equal sign to set its value.

fstin+="""
Let [write]=trials_per_write=1e4 output_file=lj
Log [write].csv
Energy [write]_en.csv
Movie [write].xyz
"""

The next code block is equivalent to the previous, but demonstrates the use of For loops. Each For executes the following lines with given string substitutions until EndFor, and then repeats those lines again for each comma-separated list of arguments. Within the same For loop, multiple string substitutions are defined with colon-separated values. For example, [Analyze]=Log and [extension]=.csv in the first iteraction of the loop, as can be seen in the output. If, Else and EndIf may also be used for optional arguments.

fstin+="""
For [an]:[ext]=Log:.csv,Energy:_en.csv,Movie:.xyz
    [an] trials_per_write=1e4 output_file=lj_with_for[ext]
EndFor
Run num_trials=1e6
"""
run(fstin)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
RandomMT19937 seed=1572362164
# Initializing random number generator with seed: 1572362164
Configuration cubic_side_length=55.03212081491043 particle_type=lj:/feasst/particle/lj.txt
Potential Model=LennardJones VisitModel=VisitModelCell min_length=max_cutoff
Potential VisitModel=LongRangeCorrections
ThermoParams beta=1.1111111111111112 chemical_potential=-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=lj weight=2
Run until_num_particles=500
Remove name=TrialAdd
Run num_trials=1e5
Remove name=Tune
Log output_file=lj.csv trials_per_write=1e4
Energy output_file=lj_en.csv trials_per_write=1e4
Movie output_file=lj.xyz trials_per_write=1e4
Log output_file=lj_with_for.csv trials_per_write=1e4
Energy output_file=lj_with_for_en.csv trials_per_write=1e4
Movie output_file=lj_with_for.xyz trials_per_write=1e4
Run num_trials=1e6

Compare the average energy to the NIST SRSW

import pandas as pd
df=pd.read_csv('lj_en.csv')
print('U_FEASST/N=', df['average'][0]/num_particles, '+/-', df['block_stdev'][0]/num_particles)
print('U_SRSW/N=', -2.9787E-02, '+/-', 3.21E-05)
import math
assert abs(-2.9787E-02-df['average'][0]/num_particles) < 3*math.sqrt(3.21E-05**2+(df['block_stdev'][0]/num_particles)**2)

U_FEASST/N= -0.030332409560752716 +/- 0.0002319803861581972
U_SRSW/N= -0.029787 +/- 3.21e-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!