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.

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
            )

%%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

<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!