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": ""}
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"""
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}
set_variable trials_per 1e3
Tune
Log trials_per_write trials_per file_name log{sim}.txt
CPUTime trials_per_write trials_per file_name cpu{sim}.txt
Movie trials_per_write trials_per file_name movie{sim}.xyz clear_file True
CheckEnergy trials_per_update trials_per tolerance 1e-8
AnalyzeCluster trials_per_write trials_per file_name cluster{sim}.txt
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.21-21-g12898a9d3b-dirty-hwh/trial_library
Configuration cubic_side_length 8 particle_type0 /home/hwh/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: 1678985917
RemoveTrial name TrialAdd
Tune
Log file_name log1.txt trials_per_write 1e3
CPUTime file_name cpu1.txt trials_per_write 1e3
Movie clear_file True file_name movie1.xyz trials_per_write 1e3
CheckEnergy tolerance 1e-8 trials_per_update 1e3
AnalyzeCluster file_name cluster1.txt trials_per_write 1e3
Run num_trials 1e5
exit: 0
# FEASST version: v0.21-21-g12898a9d3b-dirty-hwh/trial_library
Configuration cubic_side_length 8 particle_type0 /home/hwh/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: 1678985918
RemoveTrial name TrialAdd
NeighborCriteria energy_maximum -0.5
TrialRigidCluster neighbor_index 0
Tune
Log file_name log2.txt trials_per_write 1e3
CPUTime file_name cpu2.txt trials_per_write 1e3
Movie clear_file True file_name movie2.xyz trials_per_write 1e3
CheckEnergy tolerance 1e-8 trials_per_update 1e3
AnalyzeCluster file_name cluster2.txt trials_per_write 1e3
Run num_trials 1e5
exit: 0
# FEASST version: v0.21-21-g12898a9d3b-dirty-hwh/trial_library
Configuration cubic_side_length 8 particle_type0 /home/hwh/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: 1678985918
RemoveTrial name TrialAdd
NeighborCriteria energy_maximum -0.5
TrialRigidCluster neighbor_index 0
Tune
Log file_name log3.txt trials_per_write 1e3
CPUTime file_name cpu3.txt trials_per_write 1e3
Movie clear_file True file_name movie3.xyz trials_per_write 1e3
CheckEnergy tolerance 1e-8 trials_per_update 1e3
AnalyzeCluster file_name cluster3.txt trials_per_write 1e3
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 13.880786
Name: average, dtype: float64
[2]:
<unittest.main.TestProgram at 0x7fe3ba374400>
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.