Initializing Molecular Simulations¶
When performing molecular simulations of solid materials, it is often useful to initialize your system in a crystal structure. parsnip makes this extremely easy!
HOOMD-Blue¶
HOOMD-Blue can operate directly on array data, so we can move data directly from parsnip to the simulation itself.
>>> import hoomd
>>> from parsnip import CifFile
>>> filename = "example_file.cif"
>>> cif = CifFile(filename)
>>> snapshot = hoomd.Snapshot()
>>> snapshot.particles.N = len(cif.build_unit_cell())
>>> snapshot.particles.position[:] = cif.build_unit_cell()
>>> snapshot.configuration.box = cif.box
>>> snapshot.particles.types = ["Particle"]
>>> snapshot.replicate(nx=2, ny=2, nz=3) # 2 x 2 x 3 supercell
>>> assert snapshot.particles.N == (2 * 2 * 3) * len(pos)
Once the snapshot is constructed, it can be attached to a simulation as follows:
>>> import hoomd
>>> simulation = hoomd.Simulation(device=hoomd.device.CPU())
>>> simulation.create_state_from_snapshot(snapshot)
If we want to extract additional data for our simulation, there are a few extra steps.
In HOOMD-Blue, particle.types are unique string identifiers that get mapped to
individual particles via the particles.typeid array. The following code extracts
_atom_site_type_label data and assigns the “Cu” atom type to all particles. For
structures with multiple atom sites, the particles.typeid array will have nonzero
entries that correspond with other type labels.
>>> from collections import defaultdict
>>> labels, pos = cif.build_unit_cell(additional_columns=["_atom_site_type_symbol"])
>>> # ... initialize the snapshot's `N`, `box`, and `position` data as above
>>> particle_type_map = defaultdict(list)
>>> for i, label in enumerate(labels.squeeze(axis=1)):
... particle_type_map[label].append(i)
>>> particle_type_map["Cu"] # Atoms 1-4 have the type symbol "Cu"
[0, 1, 2, 3]
>>> # Construct the TypeIDs that map our atomic symbol to the corresponding position
>>> typeid_array = np.ones(len(snapshot.particles.position), dtype=int)
>>> for typeid, label in enumerate(particle_type_map.keys()):
... typeid_array[particle_type_map[label]] = typeid
>>> snapshot.particles.typeid[:] = typeid_array
>>> snapshot.particles.typeid
array([0, 0, 0, 0])
>>> snapshot.particles.types = [str(key) for key in particle_type_map.keys()]
>>> snapshot.particles.types
['Cu']
>>> assert len(snapshot.particles.types) == len(cif["_atom_site_type_symbol"])
LAMMPS¶
In contrast to HOOMD-Blue, LAMMPS typically requires us to write out structure data to a LAMMPS Data File before simulations can begin. Although topology data is not commonly stored in CIF files, parsnip makes it simple to reconstruct atomic crystals in LAMMPS.
>>> from collections import defaultdict
>>> def write_lammps_data(cif: CifFile):
... """Convert a CIF file into a LAMMPS data file."""
... data = "(LAMMPS Data File, written with parsnip)\n\n"
...
... fractional_coordinates = cif.build_unit_cell()
... atomic_positions = fractional_coordinates @ cif.lattice_vectors.T
...
... atom_types = cif["_atom_site_type_symbol"]
... particle_type_map = defaultdict(list)
...
... # Write out the number of atoms and atom types
... data += f"{len(atomic_positions)} atoms\n"
... data += f"{len(atom_types)} atom types\n\n"
...
... # Write out the box, including the (zero, in this case) tilt factors
... lx, ly, lz, xy, xz, yz = cif.box
... data += f"0.0 {lx:.12f} xlo xhi\n"
... data += f"0.0 {ly:.12f} ylo yhi\n"
... data += f"0.0 {lz:.12f} zlo zhi\n"
... data += f"{xy:.12f} {xz:.12f} {yz:.12f} xy xz yz\n\n"
...
... # Write out the atomic position data -- note the similarities with typeid!
... data += f"Atoms # atomic\n\n"
...
... for i, label in enumerate(labels.squeeze(axis=1)):
... particle_type_map[label].append(i)
...
... # Construct the TypeIDs that map our atomic symbol to an index
... atom_type_array = np.ones(len(atomic_positions), dtype=int)
... for typeid, label in enumerate(particle_type_map.keys()):
... atom_type_array[particle_type_map[label]] = typeid
...
... for i, coordinate in enumerate(atomic_positions):
... coord_str = " ".join([f"{xyz:.12f}" for xyz in coordinate])
... data += f" {i} {atom_type_array[i]} {coord_str}\n"
...
... return data
>>> print(write_lammps_data(cif))
(LAMMPS Data File, written with parsnip)
4 atoms
1 atom types
0.0 3.600000000000 xlo xhi
0.0 3.600000000000 ylo yhi
0.0 3.600000000000 zlo zhi
0.000000000000 0.000000000000 0.000000000000 xy xz yz
Atoms # atomic
0 0 0.000000000000 0.000000000000 0.000000000000
1 0 0.000000000000 1.800000000000 1.800000000000
2 0 1.800000000000 0.000000000000 1.800000000000
3 0 1.800000000000 1.800000000000 0.000000000000