Simulating a System of Patchy Particles

Overview

Questions

• How do I implement pairwise energetic interactions in HPMC?

Objectives

• Run a simulation of particles interacting through the Kern–Frenkel pair potential.

• Log the potential energy during an HPMC simulation.

Boilerplate code

import itertools
import math

import gsd.hoomd
import hoomd
import matplotlib
import matplotlib_inline
import numpy

%matplotlib inline
matplotlib.style.use('ggplot')
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')

In the previous section, you learned about patchy particles and the the Kern–Frenkel model. In this section, we will initialize a system of particles, add the Kern–Frenkel pair potential to the simulation, and run the simulation while logging the energy.

Constructing the system

Construct the system using the same code you used in Introducing HOOMD-blue tutorial and then initialize a simulation.

m = 4
N_particles = 2 * m**3
spacing = 1.2
K = math.ceil(N_particles ** (1 / 3))
L = K * spacing
x = numpy.linspace(-L / 2, L / 2, K, endpoint=False)
position = list(itertools.product(x, repeat=3))
position = position[0:N_particles]
orientation = [(1, 0, 0, 0)] * N_particles

# gsd snapshot
snapshot = gsd.hoomd.Frame()
snapshot.particles.N = N_particles
snapshot.particles.position = position
snapshot.particles.orientation = orientation
snapshot.particles.typeid = [0] * N_particles
snapshot.particles.types = ['A']
snapshot.configuration.box = [L, L, L, 0, 0, 0]
with gsd.hoomd.open(name='initial.gsd', mode='x') as f:
    f.append(snapshot)

# build simulation
cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=0)
simulation.create_state_from_gsd(filename='initial.gsd')

Define interaction parameters

Define the parameters \(\delta\), \(\varepsilon\), \(\lambda\), and \(\sigma\) of the Kern–Frenkel model to simulate.

kf_delta_deg = 45  # half-opening angle of patches
kf_epsilon = 1.0  # strength of patchy interaction in kT
kf_lambda = 1.2  # range of patchy interaction
sigma = 1.0  # hard core diameter

Add HPMC integrator

Now add an HPMC integrator to the system. Add the Sphere integrator because the cores of the particles are spherical. These particles also have an orientation, so tell HPMC to perform rotation moves on the particles via the orientable key of the shape dictionary.

mc = hoomd.hpmc.integrate.Sphere()
mc.shape['A'] = dict(diameter=sigma, orientable=True)
simulation.operations.integrator = mc

Add the pair potential

Write patch code

The hoomd.hpmc.pair.user classes use strings containing C++ code to evaluate pair potentials. The code you pass in is the body of the function with the following signature:

float eval(const vec3<float>& r_ij,
           unsigned int type_i,
           const quat<float>& q_i,
           float d_i,
           float charge_i,
           unsigned int type_j,
           const quat<float>& q_j,
           float d_j,
           float charge_j
)

The Kern–Frenkel potential only depends on the center-to-center vector r_ij and the orientations of the two particles, represented as the quaternions q_i and q_j. Users have access to other variables, which can be useful for more complicated potentials. Note the type of r_ij, q_i, and q_j: vec3 and quat are classes that represent 3D vectors and quaternions, respectively, and are defined in the HOOMD vector math library. The HOOMD vector math library also contains operations you can use on vectors, e.g., you will use rotate to rotate the particles into the global frame of reference to evaluate the pair potential.

At the beginning of the body of the function, define some constants:

patch_code = f"""
const float delta = {kf_delta_deg} * M_PI / 180;  // delta in radians
const float epsilon = {kf_epsilon:f};
const float lambda = {kf_lambda:f};
const float sigma = {sigma:f};  // hard core diameter
"""

Now set the strength \(\varepsilon\) of the patchy interactions. To do so, first define the temperature \(\beta^{-1} = kT\) of the system. HPMC assumes \(\beta = 1\), so to change the temperature of a system, you must change the relative strength of the energetic interactions (\(\varepsilon\) in this case).

Use param_array, which is an array of floats that is declared in the patch kernel that you can use to change the value of variables in the patch code. Later in this section, you will change the temperature of the system by changing the value of param_array[0].

patch_code += """
const float kT = param_array[0];
const float beta_epsilon = epsilon/kT;
"""

Now compute the patch director in the particle frame of reference (ehat_particle_reference_frame). Use rotate to compute \(\hat{e}_i\) and \(\hat{e}_j\) by rotating ehat_particle_reference_frame by the quaternions q_i and q_j, respectively.

patch_code += """
const vec3<float> ehat_particle_reference_frame(1, 0, 0);
vec3<float> ehat_i = rotate(q_i, ehat_particle_reference_frame);
vec3<float> ehat_j = rotate(q_j, ehat_particle_reference_frame);
"""

Check to see if the patches on each particle are “facing” the other particle. In other words, evaluate the orientational masking function \(f(\Omega_i, \Omega_j)\) defined in the previous section.

patch_code += """
vec3<float> r_hat_ij = r_ij / sqrtf(dot(r_ij, r_ij));
bool patch_on_i_is_aligned_with_r_ij = dot(ehat_i, r_hat_ij) >= cos(delta);
bool patch_on_j_is_aligned_with_r_ji = dot(ehat_j, -r_hat_ij) >= cos(delta);
"""

