Reference configuration of SPC/E water in non-cuboid domain

Compare the potential energy against the NIST SRSW values. This tutorial requires the use of Ewald.

[1]:

script="""
MonteCarlo
Configuration side_length=30.0,28.97777478867205,29.51512917398008 \
    xy=7.764571353075622 xz=-2.6146722824297473 yz=-4.692615336756641 \
    xyz_file=/feasst/plugin/charge/test/data/spce_triclinic_sample_periodic1.xyz \
    particle_type=spce:/feasst/particle/spce_new.txt
Potential Model=LennardJones
Potential VisitModel=LongRangeCorrections
Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7
Potential Model=ChargeScreened
Potential Model=ChargeScreenedIntra VisitModel=VisitModelBond
Potential Model=ChargeSelf
ThermoParams beta=1000000
Metropolis
Log output_file=en.csv max_precision=true clear_file=true
Run num_trials=1
"""

with open('script.txt', 'w') as file: file.write(script)
import subprocess
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)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration particle_type=spce:/feasst/particle/spce_new.txt side_length=30.0,28.97777478867205,29.51512917398008 xy=7.764571353075622 xyz_file=/home/user/feasst/plugin/charge/test/data/spce_triclinic_sample_periodic1.xyz xz=-2.6146722824297473 yz=-4.692615336756641
Potential Model=LennardJones
Potential VisitModel=LongRangeCorrections
Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7
# Ewald alpha: 0.285
# Ewald kmax_squared 2.3037
# Ewald alpha: 0.285
# Ewald kmax_squared 2.3037
Potential Model=ChargeScreened
# Ewald alpha: 0.285
# Ewald kmax_squared 2.3037
Potential Model=ChargeScreenedIntra VisitModel=VisitModelBond
# Ewald alpha: 0.285
# Ewald kmax_squared 2.3037
Potential Model=ChargeSelf
# Ewald alpha: 0.285
# Ewald kmax_squared 2.3037
ThermoParams beta=1000000
Metropolis
# Ewald alpha: 0.285
# Ewald kmax_squared 2.3037
Log clear_file=true max_precision=true output_file=en.csv
Run num_trials=1
#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:790] No Trials to attempt.

 exit: 0

If the test passes, the energy is within the tolerance of the SRSW value and the two ensemble average methods agreed.

[2]:

import numpy as np
import pandas as pd
df=pd.read_csv('en.csv')
assert np.abs(df['LennardJones'][0] - 931.15451) < 1e-4
assert np.abs(df['LongRangeCorrections'][0] + 34.16569) < 1e-4
assert np.abs(df['Ewald'][0] - 371.46525) < 1e-4
assert np.abs(df['ChargeScreened'][0] + 6046.43627) < 1e-4
assert np.abs(df['ChargeScreenedIntra'][0] - 95078.89447) < 1e-4
assert np.abs(df['ChargeSelf'][0] + 96297.75579) < 1e-4

[3]:

script="""
MonteCarlo
Configuration side_length=36.0,36.0,31.17691453623979 \
    xy=0.0 xz=18.0 yz=0.0 \
    xyz_file=/feasst/plugin/charge/test/data/spce_monoclinic_sample_periodic4.xyz \
    particle_type=spce:/feasst/particle/spce_new.txt
Potential Model=LennardJones
Potential VisitModel=LongRangeCorrections
Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7
Potential Model=ChargeScreened
Potential Model=ChargeScreenedIntra VisitModel VisitModelBond
Potential Model=ChargeSelf
ThermoParams beta=1000000
Metropolis
Log output_file=en.csv max_precision=true clear_file=true
Run num_trials=1
"""

with open('script.txt', 'w') as file: file.write(script)
import subprocess
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)

# Usage: /home/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration particle_type=spce:/feasst/particle/spce_new.txt side_length=36.0,36.0,31.17691453623979 xy=0.0 xyz_file=/home/user/feasst/plugin/charge/test/data/spce_monoclinic_sample_periodic4.xyz xz=18.0 yz=0.0
Potential Model=LennardJones
Potential VisitModel=LongRangeCorrections
Potential VisitModel=Ewald alpha=0.2850 kxmax=7 kymax=7 kzmax=7
# Ewald alpha: 0.285
# Ewald kmax_squared 1.99017
# Ewald alpha: 0.285
# Ewald kmax_squared 1.99017
Potential Model=ChargeScreened
# Ewald alpha: 0.285
# Ewald kmax_squared 1.99017
Potential Model=ChargeScreenedIntra VisitModel=VisitModelBond
# Ewald alpha: 0.285
# Ewald kmax_squared 1.99017
Potential Model=ChargeSelf
# Ewald alpha: 0.285
# Ewald kmax_squared 1.99017
ThermoParams beta=1000000
Metropolis
# Ewald alpha: 0.285
# Ewald kmax_squared 1.99017
Log clear_file=true max_precision=true output_file=en.csv
Run num_trials=1
#Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:790] No Trials to attempt.

 exit: 0

[4]:

df=pd.read_csv('en.csv')
assert np.abs(df['LennardJones'][0] - 208.07025846554663) < 1e-4
assert np.abs(df['LongRangeCorrections'][0] + 1.35601402285723) < 1e-4
# Note that the Ewald value is slightly different than SRSW due to change in kmax
assert np.abs(df['Ewald'][0] - 185.60955025935957) < 1e-4
assert np.abs(df['ChargeScreened'][0] + 1425.6143894337906) < 1e-2
assert np.abs(df['ChargeScreenedIntra'][0] - 23769.71831128121) < 1e-2
assert np.abs(df['ChargeSelf'][0] + 24074.443366236523) < 1e-2

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