Monte Carlo with rigid cluster moves
In this example, a low density and low temperature simulation of Lennard-Jones particles is conducted with and without rigid cluster moves and with and without single particle translations. A droplet is expected to form at these conditions. Without cluster moves, this droplet is unable to translate and rotate freely once it is formed. Single particle translations are still important. In order to obey detailed balance with rigid cluster moves, clusters cannot be created or destroyed. Thus, single particle translations are still required to form a droplet. Without single particle translations, particle randomly dispersed will not be allowed to coalesce.
The three possible combinations are listed as follows: - sim 1 - single particle translations only - sim 2 - rigid cluster moves only - sim 3 - both
[1]:
for sim in [1, 2, 3]: # see combinations described above
params={"sim": sim,
"potential": "Potential Model LennardJones",
"translate": "TrialTranslate tunable_param 2 tunable_target_acceptance 0.2",
"rigid": "",
"trials_per": int(1e3)}
if sim != 1:
params["potential"]="Potential Model LennardJones EnergyMap EnergyMapAll"
if sim == 2:
params["translate"]=""
if sim == 2 or sim == 3:
params["rigid"] = """
NeighborCriteria energy_maximum -0.5
TrialRigidCluster neighbor_index 0
AnalyzeCluster trials_per_write {trials_per} output_file cluster{sim}.txt""".format(**params)
script="""
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /feasst/particle/lj.fstprt
{potential}
ThermoParams beta 4 chemical_potential0 1
Metropolis
{translate}
TrialAdd particle_type 0
Run until_num_particles 25
RemoveTrial name TrialAdd
{rigid}
Tune
Log trials_per_write {trials_per} output_file log{sim}.txt
CPUTime trials_per_write {trials_per} output_file cpu{sim}.txt
Movie trials_per_write {trials_per} output_file movie{sim}.xyz clear_file True
CheckEnergy trials_per_update {trials_per} tolerance 1e-8
Run num_trials 1e5
""".format(**params)
with open('script.txt', 'w') as file: file.write(script)
import subprocess
syscode = subprocess.call("../../../build/bin/fst < script.txt > script.log", shell=True, executable='/bin/bash')
with open('script.log', 'r') as file: print(file.read(), '\n', 'exit:', syscode)
# FEASST version: v0.23.0-52-g781fe2b0e3-user/widom3
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/user/open_mab_cg/feasst/particle/lj.fstprt
Potential Model LennardJones
ThermoParams beta 4 chemical_potential0 1
Metropolis
TrialTranslate tunable_param 2 tunable_target_acceptance 0.2
TrialAdd particle_type 0
Run until_num_particles 25
# initializing random number generator with seed: 1695991527
RemoveTrial name TrialAdd
Tune
Log output_file log1.txt trials_per_write 1000
CPUTime output_file cpu1.txt trials_per_write 1000
Movie clear_file True output_file movie1.xyz trials_per_write 1000
CheckEnergy tolerance 1e-8 trials_per_update 1000
Run num_trials 1e5
exit: 0
# FEASST version: v0.23.0-52-g781fe2b0e3-user/widom3
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/user/open_mab_cg/feasst/particle/lj.fstprt
Potential EnergyMap EnergyMapAll Model LennardJones
ThermoParams beta 4 chemical_potential0 1
Metropolis
TrialAdd particle_type 0
Run until_num_particles 25
# initializing random number generator with seed: 1695991527
RemoveTrial name TrialAdd
NeighborCriteria energy_maximum -0.5
TrialRigidCluster neighbor_index 0
AnalyzeCluster output_file cluster2.txt trials_per_write 1000
Tune
Log output_file log2.txt trials_per_write 1000
CPUTime output_file cpu2.txt trials_per_write 1000
Movie clear_file True output_file movie2.xyz trials_per_write 1000
CheckEnergy tolerance 1e-8 trials_per_update 1000
Run num_trials 1e5
exit: 0
# FEASST version: v0.23.0-52-g781fe2b0e3-user/widom3
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/user/open_mab_cg/feasst/particle/lj.fstprt
Potential EnergyMap EnergyMapAll Model LennardJones
ThermoParams beta 4 chemical_potential0 1
Metropolis
TrialTranslate tunable_param 2 tunable_target_acceptance 0.2
TrialAdd particle_type 0
Run until_num_particles 25
# initializing random number generator with seed: 1695991527
RemoveTrial name TrialAdd
NeighborCriteria energy_maximum -0.5
TrialRigidCluster neighbor_index 0
AnalyzeCluster output_file cluster3.txt trials_per_write 1000
Tune
Log output_file log3.txt trials_per_write 1000
CPUTime output_file cpu3.txt trials_per_write 1000
Movie clear_file True output_file movie3.xyz trials_per_write 1000
CheckEnergy tolerance 1e-8 trials_per_update 1000
Run num_trials 1e5
exit: 0
[2]:
import unittest
# #import math
import pandas as pd
class TestRigidCluster(unittest.TestCase):
def test(self):
# end_to_end=pd.read_csv('end_to_end.txt')
# self.assertAlmostEqual(end_to_end['average'][0], math.sqrt(params['num_monomers']), delta=20)
# rg=pd.read_csv('rg.txt')
# self.assertAlmostEqual(rg['average'][0], params['num_monomers']/6, delta=200)
# With rigid cluster moves only, clusters cannot coalesce or break up.
# So the clusters in the original configuration remain.
for sim in [3]:
if sim == 3:
cluster=pd.read_csv('cluster3.txt')
print(cluster['average'])
assert cluster['average'][0] > 8
unittest.main(argv=[''], verbosity=2, exit=False)
test (__main__.TestRigidCluster) ... ok
----------------------------------------------------------------------
Ran 1 test in 0.002s
OK
0 24.602957
Name: average, dtype: float64
[2]:
<unittest.main.TestProgram at 0x7f0014f0fcd0>
Take a look at the cpu files. Sim 1 is faster than sim 2 because cluster moves are more expensive. Even when clusters are small, energy map book keeping and cluster calculations can be slow compared to a simple LJ potential.
Sim 3 is even slower than sim 2 because sim 3 forms larger clusters. Each step moves many particles and is more expensive.
Note that, even if sim 3 takes the longest time to complete a given number of steps, it can still be considered the most efficient sim due to improved sampling. Essentially, sims 1 and 2 are non ergodic and get trapped in a particular region of phase space, while sim 3 can sample all of phase space.