Boundary conditions (BCs)

Learning goals:

• To understand the concept of boundary conditions (PBC) and practical aspects related to PBC

• To understand the process of initialization in MD simulations

Keywords: Molecular dynamics, boundary conditions, simulation box, minimum image convention

We now discuss boundary conditions in more detail. The focus is on periodic boundary conditions (PBC) since they are the most used ones in molecular simulations.

Let’s recall the figure we had before and consider the three different boundary conditions in more detail.

Fig. 30 The most common boundary conditions. Left: Open boundary. The particles a free to move anywhere in space without limitations. Middle: Particles in a box with walls. Here, hard walls are assumed meaning the particles undergo elastic collisions and are reflected back from the boundary. Right: Periodic boundaries. In this case there is no physical wall, but since the system is periodic (not open), the particles have to be wrapped back into the box when they cross the border.

Periodic boundary conditions

Let’s first elaborate on the periodic boundary conditions. For simplicity, we assume that the system is two-dimensional and that the two sides are of equal length, that is, \(L_x = L_y\). Generalizations to three (or even more) dimensions and for any rectangular or parallelepiped (rhomboid) shape is straightforward.

Geometry

To be able to use PBC, the shape of the simulation box must be space filling. That a shape is space filling means that we can replicate it in different directions and as a result it fills the space. This requirement restricts the shapes of possible simulation boxes.

Fig. 31 Examples of space filling and non-space filling structures in two dimensions.

In 2-d the situation is quite intuitive (see the Figure below), but in 3-d the situation is different. In 3-d the cubic box is the most obvious one and hence a very common choice. Other possible choices are, for example, the hexagonal prism, truncated octahedron and rhombic dodecahedron. Note that a cubic (and any parallelepiped) box has 6 faces, but dodecahedron has 12 and octahedron has 14 faces, and it is possible to have space filling structures with faces up to 38. It would be tempting to think that platonic polyhedra would all be space filling, but that is not so. The cube is the only one. Shapes that can fill 3-d space are called pleisiohedra.

As it is clear from the above, it is possible to use quite complex geometries to satisfy the requirement that the structure is space filling. It is also very clear, that implementing the formulas that are needed to take care of periodicity can become quite complicated indeed. The obvious question is: Is anything else needed but the simplest structures? The answer is yes. For example, Gromacs offers cubic, rhombic dodecahedron, and truncated octahedron as possible choices. What is the reason for this? The answer is two-fold: First, using the latter two allows significant saving in the amount of solvent (water, in most cases). This results from the fact that the unit cell of a cube is much larger than that of a rhombic dodecahedron or truncated octahedron. The rhombic dodecahedron has the smallest unit cell of three, only 71% of the volume of the cubic unit cells. This fact allows a very large saving in terms of the number of water molecules and hence a large saving in computer time. The second reason is that the rhombic dodecahedron and the truncated octahedron have shape that is closer to a sphere than the cube has thus making the space more symmetric. Due to these two reasons, rhombic dodecahedron or truncated octahedron are very often used in simulations of solvated proteins and polymers.

Fig. 32 3D space filling stuctures used commonly in simulations.

Computational systems have a very limited number of particles compared to macroscopic systems. This is one of the reasons why it is very difficult and impractical to simulate systems that have real walls or boundaries. In addition, one is often interested in bulk systems, such as bulk solutions of polymers or colloids, or proteins in solution. One additional difficulty in having walls is that they tend to give rise to finite-size effects.

The most common solution is to apply Periodic Boundary Conditions (PBC). The figures below explain the concept.

PBC practicals

Periodic boundary conditions bring in several practical aspects that need to be considered. We now discuss some of them. The figure below shows a 2-d system with the primary simulation box, several image boxes with some of the particles marked. The figure on the left shows the general idea: The primary simulation box is replicated in all directions such that the space is filled: It is in this sense that we can use a small and finite simulation box to model an infinite, or bulk, system. This construction removes the need for physical walls. Note that when the system is initialized, all the particles must be placed in the primary simulation box. Let’s now focus on the right hand side of the figure.

Fig. 33 Particle motions in system with periodic boundary conditions. Left: The primary simulation box is shown in white, its nearest images on violet and the next-nearest images in purple. The system in the primary simulation box is replicated in all of the image boxes. Right: A closer look at what happens when particles leave the primary simulation box. Three particles are marked to illustrate the points.

The right hand side shows the primary simulation box with three particles marked in red, blue and gray. The arrows show the directions of their motions. Let’s first focus on the blue particle in the primary box. As the arrow shows, the particle will cross the boundary and move out of the primary simulation box. When that happens, it needs to wrapped back to the simulation box. In terms of the images this means that the blue particle in the image box to the left of the primary box will move in through the left boundary of the primary box: Since the box is periodic, the particle has been wrapped back into to the primary box. The red and gray particles shows similar processes. The question is how should this wrapping be done?

The minimum image convention: The concept

To be able to discuss the process of wrapping better, let’s focus on the figure below. The blue circle shows the neighborhood of the blue particle up to a cutoff \(r_\mathrm{cut}\). When computing the interactions between the particles, only those particles that are within \(r_\mathrm{cut}\) are considered. The gray particle that is in the image box to the right is within the circle - the gray particle that is inside the simulation box is much further away: The gray particle that is in the image box is the minimum image of that particle and hence, it is used when computing the interactions. The minimum image convention means that each particle interacts with the nearest image of any other given particle. This is specific to periodic boundary conditions.

