While calculating various quantities that involves particle positions ( for example, radial distribution function, mean square displacement etc.) using Molecular Dynamics simulation, the effects of Periodic Boundary Conditions ( PBC ) are removed by adopting some unwrapping procedure. I don't understand how it is done. Can anybody explain the concept ?
When you are doing an MD simulation, you should ideally be using a large simulation box so that you can study the properties of the system. However, this is easier said than done because of the computational resources that would be required to study the dynamics of very large systems. So we instead take simulation boxes of trackable size, and assume periodic boundary conditions so that you are essentially simulating a very large system.
Now to the question of unwrapping the trajectory to remove the PBC. When you are doing a simulation using PBC, when a particle leaves the simulation box at one of the ends, an equivalent particle enters the box from the other end. What unwrapping does is essentially reversing this transformation, so that the PBC image is shifted back to where it should have been, if no PBC were applied.
Let us consider a collection of helium atoms in a box of dimensions 16 x 16 x 16 cubic angstroms. Now, in the simulation if one of the He atom leaves the box, a corresponding periodic image will replace it in the box. This means that the He atom was folded into the simulation box. The transformation works in this way:
check if the x-coordinate of the atom is greater than positive x-end: if yes set x-coordinate as (x-coordinate - box length in x) check if the y-coordinate of the atom is greater than positive y-end: if yes set y-coordinate as (y-coordinate - box length in y) check if the z-coordinate of the atom is greater than positive z-end: if yes set z-coordinate as (z-coordinate - box length in z) check if the x-coordinate of the atom is less than negative x-end: if yes set x-coordinate as (x-coordinate + box length in x) check if the y-coordinate of the atom is less than negative y-end: if yes set y-coordinate as (y-coordinate + box length in y) check if the z-coordinate of the atom is less than negative z-end: if yes set z-coordinate as (z-coordinate + box length in z)
This is how PBC wrapping works. So when you are doing a PBC unwrapping you negate these transformations.
Just to clarify a few points: You generally unwrap a trajectory if you have bonds that span across the simulation box. In such cases, you add or subtract the corresponding cell vectors to move the wrapped atom to the position where it should have been, if not wrapped. Consider the example of a water molecule, with one of the hydrogen atom getting wrapped to the negative x-end of the box, while the other hydrogen and oxygen atoms are at the positive x-end of the box. Now when you unwrap, the hydrogen atom that got wrapped will be shifted to the original position, which was (positive x-end + wrapped x-coordinate).
Please read the transcript starting from here for a detailed walkthrough on how is PBC unwrapping done: