{"id":371,"date":"2025-05-20T13:00:54","date_gmt":"2025-05-20T05:00:54","guid":{"rendered":"https:\/\/www.shirui.me\/blog\/?p=371"},"modified":"2025-05-21T01:47:15","modified_gmt":"2025-05-20T17:47:15","slug":"unwrap-big-molecule-in-pbc","status":"publish","type":"post","link":"https:\/\/www.shirui.me\/blog\/2025\/05\/20\/unwrap-big-molecule-in-pbc\/","title":{"rendered":"Unwrap big molecules & self-assemblies in PBC"},"content":{"rendered":"\n

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:<\/p>\n\n\n\n\n\n\n\n

    \n
  1. For non\u2010percolated molecules: apply a differential\u2013integral approach based on the idea that any bond vector between adjacent particles never exceeds half the box length. For a simple linear polymer, you \u201cunwrap\u201d each bond vector (the differential step) and then cumulatively sum them (the integral step), placing the first bead at the origin:
    \n
    # for x -> (n_particles, n_dim) and box -> (n_dim,)\nx_p = np.cumsum(pbc(np.diff(x, prepend=x[:1], axis=0), box), axis=0)<\/code><\/pre>
    \nFor more complex topologies, the principle is the same; instead of using np.diff<\/code> and np.cumsum<\/code>, you perform the differential and integral operations along a BFS or DFS traversal of the molecule\u2019s bonding graph.<\/li>\n\n\n\n
  2. For percolated networks: the differential\u2013integral 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\nx_com = np.asarray([circmean(x.T[d], high=box.T[d], low=box.T[d]) for d in range(x.shape[1])])\nx_p = pbc(x - x_com, box)<\/code><\/pre>
     The circular mean works by mapping each data point $s_n$ to an angle $\\alpha_n \u2208 [0, 2\u03c0)$, 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\u2010percolated molecules. However, if the molecule is very elongated (for example, pipette-bulb\u2010shaped), some atoms may lie more than half a box length from the circular mean. In that case, it\u2019s 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\nbin_size = 2 * np.pi \/ bins\nd = np.linspace(0, 2 * np.pi, bins, endpoint=False)\np = np.bincount((x \/ bin_size).astype(np.int), minlength=bins) # x is scaled into [0, 2\u03c0)\nx_mid = circmean(d[p > 0])<\/code><\/pre><\/li>\n<\/ol>\n\n\n\n
    In summary, for simple polymers and other non\u2010percolated molecules, the differential\u2013integral method is the most straightforward. For fully percolated networks and self\u2010assembled 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).<\/p>\n\n\n\n
    <\/p>\n","protected":false},"excerpt":{"rendered":"
    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:<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[3,9],"tags":[16],"class_list":["post-371","post","type-post","status-publish","format-standard","hentry","category-misc","category-notes","tag-molecular-dynamics"],"_links":{"self":[{"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/posts\/371","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/comments?post=371"}],"version-history":[{"count":8,"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/posts\/371\/revisions"}],"predecessor-version":[{"id":395,"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/posts\/371\/revisions\/395"}],"wp:attachment":[{"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/media?parent=371"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/categories?post=371"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.shirui.me\/blog\/wp-json\/wp\/v2\/tags?post=371"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}