The minimum image convenetion was first used by Metropolis et al. in their 1952 Monte Carlo simulations [10].

Before discuss how to implement the minimum image convention, we need to consider the cutoff length, \(r_\mathrm{cut}\), in more detail.

Fig. 34 The dotted blue circle shows the neighborhood of the blue particle. In MD simulations, the nighborhood is determined in terms of a cutoff length, \(r_\mathrm{cut}\). When interactions are computed, they are computed only up to \(r_\mathrm{cut}\). Here, the gray particle that is in the image box to the right of the primary simulation box: it is the image that is nearest to the blue particle in the primary box, that is, the nearest image. In a simulation, each particle interacts with the nearest image of any other given particle.

Cutoff length, \(r_\mathrm{cut}\)

From the practical perspective, PBC require the use of a cutoff length and set a condition for \(r_\mathrm{cut}\); interactions between pairs of particles are computed up to the distance of \(r_\mathrm{cut}\) from any given particle. It is immediately clear from the figure that as long as \(r_\mathrm{cut} < L/2\) (square box of length \(L\) is assumed), every pair of particles (for a given particle) is counted only once; this easy to test, simply draw a small system, apply the minimum image convention and list all the pairs. However, if \(r_\mathrm{cut} > L/2\), it is possible that a given particle interacts with a particle and the image of the same particle. That scenario implies double counting and it constitutes a serious artifact. Thus, the condition for the upper limit for the cutoff is \(r_\mathrm{cut} < L/2\).

We will later discuss how to select the cutoff length based on physical principles, but the condition

\[ r_\mathrm{cut} < L/2 \]

sets the absolute upper limit for it. Physical arguments allow to make the cutoff reasonable small which greatly improves computational efficiency, but we will return to that aspect later.

Fig. 35 The minimum distance between the gray particle is the distance from the blue particle to the image of the gray particle on the right. That distance is also within the cutoff distance, \(r_\mathrm{cut}\), shown by the dotted circle and the blue arrow.

The minimum image convention: How to use it

Let’s see how we can put the minimum image condition into practice.

Avoid wrapping particles back continuously

While it may be tempting to wrap the particles back to the primary simulation box every time they cross the boundary, it should be avoided. There are several reasons, but the main reason is this: this process introduces unnecessary numerical error, and this error can grow to be significant.

Instead, do it on the fly

Important: The coordinates are not wrapped back, we simply use the box length (\(L_x\), \(L_y\), \(L_z\)) and map the position of any given particle in primary box only for the calculation. Note also that the calculation becomes tricker for boxes that are not rectangular.

The code shows how it works, cutoff \(r_\mathrm{cut}\). You can easily check on paper and pen that the code works. There are several possible ways to implement PBC and minimum image condition.

# Loop over all pairs, ensure that each pair is taken into account only once. 
# Then, for each pair compute

x_dist = x_i - x_j 
y_dist = y_i - y_j 
z_dist = z_i - z_j 

dx     = x_dist - L_x * round (x_dist / L_x)  # round: rounds to the nearest integer
dy     = y_dist - L_y * round (y_dist / L_y)  
dz     = z_dist - L_z * round (z_dist / L_z)  

rij2   = dx*dx + dy*dy + dz+dz                # Better to use squares
if rij2 < cutoff*cutoff:
   compute interactions

We can make this even more concrete:

# Loop over all pairs, ensure that each pair is taken into account only once. 
# Then, for each pair compute

Loop i from particle 1 to particle N-1         # We assume N particles
   Loop j from particle (i+1) to particle N    # This ensures that we compute each pair just once.
      x_dist = x_i - x_j 
      y_dist = y_i - y_j 
      z_dist = z_i - z_j 

      dx     = x_dist - L_x * round (x_dist / L_x)  # round: rounds to the nearest integer
      dy     = y_dist - L_y * round (y_dist / L_y)  
      dz     = z_dist - L_z * round (z_dist / L_z)  

      rij2   = dx*dx + dy*dy + dz+dz                # Better to use squares
         if rij2 < cutoff*cutoff:
         compute interactions

There are also many other ways to implement PBC and the minimum image convention.

Possible problems

There are several possible problems that can arise from the use of boundary conditions. While the use of PBC may help to reduce the size of simulated systems, problems can arise from finite size effects. It may also occur that the boundaries have some artificial unwanted effects on the system. That is one reason for using the complex rhombic dodecahedrons or truncated octahedrons instead of simple cubic geometry for the system - those complex geometries make the space more symmetric. It is also possible, that poorly chosen, or poorly implemented, boundary conditions may give rise to articial correlations in the system. Such effects may be difficult detect and identify to be due boundaries. Boundary conditions may also have the opposite effect: small systems size and poorly chosen BCs may lead to suppression of long-range correlations. While we have not yet discussed the details of interactions in molecular systems, the choice of boundary conditions has a significant influence on how we can compute some of them, especially the electrostatic interactions. We will return to that matter when we discuss electrostatics.

Distance between images

One common problem when simulating single molecules solvated in, say, water, is that the distance between the images is too short. Such a situation leads to artificial interactions of the molecule with itself, an obvious artifact. This is easy to avoid when planning the simulations and estimating the size of the simulation box and the amount of solvent needed.

Fig. 36 One important matter to consider especially when simulating single molecules such as proteins or peptides in solution is that the distance between the images is large enough. If that is not the case, the molecule interacts with itself.

Bibliography

1

Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087, 1953. doi:10.1063/1.1699114.