Molecular Dynamics: Periodic Boundary Conditions
Imagine you’re building a simulation of atoms. Think pool-ball type atoms packed into a cube. You start by creating a cube structure of 100x100x100 atoms. Of the entire 1,000,000 atom volume, roughly 58,800 atoms will lie on the surface of the cube. That’s about 6% of all the atoms in the simulation. Since these surface atoms aren’t surrounded by the ideal number of neighboring atoms they typically have higher energy than those in the bulk and will exhibit surface effects. This is true for crystals, liquids, and gasses alike in that atoms on the boundary of a phase may behave differently than those in the bulk.
If your goal is to model bulk properties of a molecule/atom and a large portion of “surface atoms” exist in a simulation, then the resultant properties will be influenced by undesirable surface effects.
One approach is to conduct simulations with a larger and larger number of particles. The more atoms there are in the simulation volume, the fewer of them will occupy the surface as opposed to the bulk. If we use the same example of a cubic structure but with a 1000x1000x1000 atom volume, the surface-to-volume ratio drops to 0.6%! But this is a dangerous route because computational power (and time/statistics/etc.) will be quickly consumed if we have to keep scaling up simulations to remove surface effects…
Periodic Boundary Condition
Imposing Periodic Boundary Conditions (PBCs) is an alternative, and preferred, way of solving the surface-effects issue. There are multiple different approaches to PBCs, and I’ll be using the minimum-image convention. To start the discussion on PBCs we have to be familiar with the concept of a unit cell. A unit cell is the most simplified representation of a system from which we can build. If we’re simulating a crystal, we might choose a small cell of several hundred atoms in the desired crystal form. If we’re simulating a gas, our unit cell might be a small volume with several hundred gas molecules. While the goal of the simulation might be to gain insight into bulk crystal or gas (easily >10¹⁰ molecules) properties, we can begin with a small unit cell volume.
Then, we can create neighboring copies of the unit cell, called images, that repeat the contents of the unit cell in the adjacent volumes.
The images are a replica of the original simulation area and serve to reduce or eliminate boundary effects by providing every atom, regardless of position in the unit cell, an equivalent surrounding environment of atoms as if it were in the bulk.
Even though there appears to now be x9 atoms to consider, we will only evolve the original unit cell with updated positions/forces/velocities/etc. The surrounding images cells will be updated as mirror replicas of the unit cell. In this way, the atoms in the image cells have no physical meaning on their own and are only constructs created for PBCs.
Image Atom Position Calculation
The procedure for finding atom’s positions in their image location can be thought of as translating their positions along the square axes(or cubic axes if working in 3D). In 2D one could expect that for a single cubic unit cell there should be 3x3–1 = images (as shown above). If the simulation is performed in 3D then it would require all combinations of x,y,z axes and, lead to at least 3x3x3–1 = 26 images.
As simple algorithm for calculating all the image atoms from a single atom in 2D with a cubic unit cell could be as follows:
for x_shift: -1 to 1
for y_shift: -1 to 1 if x_shift == 0 and y_shift == 0
continue
else
image_atom_y_pos = atom_y_pos + unit_celll_length*y_shift
image_atom_x_pos = atom_x_pos + unit_celll_length*x_shift
Implications of PBCs: Image Atom Interactions
Simulations in MD almost always require using a cutoff radius for calculating interatomic forces — restricting interactions between atoms if they are separated by further than a pre-defined cutoff distance. When using PBCs, this often means that atoms at the edge of the unit cell interact with an image atom rather than an atom in the unit cell. For example, if we take a minimized representation of one of the above illustrations with only two atoms and set the cutoff radius to roughly half the unit volume, we can create the situation shown below:
Here, the blue atom’s cutoff radius would restrict it from interacting with the original orange atom in the unit cell. However, once PBC’s are applied the blue atom will interact with the image atom in the right-most middle image because it lies within the cutoff radius.
But what if the cutoff radius were three times as large? Then the cutoff radius might not exclude the blue atom’s images in the neighbor cells, leading to an unphysical self-interaction, sometimes called an artifact. So when applying PBCs it is extremely important to consider the size of the unit cell and the cutoff radius for calculating intermolecular interactions to avoid self-interactions and artifacts. A safe precaution is to ensure that the cutoff radius is at most half of the unit cell length.
Lastly, you are familiar with neighbor lists as a means to speeding up molecular dynamics simulations, the use of PBCs and image atoms simply means that the neighbor-list needs to include either: (1) the explicit image atom position information, or (2) the original unit atom and the transformation needed to reach its image.
Implications of PBCs: Removal of Surface Atoms
Now that the basics of PBCs have been covered, let’s return to how exactly PBCs remove the previously mentioned surface effects. We’ll start with a 2D representation of a simple cubic crystal. In the unit cell of this 2D projection, there are 4 atoms packed into a square area. In a bulk representation of this cubic structure (left image below), we can see that in the bulk we would except each atom to have exactly 4 nearest neighbor atoms. However if we ran a simulation with only the unit cell, each atom would have only two neighbors so we would have four atoms displaying surface effects…
However, if we apply PBCs and create a series of image atoms in adjacent cells, each atom immediately has the ideal number of 4 nearest neighbors, eliminating any “surface” atoms!
Implications of PBCs: Perfect Tiling
One of the simple and beautiful effects of PBCs is the seeming “never-ending” nature of the unit cell, that is, when an atom exits through a wall in the unit cell it subsequently enters on the opposite side of the unit cell with the same velocity. This continuity is supported by the layout of the image cells because as an atom leaves the unit cell, that same atom’s image can be observed to enter the unit cell from an image. The effect is shown below:
This “perfect tiling” effect can be ensured by translating the positions of the atoms by the unit cell length as they cross the unit cell wall boundaries, so their new position is inside the unit cell. This type of wall is enforced using a simple algorithm like the one written below— assuming we have a square unit cell.
for x,y,z position in atom_coordinates if ( xyz_position > wall_max)
position = position - unit_cell_length if ( xyz_position < wall_min)
position = position + unit_cell_length
If we initialize a group of non-interacting atoms with random velocities, we can see how the atoms enter and exit the central unit cell, creating a seemingly “infinite” flow. The atoms in the unit cell are colored with blue while the image atoms are green.
Using the same pattern, we can highlight how each image is a replica of the unit cell. Note the matching colors of the image atoms translated exactly the distance of the unit cell in the x,y, directions.
Code for tiling effect demonstration and gif
Thanks for reading — any suggestions or corrections are welcome! Please reach out to me at mediumCameron@gmail.com
References
- Frenkel, Daan, and Berend Smit. Understanding Molecular Simulation from Algorithms to Applications. Academic Press, 2002.
- Braun, Efrem, et al. “Best Practices for Foundations in Molecular Simulations [Article v1.0].” Living Journal of Computational Molecular Science, vol. 1, no. 1, 2019, doi:10.33011/livecoms.1.1.5957.
