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
Remove 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} 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: ./fst < file.txt
FEASST version 0.25.8
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/hwh/nightly/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: 1739454276
Remove 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 decimal_places 6 trials_per_update 1000
Run num_trials 1e5

 exit: 0
# Usage: ./fst < file.txt
FEASST version 0.25.8
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/hwh/nightly/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: 1739454277
Remove 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 decimal_places 6 trials_per_update 1000
Run num_trials 1e5

 exit: 0
# Usage: ./fst < file.txt
FEASST version 0.25.8
MonteCarlo
Configuration cubic_side_length 8 particle_type0 /home/hwh/nightly/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: 1739454277
Remove 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 decimal_places 6 trials_per_update 1000
Run num_trials 1e5

 exit: 0

[2]:

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

0    24.638376
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.