Equilibrating the System

Overview

Questions

  • What is equilibration?

  • How do I save simulation results?

Objectives

  • Explain the process of equilibration.

  • Demonstrate using GSD to write the simulation trajectory to a file.

  • Demonstrate best practices for move size tuning using Before and And Triggers.

Boilerplate code

[1]:
import hoomd
import math

The render function in the next (hidden) cell will render a snapshot using fresnel.

This is not intended as a full tutorial on fresnel - see the fresnel user documentation if you would like to learn more.

Equilibration

So far, this tutorial has placed N non-overlapping octahedra randomly in a box and then compressed it to a moderate packing fraction. The resulting configuration of particles is valid, but strongly dependent on the path taken to create it. There are many more equilibrium configurations in the set of possible configurations that do not depend on the path. Equilibrating the system is the process of taking an artificially prepared state and running a simulation. During the simulation run, the system will relax to equilibrium. Initialize the Simulation first:

[3]:
cpu = hoomd.device.CPU()
sim = hoomd.Simulation(device=cpu)
mc = hoomd.hpmc.integrate.ConvexPolyhedron(seed=2)
mc.shape['octahedron'] = dict(vertices=[(-0.5, 0, 0),
                                         (0.5, 0, 0),
                                         (0, -0.5, 0),
                                         (0, 0.5, 0),
                                         (0, 0, -0.5),
                                         (0, 0, 0.5)])
sim.operations.integrator = mc

The previous section of this tutorial wrote the compressed system to compressed.gsd. Initialize the system state from this file:

[4]:
sim.create_state_from_gsd(filename='compressed.gsd')

Here is what the compressed random initial state looks like:

[5]:
render(sim.state.snapshot)
[5]:
../../_images/tutorial_00-Introducing-HOOMD-blue_06-Equilibrating-the-System_9_0.png

Writing simulation trajectories

Save the system state to a file periodically so that you can observe the equilibration process. This tutorial previously used GSD files to store a single frame of the system state using either the GSD Python package or GSD.write. The GSD Analyzer (another operation) will create a GSD file with many frames in a trajectory.

[6]:
gsd = hoomd.write.GSD(filename='trajectory.gsd',
                      trigger=hoomd.trigger.Periodic(10000),
                      overwrite=True)
sim.operations.writers.append(gsd)

Tuning the trial move size

The previous section used the MoveSize tuner regularly during compression to adjust d and a to achieve a target acceptance ratio while the system density changed rapidly. Use it again during the equilibration run to ensure that HPMC is working optimally.

Move sizes should be tuned briefly at the beginning, then left constant for the duration of the run. Changing the move size throughout the simulation run violates detailed balance and can lead to incorrect results. Trigger the tuner every 100 steps but only for the first 5000 steps of the simulation by combining a Periodic and Before trigger with an And operation. Before returns True for all time steps t < value and the And trigger returns True when all of its child triggers also return True.

[7]:
tune = hoomd.hpmc.tune.MoveSize.scale_solver(moves=['a', 'd'],
                                             target=0.2,
                                             trigger=hoomd.trigger.And(
                                                 [hoomd.trigger.Periodic(100),
                                                  hoomd.trigger.Before(sim.timestep + 5000)]))
sim.operations.tuners.append(tune)
[8]:
sim.run(5000)

Check the acceptance ratios over the next 100 steps to verify that the tuner achieved the target acceptance ratios:

[9]:
sim.run(100)
[10]:
rotate_moves = mc.rotate_moves
mc.rotate_moves[0] / sum(mc.rotate_moves)
[10]:
0.20421552205471805
[11]:
translate_moves = mc.translate_moves
mc.translate_moves[0] / sum(mc.translate_moves)
[11]:
0.19965096686191125

Equilibrating the system

To equilibrate the system, run the simulation. The length of the run needed is strongly dependent on the particular model, the system size, the density, and many other factors. Hard particle Monte Carlo self-assembly often takes tens of millions of time steps for systems with ~10,000 particles. This system is much smaller, but it still takes approximately 1-2 million steps.

This cell may require hours to complete.

[12]:
sim.run(2e6)

Here is the final state of the system after the run.

[13]:
render(sim.state.snapshot)
[13]:
../../_images/tutorial_00-Introducing-HOOMD-blue_06-Equilibrating-the-System_22_0.png

Is the final state an equilibrium state? The next section in this tutorial shows you how to analyze the trajectory and answer this question.