Grand canonical ensemble transition-matrix Monte Carlo with RPMΒΆ
This example is similar to the SPC/E tutorial, except this time we simulate the RPM model.
[1]:
import unittest
import feasst as fst
class TestFlatHistogramRPM(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 previously
published values: https://doi.org/10.1063/1.5123683
"""
monte_carlo = fst.MonteCarlo()
monte_carlo.set(fst.rpm(fst.args({
"cubic_box_length": "12",
"cutoff": "4.891304347826090",
"alphaL": "6.87098396396261"})))
monte_carlo.add_to_reference(fst.MakePotential(fst.MakeDontVisitModel()))
temperature = 0.047899460618081;
beta_mu = -13.94;
monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": str(1/temperature),
"chemical_potential0": str(beta_mu*temperature),
"chemical_potential1": str(beta_mu*temperature)})))
criteria = fst.MakeFlatHistogram(
fst.MakeMacrostateNumParticles(
fst.Histogram(fst.args({"width": "1", "max": "2", "min": "0"})),
fst.args({"particle_type": "0"})),
fst.MakeTransitionMatrix(fst.args({"min_sweeps": "100"})))
monte_carlo.set(criteria)
monte_carlo.add(fst.MakeTrialTranslate(fst.args({"weight": "0.25", "tunable_param": "0.1"})))
monte_carlo.add(fst.MakeTrialTransferMultiple(fst.args({
"weight": "4",
"particle_type0": "0",
"particle_type1": "1",
"reference_index": "0"})))
steps_per = str(int(1e5))
monte_carlo.add(fst.MakeCriteriaUpdater(fst.args({"steps_per": steps_per})))
monte_carlo.add(fst.MakeCriteriaWriter(fst.args({"steps_per": steps_per, "file_name": "crit.txt"})))
monte_carlo.add(fst.MakeLogAndMovie(fst.args({"steps_per": steps_per, "file_name": "rpm"})))
monte_carlo.add(fst.MakeCheckEnergyAndTune(fst.args({"steps_per": steps_per})))
energy = fst.MakeEnergy(fst.args({"file_name": "rpm_fh_energy",
"steps_per_update": "1",
"steps_per_write": steps_per,
"multistate": "true"}))
monte_carlo.add(energy)
monte_carlo.run_until_complete()
lnpi_previous = [
[-1.2994315780357, 0.05],
[-1.08646312498868, 0.05],
[-0.941850889679828, 0.05]]
energy_previous = [
[0, 1e-14],
[-0.939408, 0.02],
[-2.02625, 0.04]]
for macro in range(criteria.num_states()):
self.assertAlmostEqual(
lnpi_previous[macro][0],
criteria.bias().ln_prob().value(macro),
delta=5.*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=10*stdev
)
[2]:
%%time
unittest.main(argv=[''], verbosity=2, exit=False)
test_serial_5max (__main__.TestFlatHistogramRPM)
Compare the free energies and potential energies with the previously ...
CPU times: user 3min 30s, sys: 253 ms, total: 3min 30s
Wall time: 3min 30s
ok
----------------------------------------------------------------------
Ran 1 test in 210.298s
OK
[2]:
<unittest.main.TestProgram at 0x7fda8c209cf8>
Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!