Compressing the System



  • How do I compress the system to a target density?


  • Describe how Variants describe values that vary with the timestep.

  • Show how to use BoxResize with a Ramp variant to slowly compress the system to the target volume.

  • Demonstrate the crystallization of Lennard-Jones particles

Boilerplate code

import math

import hoomd
import matplotlib

%matplotlib inline'ggplot')
import matplotlib_inline


The render function in the next (hidden) cell will render the system state using fresnel. Find the source in the hoomd-examples repository.


random.gsd, generated by the previous section of this tutorial has particles randomly placed in a low density fluid state. To find the crystalline state, you need to compress the system to a higher density. If you were to immediately scale the system to the target density, the particle’s hard cores would overlap and the integrator would be become unstable. A more effective strategy gradually resizes the simulation box a little at a time over the course of a simulation, allowing any slight overlaps to relax.

Initialize the simulation

Configure this simulation to run on the CPU.

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=1)

Set up the molecular dynamics simulation, as discussed in the previous sections of this tutorial:

integrator =
cell =
lj =
lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)
lj.r_cut[('A', 'A')] = 2.5
nvt =
simulation.operations.integrator = integrator


The BoxResize class performs scales the simulation box. You provide a Variant that tells BoxResize how to interpolate between two boxes, box1 and box2, over time. A Variant is a scalar-valued function of the simulation timestep. BoxResize sets box1 when the function is at its minimum, box2 when the function is at its maximum, and linearly interpolates between the two boxes for intermediate values.

Use the Ramp variant to define a function that ramps from one value to another over a given number of steps. Note that this simulation does not start at timestep 0, random.gsd was saved after running a number of steps to randomize the system in the previous section:


Create a Ramp that starts at this step and ends 100,000 steps later:

ramp = hoomd.variant.Ramp(A=0, B=1, t_start=simulation.timestep, t_ramp=20000)

Here is what the ramp looks like as a function of timestep:

steps = range(0, 40000, 20)
y = [ramp(step) for step in steps]

fig = matplotlib.figure.Figure(figsize=(5, 3.09))
ax = fig.add_subplot()
ax.plot(steps, y)
ax.tick_params('x', labelrotation=30)

Gradually resize the simulation box

Now, you need to determine the initial simulation box box1 and final simulation box at the target density box2. This is the initial number density \(\rho\) of the current simulation box:

rho = simulation.state.N_particles /

Create a new Box at the target number density:

initial_box =
final_box = hoomd.Box.from_box(initial_box)  # make a copy of initial_box
final_rho = 1.2
final_box.volume = simulation.state.N_particles / final_rho

To avoid destabilizing the integrator with large forces due to large box size changes, scale the box with a small period.

box_resize_trigger = hoomd.trigger.Periodic(10)

Construct the BoxResize updater and add it to the simulation.

box_resize = hoomd.update.BoxResize(
    box1=initial_box, box2=final_box, variant=ramp, trigger=box_resize_trigger

Run the simulation to compress the box.


At the half way point of the ramp, the variant is 0.5 and the current box lengths are half way between the initial and final ones:

ramp(simulation.timestep - 1)
current_box =
(current_box.Lx - initial_box.Lx) / (final_box.Lx - initial_box.Lx)

Finish the rest of the compression:


The current box matches the final box:

[16]: == final_box

The system is at a higher density:


Now that the system has reached the target volume, remove the BoxResize updater from the simulation:


Equilibrating the system

Run the simulation for a longer time and see if it self assembles the ordered fcc structure:

This cell may require several minutes to complete.


Here is what the final state of the system looks like:


Analyzing the simulation results

Is this the fcc structure? Did the system reach equilibrium? These are topics covered in the Introducing HOOMD-blue tutorial. Consider it an exercise to modify this example and implement that analysis here. You will need to save a GSD file with the simulation trajectory, then use freud’s SolidLiquid order parameter to identify when the ordered structure appears.

This is the end of the introducing molecular dynamics tutorial! It described molecular dynamics simulations, demonstrated how to initialize, randomize, compress, and equilibrate a system of Lennard-Jones particles. See the other HOOMD-blue tutorials to learn about other concepts, or browse the reference documentation for more information.