Canonical ensemble Monte Carlo of bulk SPC/E waterΒΆ

In this example, we reproduce the average energy reported in https://doi.org/10.1063/1.476834

[1]:
import unittest
import feasst as fst

class TestEwald1SPCENVT(unittest.TestCase):
    """Test a canonical ensemble Lennard Jones Monte Carlo simulation"""
    def test_srsw(self, num_particles=512, box_length=24.8586887, steps_per=1e5, temperature=298):
        """Compare with the reported average energy from https://doi.org/10.1063/1.476834.

        num_particles -- number of water molecules
        box_length -- box length in angstroms
        steps_per -- steps between each Anaylze/Modify
        temperature -- in Kelvin
        """

        monte_carlo = fst.MonteCarlo()
        # monte_carlo.set(fst.MakeRandomMT19937(fst.args({"seed": "1234"})))
        # Original manuscript used rcut = L/2, but the default rcut in spce used here is 10
        monte_carlo.set(fst.spce(fst.args({
            "cubic_box_length": str(box_length),
            "physical_constants": "CODATA2010",
            "alphaL": "5.6",
            "kmax_squared": "38",
        })))
        #monte_carlo.set(spce.system(
        #    config,
        #    rcut=config.domain().min_side_length()/2.))
        R = monte_carlo.configuration().physical_constants().ideal_gas_constant()
        monte_carlo.set(fst.MakeThermoParams(fst.args(
            {"beta": str(1./(R*temperature/1e3)),
             "chemical_potential": "1."})))
        monte_carlo.set(fst.MakeMetropolis())
        monte_carlo.add(fst.MakeTrialTranslate(fst.args({"weight": "1.", "tunable_param": "0.275"})))
        monte_carlo.add(fst.MakeTrialRotate(fst.args({"weight": "1.", "tunable_param": "50"})))
        fst.SeekNumParticles(num_particles).with_trial_add().run(monte_carlo)
        monte_carlo.add(fst.MakeLogAndMovie(fst.args({
            "steps_per" : str(steps_per),
            "file_name": "spce",
            "clear_file": "true"})))
        monte_carlo.add(fst.MakeCheckEnergyAndTune(fst.args({"steps_per" : str(steps_per),
                                                             "tolerance": str(1e-6)})))

        # equilibrate
        monte_carlo.attempt(int(1e6))

        # compute average energy using a stepper/analysis and output into file
        energy = fst.MakeEnergy(fst.args({
            "steps_per_update": "1",
            "steps_per_write": str(steps_per),
            "file_name": "spce_nvt_energy.txt"}))
        monte_carlo.add(energy)

        # compute average using this script
        energy_alt = fst.Accumulator()

        # production
        for _ in range(int(1e6)):
            monte_carlo.attempt()
            energy_alt.accumulate(monte_carlo.criteria().current_energy())

        # test that the two methods to compute average energy are the same
        self.assertAlmostEqual(energy.energy().average(), energy_alt.average(), delta=1e-6)

        # test the average against the NIST SRSW
        num = monte_carlo.configuration().num_particles()
        stdev = (energy.energy().block_stdev()**2 + (0.02*num)**2)**(1./2.)
        print(energy.energy().average(), energy_alt.average(), stdev)
        self.assertAlmostEqual(-46.82*num, energy.energy().average(), delta=8*stdev)

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

[2]:
%%time
unittest.main(argv=[''], verbosity=2, exit=False)
test_srsw (__main__.TestEwald1SPCENVT)
Compare with the reported average energy from https://doi.org/10.1063/1.476834. ...
-23844.03832883373 -23844.03832883373 29.44706507790323
CPU times: user 20min 14s, sys: 964 ms, total: 20min 15s
Wall time: 20min 13s
ok

----------------------------------------------------------------------
Ran 1 test in 1213.842s

OK
[2]:
<unittest.main.TestProgram at 0x7f254c596198>

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