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

[2]:
import unittest
import feasst as fst

for sim in [1, 2, 3]: # see combinations described above
    monte_carlo = fst.MonteCarlo()
    monte_carlo.add(fst.Configuration(fst.MakeDomain(fst.args({"cubic_box_length": "8"})),
         fst.args({"particle_type": fst.install_dir() + "/forcefield/data.lj"})))
    if sim == 1:
        monte_carlo.add(fst.MakePotential(fst.MakeLennardJones()))
    else:
        monte_carlo.add(fst.MakePotential(fst.MakeLennardJones(),
                        fst.MakeVisitModel(fst.MakeVisitModelInner(fst.MakeEnergyMapAll()))))
    monte_carlo.set(fst.MakeThermoParams(fst.args({"beta": "4", "chemical_potential": "1."})))
    monte_carlo.set(fst.MakeMetropolis())
    if sim == 1 or sim == 3:
        monte_carlo.add(fst.MakeTrialTranslate(fst.args(
            {"tunable_param": "2.", "tunable_target_acceptance": "0.2"})))
    fst.SeekNumParticles(25).with_trial_add().run(monte_carlo)

    # record the size of the cluster of the first particle for comparison
    if sim == 2 or sim == 3:
        monte_carlo.add(fst.MakeNeighborCriteria(fst.args({'energy_maximum': '-0.5'})))
        select_cluster = fst.SelectCluster(fst.args({"neighbor_index": "0"}))
        select_cluster.select_cluster(0, monte_carlo.system())
        cluster_size = select_cluster.mobile().num_particles()
        monte_carlo.add(fst.MakeTrialRigidCluster(fst.args({"neighbor_index": "0"})))
    steps_per = int(1e3)
    monte_carlo.add(fst.MakeTuner(fst.args({"steps_per" : str(steps_per)})))
    monte_carlo.add(fst.MakeLog(fst.args({"steps_per" : str(steps_per)})))
    monte_carlo.add(fst.MakeCPUTime(fst.args({"steps_per" : str(steps_per),
                                                    'file_name': 'cpu'+str(sim)+'.txt'})))
    moviename='movie'+str(sim)+'.xyz'
    monte_carlo.add(fst.MakeMovie(fst.args(
        {"steps_per" : str(steps_per), "file_name" : moviename, "clear_file": "True"})))
    monte_carlo.add(fst.MakeCheckEnergy(fst.args(
        {"steps_per" : str(steps_per), "tolerance" : "1e-8"})))
    monte_carlo.attempt(int(1e5))

    # With rigid cluster moves only, clusters cannot coalesce or break up.
    # So the clusters in the original configuration remain.
    if sim == 2 or sim == 3:
        select_cluster.select_cluster(0, monte_carlo.system())
        if sim == 2:
            assert(select_cluster.mobile().num_particles() == cluster_size)
        # otherwise, we expect the formation of a droplet
        if sim == 3:
            cluster_size_av = monte_carlo.trial(1).stage(0).trial_select().printable("cluster_size").average()
            assert(cluster_size_av > 20)

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.