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

In this example, we reproduce the average energy reported in

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

        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
            "cubic_box_length": str(box_length),
            "physical_constants": "CODATA2010",
            "alphaL": "5.6",
            "kmax_squared": "38",
        #    config,
        #    rcut=config.domain().min_side_length()/2.))
        R = monte_carlo.configuration().physical_constants().ideal_gas_constant()
            {"beta": str(1./(R*temperature/1e3)),
             "chemical_potential": "1."})))
        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"})))
            "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

        # 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"}))

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

        # production
        for _ in range(int(1e6)):

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

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

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

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

Ran 1 test in 1213.842s

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