All that’s left now is to actually evaluate the pair potential. In HOOMD, the patch kernel is evaluated after any overlap checks, so we do not need to check for the hard sphere overlap (\(r_{ij} < \sigma_{ij}\)). Instead, we only need to check if they are within range of the patchy interaction; this condition simplifies to \(|\vec{r}_{ij}| < \lambda\sigma\). For clarity, evaluate the expression as it is written in the previous section, but note that you can avoid the square root calculation by evaluating \(|\vec{r}_{ij}|^2 = \vec{r}_{ij} \cdot \vec{r}_{ij} < (\lambda\sigma)^2\).

Note the use of double curly brackets {{..}} for the if..else code block; Python evaluates expressions inside of curly brackets in f-strings, so we must include a second set of curly brackets so that one set of curly brackets remain after evaluation of the f-string.

patch_code += """
float rsq = dot(r_ij, r_ij);
float r_ij_length = sqrtf(rsq);
if (patch_on_i_is_aligned_with_r_ij
    && patch_on_j_is_aligned_with_r_ji
    && r_ij_length < lambda*sigma)
    {
    return -beta_epsilon;
    }
else
    {
    return 0.0;
    }
"""

Set cutoff distance for the patchy interaction

Before creating the pair potential object to add to the Simulation, calculate r_cut. Similar to how HOOMD uses neighbor lists in MD to increase performance by only calculating interactions between nearby particles, the hoomd.hpmc.pair.user classes take an r_cut parameter, which is the center to center distance cutoff beyond which all pair interactions are 0. Here, that distance is the hard core diameter \(\sigma\) plus the radial extent of the patchy interaction \(\sigma(\lambda-1)\).

r_cut = sigma + sigma * (kf_lambda - 1)

Set variable parameters for patch kernel

initial_kT = 3.0
patch_param_array = [initial_kT]

Adding the patchy potential to the simulation

Now create a hoomd.hpmc.pair.user.CPPPotential object that will evaluate the Kern–Frenkel potential in the simulation. Use the patch_code string that you built up that serves as the body of the function that computes patch interaction energy. Initialize the temperature to \(kT = 3\) to randomize the system on the first call to Simulation.run().

initial_kT = 3.0
patch_potential = hoomd.hpmc.pair.user.CPPPotential(
    r_cut=r_cut, code=patch_code, param_array=[initial_kT]
)

Finally, attach patch_potential to the HPMC integrator so that the integrator evaluates the pair potential.

mc.pair_potential = patch_potential

You can query the energy associated with the patch potential via the energy property of the patch potential object (after running for 0 or more steps). The particles are not interacting in the initial configuration, and hence the initial value of the energy is zero.

simulation.run(0)
patch_potential.energy

0.0

Log the patch energy

Create a logger and add the patchy potential to it, and then create a GSD writer to write the log values to file during the simulation. Also add the temperature to the logger so that later you can compare energies from parts of the simulation at different temperatures.

logger = hoomd.logging.Logger()
logger.add(
    patch_potential,
    quantities=[
        'energy',
    ],
)
logger[('kT',)] = (lambda: patch_potential.param_array[0], 'scalar')
gsd_writer = hoomd.write.GSD(
    filename='log.gsd',
    trigger=hoomd.trigger.Periodic(10),
    mode='xb',
    filter=hoomd.filter.Null(),
    logger=logger,
)
simulation.operations.writers.append(gsd_writer)

Run the simulation

Randomize the system

Run for a bit at the high initial temperature to randomize the system.

simulation.run(500)

Change the temperature

Now run for a bit at a lower temperature so that the patchy interactions cause the particles to cluster. Recall, the param_array property of the patch potential is the effective temperature of the system.

patch_potential.param_array[0] = 0.1
simulation.run(2000)

Now, increase the temperature again so that the clusters break apart.

patch_potential.param_array[0] = initial_kT
simulation.run(2000)

As discussed in the Logging tutorial, flush the write buffer (this is not necessary in typical workflows).

gsd_writer.flush()

Analyzing the results

Plot the energy of the system to analyze the results.

Use the gsd package to read just the log values from the GSD file.

data = gsd.hoomd.read_log('log.gsd')
step = data['configuration/step']

Recall, the the value of patch_potential.energy is really the energy divided by the thermal energy scale \(kT\), so multiply the energy by \(kT\) to get the actual energy.

beta_potential_energy = data['log/hpmc/pair/user/CPPPotential/energy']
kT = data['log/kT']
potential_energy = beta_potential_energy * kT

fig = matplotlib.figure.Figure(figsize=(5, 3.09))
ax = fig.add_subplot()
(ll,) = ax.plot(step, potential_energy)
ax.set_xlabel('timestep')
ax.set_ylabel('potential energy', c=ll.get_c())

ax2 = ax.twinx()
ax2.grid(False)
ax2._get_lines.get_next_color()
(ll,) = ax2.plot(step, kT)
ylabel = ax2.set_ylabel('kT', c=ll.get_c())

fig

Notice the decrease in energy when the temperature was at the lower setting, indicative of particle clustering.

Summary

In this tutorial, you learned about patchy particles and the Kern–Frenkel model of patchy particles, and you implemented the Kern–Frenkel model in HPMC. You observed the change in potential energy upon changing the temperature in a patchy particle simulation.

For more information on simulating patchy particles, or more generally HPMC simulations with energetic interactions, see the HOOMD-blue documentation on pair potentials in HPMC.