Lennard Jones potential test

In this example, an instantaneous configuration is compared against published values for the potential energy.

First, the potential energy between two particles is compared to the analytical value. Next, the energy of an entire configuration is compared against the NIST SRSW reference calculations.

import math
import sys
import subprocess
import numpy as np
import pandas as pd

# Run a feasst simulation with the given parameters
def run_fst(params):
    with open("launch0.txt", "w") as myfile: myfile.write("""
MonteCarlo
Configuration {config_params}
Potential Model LennardJones
Potential VisitModel LongRangeCorrections
ThermoParams beta 1000000
Metropolis
Log output_file lj.csv max_precision true clear_file true
Run num_trials 1
""".format(**params))
    syscode = subprocess.call("../../../build/bin/fst < launch0.txt", shell=True, executable='/bin/bash')
    if syscode > 0: sys.exit(1)

"""Test the LJ potential against analytical calculation of two particles"""
params = {"displacement": 1.2345}
with open("two.xyz", "w") as myfile: myfile.write(
"""2
-1 8 8 8
0 0 0 0
1 0 0 {displacement}""".format(**params))
run_fst({"config_params": "particle_type0 /feasst/particle/lj.fstprt xyz_file two.xyz"})
df = pd.read_csv('lj.csv')
assert 2 == df['num_particles_of_type0'][0]

# compute the expected analytical LJ and LRC energies
enlj = 4*(params["displacement"]**(-12) - params["displacement"]**(-6))
rcut = 3 #mc.system().configuration().model_params().select("cutoff").value(0)
enlrc = (8./3.)*math.pi*2**2/ 8**3*((1./3.)*rcut**(-9) - rcut**(-3))

# Compare the analytical results with the FEASST computed energies.
assert np.abs(enlj - df['LennardJones'][0]) < 12
assert np.abs(enlrc - df['LongRangeCorrections'][0]) < 12
assert np.abs(enlj + enlrc - df['energy'][0]) < 12

"""Test the LJ potential against a configuration of 30 particles.
In particular, the 4th configuration of the LJ SRSW reference:
https://www.nist.gov/mml/csd/chemical-informatics-research-group/lennard-jones-fluid-reference-calculations
"""
run_fst({"config_params": "cubic_side_length 8 particle_type0 /feasst/particle/lj.fstprt \
          xyz_file /feasst/plugin/configuration/test/data/lj_sample_config_periodic4.xyz"})
df = pd.read_csv('lj.csv')
assert df['num_particles_of_type0'][0] == 30
enlj = -16.790321304625856
enlrc = -0.5451660014945704
assert np.abs(enlj - df['LennardJones'][0]) < 12
assert np.abs(enlrc - df['LongRangeCorrections'][0]) < 12
assert np.abs(enlj + enlrc - df['energy'][0]) < 12

"""Test the LJ potential against a configuration of 300 particles in a trinclinic cell.
In particular, the 3th configuration of the triclinic LJ SRSW reference:
https://www.nist.gov/mml/csd/chemical-informatics-group/lennard-jones-fluid-reference-calculations-non-cuboid-cell
"""
run_fst({"config_params": "side_length0 10.0 side_length1 9.84807753012208 side_length2 9.64974312607518 \
    xy 1.7364817766693041 xz 2.5881904510252074 yz 0.42863479791864567 \
    particle_type0 /feasst/particle/lj.fstprt \
    xyz_file /feasst/plugin/configuration/test/data/lj_triclinic_sample_config_periodic3.xyz"})
df = pd.read_csv('lj.csv')
assert df['num_particles_of_type0'][0] == 300
enlj = -505.78567945268367
enlrc = -29.37186430697248
assert np.abs(enlj - df['LennardJones'][0]) < 12
assert np.abs(enlrc - df['LongRangeCorrections'][0]) < 12
assert np.abs(enlj + enlrc - df['energy'][0]) < 12

# Usage: ./fst < file.txt
FEASST version 0.25.6
MonteCarlo
Configuration particle_type0 /home/user/feasst/particle/lj.fstprt xyz_file two.xyz
Potential Model LennardJones
Potential VisitModel LongRangeCorrections
ThermoParams beta 1000000
Metropolis
Log clear_file true max_precision true output_file lj.csv
Run num_trials 1
#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:744] No Trials to attempt.
# Usage: ./fst < file.txt
FEASST version 0.25.6
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/user/feasst/particle/lj.fstprt xyz_file /home/user/feasst/plugin/configuration/test/data/lj_sample_config_periodic4.xyz
Potential Model LennardJones
Potential VisitModel LongRangeCorrections
ThermoParams beta 1000000
Metropolis
Log clear_file true max_precision true output_file lj.csv
Run num_trials 1
#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:744] No Trials to attempt.
# Usage: ./fst < file.txt
FEASST version 0.25.6
MonteCarlo
Configuration particle_type0 /home/user/feasst/particle/lj.fstprt side_length0 10.0 side_length1 9.84807753012208 side_length2 9.64974312607518 xy 1.7364817766693041 xyz_file /home/user/feasst/plugin/configuration/test/data/lj_triclinic_sample_config_periodic3.xyz xz 2.5881904510252074 yz 0.42863479791864567
Potential Model LennardJones
Potential VisitModel LongRangeCorrections
ThermoParams beta 1000000
Metropolis
Log clear_file true max_precision true output_file lj.csv
Run num_trials 1
#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:744] No Trials to attempt.

Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!