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

for sim in [1, 2, 3]: # see combinations described above
    params={"sim": sim,
            "potential": "Potential Model=LennardJones",
            "translate": "TrialTranslate tunable_param=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}_cluster.csv""".format(**params)
    script="""
MonteCarlo
Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt
{potential}
ThermoParams beta=4 chemical_potential=1
Metropolis
{translate}
TrialAdd particle_type=lj
Run until_num_particles=25
Remove name=TrialAdd
{rigid}
Tune
Let [write]=trials_per_write {trials_per} output_file cluster{sim}
Log [write].csv
CPUTime [write]_cpu.txt
Movie [write].xyz clear_file True
CheckEnergy trials_per_update {trials_per} decimal_places 6
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)

# Usage: /Users/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt
Potential Model=LennardJones
ThermoParams beta=4 chemical_potential=1
Metropolis
TrialTranslate tunable_param=2
TrialAdd particle_type=lj
Run until_num_particles=25
# Initializing random number generator with seed: 1749409027
Remove name=TrialAdd
Tune
Log output_file=cluster1.csv trials_per_write=1000
CPUTime output_file=cluster1_cpu.txt trials_per_write=1000
Movie clear_file=True output_file=cluster1.xyz trials_per_write=1000
CheckEnergy decimal_places=6 trials_per_update=1000
Run num_trials=1e5

 exit: 0
# Usage: /Users/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt
Potential EnergyMap=EnergyMapAll Model=LennardJones
ThermoParams beta=4 chemical_potential=1
Metropolis
TrialAdd particle_type=lj
Run until_num_particles=25
# Initializing random number generator with seed: 1749409028
Remove name=TrialAdd
NeighborCriteria energy_maximum=-0.5
TrialRigidCluster neighbor_index=0
AnalyzeCluster output_file=cluster2_cluster.csv trials_per_write=1000
Tune
Log output_file=cluster2.csv trials_per_write=1000
CPUTime output_file=cluster2_cpu.txt trials_per_write=1000
Movie clear_file=True output_file=cluster2.xyz trials_per_write=1000
CheckEnergy decimal_places=6 trials_per_update=1000
Run num_trials=1e5

 exit: 0
# Usage: /Users/user/feasst/build/bin/fst < file.txt
FEASST version 0.25.13
MonteCarlo
Configuration cubic_side_length=8 particle_type=lj:/feasst/particle/lj_new.txt
Potential EnergyMap=EnergyMapAll Model=LennardJones
ThermoParams beta=4 chemical_potential=1
Metropolis
TrialTranslate tunable_param=2
TrialAdd particle_type=lj
Run until_num_particles=25
# Initializing random number generator with seed: 1749409028
Remove name=TrialAdd
NeighborCriteria energy_maximum=-0.5
TrialRigidCluster neighbor_index=0
AnalyzeCluster output_file=cluster3_cluster.csv trials_per_write=1000
Tune
Log output_file=cluster3.csv trials_per_write=1000
CPUTime output_file=cluster3_cpu.txt trials_per_write=1000
Movie clear_file=True output_file=cluster3.xyz trials_per_write=1000
CheckEnergy decimal_places=6 trials_per_update=1000
Run num_trials=1e5

 exit: 0

import numpy as np
# #import math
import pandas as pd

#         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_cluster.csv')
        print(cluster['average'])
        assert cluster['average'][0] > 8

0    24.298629
Name: average, dtype: float64

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.