This Blog is Untitled

Periodic boundary conditions (PBC) are ubiquitous in molecular simulations. However, when visualizing results, one often prefers to see whole molecules rather than fragments cut by the simulation box. If the molecule percolates through the box, minimizing artificial cuts becomes essential. Below are two methods for reconstructing intact structures under PBC:

1. For non‐percolated molecules: apply a differential–integral approach based on the idea that any bond vector between adjacent particles never exceeds half the box length. For a simple linear polymer, you “unwrap” each bond vector (the differential step) and then cumulatively sum them (the integral step), placing the first bead at the origin:
   

   # for x -> (n_particles, n_dim) and box -> (n_dim,)
   x_p = np.cumsum(pbc(np.diff(x, prepend=x[:1], axis=0), box), axis=0)

   
   For more complex topologies, the principle is the same; instead of using np.diff and np.cumsum, you perform the differential and integral operations along a BFS or DFS traversal of the molecule’s bonding graph.
2. For percolated networks: the differential–integral method fails once the molecule spans the entire box. Instead, you can compute a circular mean, which is designed to find the average of periodic data. (A familiar example is averaging sleep times that wrap around midnight.) First, calculate the circular mean of each coordinate to estimate the periodic center of mass, then shift and unwrap:
   

   from scipy.stats import circmean
   x_com = np.asarray([circmean(x.T[d], high=box.T[d], low=box.T[d]) for d in range(x.shape[1])])
   x_p = pbc(x - x_com, box)

   
   The circular mean works by mapping each data point $s_n$ to an angle $\alpha_n ∈ [0, 2π)$, projecting it onto the unit circle $z_n = e^{\mathrm{i}\alpha_n}$, and then taking the argument of the mean of all $z_n$: $$\arctan2(\sum \sin\alpha_n, \sum \cos\alpha_n)$$ This approach also applies to non‐percolated molecules. However, if the molecule is very elongated (for example, pipette-bulb‐shaped), some atoms may lie more than half a box length from the circular mean. In that case, it’s better to use a circular midpoint (an unweighted center). One simple way to compute it is by histogramming the angles and taking the mean over occupied bins:
   

   bins = 500
   bin_size = 2 * np.pi / bins
   d = np.linspace(0, 2 * np.pi, bins, endpoint=False)
   p = np.bincount((x / bin_size).astype(np.int), minlength=bins) # x is scaled into [0, 2π)
   x_mid = circmean(d[p > 0])

In summary, for simple polymers and other non‐percolated molecules, the differential–integral method is the most straightforward. For fully percolated networks and self‐assembled structures, using the circular mean is ideal. The circular midpoint method works in all cases but involves extra steps and one hyperparameter (the histogram bin count).