Lennard Jones Alpha potential test

In this example, the potential energy between two particles is compared to the analytical value.

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 LennardJonesAlpha alpha 12
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/plugin/models/particle/ljdelta.fstprt xyz_file two.xyz delta_sigma 0.5"})
df = pd.read_csv('lj.csv')
assert 2 == df['num_particles_of_type0'][0]

# compute the expected analytical LJ energy
enlj = 4*((params["displacement"]+0.5)**(-24) - (params["displacement"]+0.5)**(-12))

# Compare the analytical results with the FEASST computed energies.
assert np.abs(enlj - df['LennardJonesAlpha'][0]) < 10**-12
assert np.abs(enlj - df['energy'][0]) < 10**-12

# Usage: ./fst < file.txt
FEASST version 0.25.7
MonteCarlo
Configuration delta_sigma 0.5 particle_type0 /home/user/feasst/plugin/models/particle/ljdelta.fstprt xyz_file two.xyz
Potential Model LennardJonesAlpha alpha 12
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:743] No Trials to attempt.

"""Similar to above, but with the lambda parameter"""
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/plugin/models/particle/ljlambda.fstprt xyz_file two.xyz"})
df = pd.read_csv('lj.csv')
assert 2 == df['num_particles_of_type0'][0]

# compute the expected analytical LJ energy
enlj = 2*((params["displacement"]+0.5)**(-24) - (params["displacement"]+0.5)**(-12))

# Compare the analytical results with the FEASST computed energies.
assert np.abs(enlj - df['LennardJonesAlpha'][0]) < 10**-12
assert np.abs(enlj - df['energy'][0]) < 10**-12

# Usage: ./fst < file.txt
FEASST version 0.25.7
MonteCarlo
Configuration particle_type0 /home/user/feasst/plugin/models/particle/ljlambda.fstprt xyz_file two.xyz
Potential Model LennardJonesAlpha alpha 12
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:743] 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!