Grand canonical ensemble transition-matrix Monte Carlo with SPC/EΒΆ
This example is similar to the Lennard Jones tutorial, except this time we simulate the SPC/E water model.
[1]:
import unittest
import feasst as fst
class TestFlatHistogramSPCE(unittest.TestCase):
"""Test flat histogram grand canonical ensemble Monte Carlo simulations"""
def test_serial_5max(self):
"""Compare the free energies and potential energies with the NIST SRSW
https://www.nist.gov/programs-projects/nist-standard-reference-simulation-website
https://mmlapps.nist.gov/srs/LJ_PURE/eostmmc.htm
"""
monte_carlo = fst.MonteCarlo()
monte_carlo.set(fst.spce(fst.args({
"cubic_box_length": "20",
"physical_constants": "CODATA2010",
"alphaL": "5.6",
"kmax_squared": "38"})))
temperature = fst.kelvin2kJpermol(525, monte_carlo.configuration())
monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1/temperature),
"chemical_potential": str(-8.14*temperature)})))
criteria = fst.MakeFlatHistogram(
fst.MakeMacrostateNumParticles(fst.Histogram(fst.args({"width": "1", "max": "5", "min": "0"}))),
fst.MakeTransitionMatrix(fst.args({"min_sweeps": "20"})))
monte_carlo.set(criteria)
steps_per = str(1e4)
monte_carlo.add(fst.MakeTrialTranslate())
monte_carlo.add(fst.MakeTrialTransfer(fst.args({"particle_type": "0"})))
monte_carlo.add(fst.MakeCriteriaUpdater(fst.args({"steps_per": steps_per})))
monte_carlo.add(fst.MakeCriteriaWriter(fst.args({"steps_per": steps_per,
"file_name": "spce_fh.txt"})))
monte_carlo.add(fst.MakeEnergy(fst.args({"steps_per_write": steps_per,
"file_name": "spce_en.txt",
"multistate": "True"})))
monte_carlo.run_until_complete()
lnpi_previous = [
[-2.7207, 0.015],
[-1.8523, 0.015],
[-1.54708, 0.016],
[-1.51786, 0.015],
[-1.6479, 0.015],
[-1.8786, 0.03]]
energy_previous = [
[0, 1e-13],
[-0.0879115, 1.1293158298007674394e-05],
[-2.326, 0.12],
[-6.806, 0.24],
[-13.499, 0.5],
[-22.27, 1.0]]
for macro in range(criteria.num_states()):
self.assertAlmostEqual(
lnpi_previous[macro][0],
criteria.bias().ln_prob().value(macro),
delta=15.*lnpi_previous[macro][1]
)
energy_analyzer = monte_carlo.analyze(monte_carlo.num_analyzers() - 1)
energy_accumulator = energy_analyzer.analyze(macro).accumulator()
stdev = (energy_previous[macro][1]**2 + energy_accumulator.block_stdev()**2)**(1./2.)
self.assertAlmostEqual(
energy_previous[macro][0],
energy_accumulator.average(),
#criteria.bias().ln_prob().value(macro),
delta=15*stdev
)
[2]:
%%time
unittest.main(argv=[''], verbosity=2, exit=False)
test_serial_5max (__main__.TestFlatHistogramSPCE)
Compare the free energies and potential energies with the NIST SRSW ...
CPU times: user 8.15 s, sys: 28.9 ms, total: 8.17 s
Wall time: 8.14 s
ok
----------------------------------------------------------------------
Ran 1 test in 8.130s
OK
[2]:
<unittest.main.TestProgram at 0x7fc355534940>
Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!