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_length0 30.0 side_length1 28.97777478867205 side_length2 29.51512917398008 \
    xy 7.764571353075622 xz -2.6146722824297473 yz -4.692615336756641 \
    xyz_file /feasst/plugin/charge/test/data/spce_triclinic_sample_periodic1.xyz \
    particle_type0 /feasst/particle/spce.fstprt
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: ./fst < file.txt
FEASST version 0.25.6
MonteCarlo
Configuration particle_type0 /home/user/feasst/particle/spce.fstprt side_length0 30.0 side_length1 28.97777478867205 side_length2 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
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 2.3037
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 2.3037
Potential Model ChargeScreened
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 2.3037
Potential Model ChargeScreenedIntra VisitModel VisitModelBond
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 2.3037
Potential Model ChargeSelf
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 2.3037
ThermoParams beta 1000000
Metropolis
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] 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:744] 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_length0 36.0 side_length1 36.0 side_length2 31.17691453623979 \
    xy 0.0 xz 18.0 yz 0.0 \
    xyz_file /feasst/plugin/charge/test/data/spce_monoclinic_sample_periodic4.xyz \
    particle_type0 /feasst/particle/spce.fstprt
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: ./fst < file.txt
FEASST version 0.25.6
MonteCarlo
Configuration particle_type0 /home/user/feasst/particle/spce.fstprt side_length0 36.0 side_length1 36.0 side_length2 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
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 1.99017
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 1.99017
Potential Model ChargeScreened
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 1.99017
Potential Model ChargeScreenedIntra VisitModel VisitModelBond
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 1.99017
Potential Model ChargeSelf
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] kmax_squared 1.99017
ThermoParams beta 1000000
Metropolis
#Info 0 [plugin/charge/src/ewald.cpp:186] alpha: 0.285
#Info 0 [plugin/charge/src/ewald.cpp:187] 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:744] 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!