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

https://www.nist.gov/mml/csd/chemical-informatics-group/spce-water-reference-calculations-non-cuboid-cell-10a-cutoff

[1]:

import unittest

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)

# FEASST version: v0.24.5-3-g7548e16503-hwh/srswmonoclinic4
MonteCarlo
Configuration particle_type0 /home/hwh/feasst/particle/spce.fstprt side_length0 30.0 side_length1 28.97777478867205 side_length2 29.51512917398008 xy 7.764571353075622 xyz_file /home/hwh/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:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 2.3037
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 2.3037
Potential Model ChargeScreened
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 2.3037
Potential Model ChargeScreenedIntra VisitModel VisitModelBond
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 2.3037
Potential Model ChargeSelf
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 2.3037
ThermoParams beta 1000000
Metropolis
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 2.3037
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 831
Log clear_file true max_precision true output_file en.csv
Run num_trials 1
# Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:574] 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 pandas as pd
class TestEwald0SPCE(unittest.TestCase):
    def test_srsw(self):
        df=pd.read_csv('en.csv')
        self.assertAlmostEqual(df['LennardJones'][0], 931.15451, delta=1e-4)
        self.assertAlmostEqual(df['LongRangeCorrections'][0], -34.16569, delta=1e-4)
        self.assertAlmostEqual(df['Ewald'][0], 371.46525, delta=1e-4)
        self.assertAlmostEqual(df['ChargeScreened'][0], -6046.43627, delta=1e-4)
        self.assertAlmostEqual(df['ChargeScreenedIntra'][0], 95078.89447, delta=1e-4)
        self.assertAlmostEqual(df['ChargeSelf'][0], -96297.75579, delta=1e-4)
unittest.main(argv=[''], verbosity=2, exit=False)

test_srsw (__main__.TestEwald0SPCE) ... ok

----------------------------------------------------------------------
Ran 1 test in 0.002s

OK

[2]:

<unittest.main.TestProgram at 0x7f184c312080>

[3]:

import unittest

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)

# FEASST version: v0.24.5-3-g7548e16503-hwh/srswmonoclinic4
MonteCarlo
Configuration particle_type0 /home/hwh/feasst/particle/spce.fstprt side_length0 36.0 side_length1 36.0 side_length2 31.17691453623979 xy 0.0 xyz_file /home/hwh/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:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 1.99017
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 1.99017
Potential Model ChargeScreened
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 1.99017
Potential Model ChargeScreenedIntra VisitModel VisitModelBond
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 1.99017
Potential Model ChargeSelf
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 1.99017
ThermoParams beta 1000000
Metropolis
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:179] alpha: 0.285
# Info 0 [plugin/charge/src/ewald.cpp:180] kmax_squared 1.99017
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
# Info 0 [plugin/charge/src/ewald.cpp:113] num vectors 1010
Log clear_file true max_precision true output_file en.csv
Run num_trials 1
# Warn 0 [plugin/monte_carlo/src/monte_carlo.cpp:574] No Trials to attempt.

 exit: 0

[7]:

import pandas as pd
class TestEwald0SPCE(unittest.TestCase):
    def test_srsw(self):
        df=pd.read_csv('en.csv')
        self.assertAlmostEqual(df['LennardJones'][0], 208.07025846554663, delta=1e-4)
        self.assertAlmostEqual(df['LongRangeCorrections'][0], -1.35601402285723, delta=1e-4)
        # Note that the Ewald value is slightly different than SRSW due to change in kmax
        self.assertAlmostEqual(df['Ewald'][0], 185.60955025935957, delta=1e-4)
        self.assertAlmostEqual(df['ChargeScreened'][0], -1425.6143894337906, delta=1e-2)
        self.assertAlmostEqual(df['ChargeScreenedIntra'][0], 23769.71831128121, delta=1e-2)
        self.assertAlmostEqual(df['ChargeSelf'][0], -24074.443366236523, delta=1e-2)
unittest.main(argv=[''], verbosity=2, exit=False)

test_srsw (__main__.TestEwald0SPCE) ... ok

----------------------------------------------------------------------
Ran 1 test in 0.002s

OK

[7]:

<unittest.main.TestProgram at 0x7f17d49b7250>

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