Packing a Dimer on a Cubic Lattice

In this tutorial, we’ll walk through filling a box with dimers using lammpsio. We’ll create a simple cubic lattice where each unit cell contains one dimer made up of two particle (one of type A and one of type B) that are bonded to each other. At the end, we will create a data file ready to be used by LAMMPS.

First, we import numpy for making the simple cubic lattice and lammpsio to handle creating the data file.

import lammpsio
import numpy

Creating the lattice

We define the core parameters of our system, including particle diameters, offset between bonded particles, and the number of repetitions of the unit cell. We choose the particles to have a unit diameter \(d\) and a bond length of \(\ell = 1.5 d\). We chose a lattice spacing of \(2 d + \ell\) so that each dimer is separated by one diameter. We also choose to place 1000 dimers total (10 in each direction).

diameter = 1.0
bond_length = 1.5
lattice_spacing = 2 * diameter + bond_length
num_repeat = [10, 10, 10]

Define the unit cell

The unit cell is a volume containing particles that can be copied along the vectors that define it in order to fill space. Our unit cell is a cube, so its vectors are the 3 Cartesian axes, scaled by the lattice spacing. Each unit cell contains two particles: type A at the origin and type B shifted along the x axis to give the proper spacing. We also specify the type IDs (A is 1, B is 2) and masses of the particles in the unit cell.

unit_cell_vectors = numpy.array([
    [1, 0, 0],
    [0, 1, 0],
    [0, 0, 1]
]) * lattice_spacing
unit_cell_coords = numpy.array([
    [0, 0, 0],
    [bond_length, 0, 0]
])
unit_cell_typeids = [1, 2]
unit_cell_mass = [1.0, 1.5]

Define the box

The simulation box is the volume obtained by repeating the unit cell the desired number of times. First, we repeat each vector by the number of repeats we specified. Then, we transpose these vectors to form the matrix [a b c] that defines a LAMMPS box. We use the from_matrix method to create our box and choose to put the lower corner of the box at the origin [0, 0, 0].

box_matrix = (unit_cell_vectors * num_repeat).T
box = lammpsio.Box.from_matrix(low=[0, 0, 0], matrix=box_matrix)

Creating the snapshot

The Snapshot holds the data about the particle’s configuration and properties, the Box, and the timestep. First, we calculate the total number of particles by multiplying the number of unit cells by the number of particles per cell.

num_cells = numpy.prod(num_repeat)
num_per_unit_cell = unit_cell_coords.shape[0]
snap = lammpsio.Snapshot(N=num_cells * num_per_unit_cell, box=box)

We then generate positions for all particles by iterating through each unit cell in our lattice. The origin of the lattice (relative to the origin of the box) is calculated and used to shift the unit cell coordinates. Finally, the origin of the box is added to give the final particle positions.

Note that lammpsio automatically allocates an array with the right data type and shape for particle data, so particle data can be assigned directly to the snapshot rather than using an intermediate array!

for i, unit_cell_idx in enumerate(numpy.ndindex(*num_repeat)):
    first = i * num_per_unit_cell
    last = first + num_per_unit_cell
    snap.position[first:last] = (
        numpy.array(unit_cell_idx) * lattice_spacing + unit_cell_coords
    )
snap.position += snap.box.low

We then create an array of type IDs by replicating our unit cell type IDs. We do the same thing to give the particles their masses.

snap.typeid = numpy.tile(unit_cell_typeids, num_cells)
snap.mass = numpy.tile(unit_cell_mass, num_cells)

To specify the bond for each dimer, we also need to add a Bonds object to the snapshot. We know that each cell contains one dimer and thus one bond. Since all of our bonds in this system are the same, we assign them all type ID 1.

snap.bonds = lammpsio.Bonds(N=num_cells, num_types=1)
snap.bonds.typeid = numpy.ones(snap.bonds.N)

Each dimer consists of consecutive particle IDs (1-2, 3-4, etc.). We create bonds by connecting these consecutive pairs.

snap.bonds.members = [[2 * i + 1, 2 * i + 2] for i in range(snap.bonds.N)]

Save to LAMMPS data file

Finally, we save the configuration to a data file ready to be used in LAMMPS. The data file is written in the molecular style since we have bonds.

lammpsio.DataFile.create(
    filename="dimer_lattice.data", snapshot=snap, atom_style="molecular"
)

<lammpsio.data.DataFile at 0x109c575f0>

Save to HOOMD-blue GSD file

HOOMD-blue is another molecular simulation tool whose user community overlaps with LAMMPS. There may come a time when you want to use HOOMD-blue or share your LAMMPS data file with someone who is more familar with it. Manually converting a LAMMPS data file to HOOMD-blue’s GSD format can be tedious — for example, HOOMD-blue requires the box to be centered at [0, 0, 0]. Luckily, lammpsio automatically handles this coordinate transformation and other format differences for you!

HOOMD-blue requires alphanumeric type names along with type IDs, so we have to add those to our snapshot using a LabelMap. By default, if you do not specify a LabelMap, lammpsio will convert type IDs to string types in the GSD file (e.g., typeID 1 becomes "1"). Here, we’re going to explicitly map particle typeID 1 -> A & 2 -> B and bond typeID 1 -> dimer. lammpsio will use the LabelMap to set the alphanumeric types in the GSD file.

Note: LAMMPS now also supports alphanumeric type labeling, but lammpsio does not currently support for this feature. It is planned as a future addition.

snap.type_label = lammpsio.LabelMap({1: "A", 2: "B"})
snap.bonds.type_label = lammpsio.LabelMap({1: "dimer"})

We can write this out to file and have the same particle and bond data ready to be used in a different simulation engine! You can similarly use Snapshot.from_hoomd_gsd to convert a HOOMD-blue GSD frame into a Snapshot!

import gsd.hoomd

with gsd.hoomd.open("dimer_lattice.gsd", "w") as f:
    f.append(snap.to_hoomd_gsd())