Molecular Dynamics: Diffusion with PBC

How can I implement the computation of the diffusion coefficient $D$ using periodic boundary conditions (PBC)?

I use molecular dynamics of a set of $nboby$ particles with positions $pos(3,nbody)$ in a box of length $length$. The implemetion of the PBC is

do k = 1, 3
pos(k,i) = modulo(pos(k,i),length)
end do

At now I'm using for $D$ the following code

do it=1,nstep
diff=0
do i = 1,nbody
pos2(:)=pos2(:)+pos(:,i)
diff=diff+dot_product(pos(:,i),pos(:,i))
end do
diff=diff/nbody-dot_product(pos2(:),pos2(:))/nbody**2
end do
diff=diff/nstep/6

which I think it corresponds to

$D=\lim_{t\to\infty}\dfrac{<x^2>-<x>^2}{6t}$

but I'm not very sure that the PBC are taking into account in the right way.

Can someone help me?

Thanks Matteo

You need to unfold the positions after the simulation is done. I use the following subroutine:

subroutine unfold_positions(L, X, Xu)
! Unwinds periodic positions. It is using the positions of particles from the
! first time step as a reference and then tracks them as they evolve possibly
! outside of this box. The X and Xu arrays are of the type X(:, j, i), which
! are the (x,y,z) coordinates of the j-th particle in the i-th time step.
real(dp), intent(in) :: L ! Box length
! Positions in [0, L]^3 with possible jumps:
real(dp), intent(in) :: X(:, :, :)
! Unwinded positions in (-oo, oo)^3 with no jumps (continuous):
real(dp), intent(out) :: Xu(:, :, :)
real(dp) :: d(3), Xj(3)
integer :: i, j
Xu(:, :, 1) = X(:, :, 1)
do i = 2, size(X, 3)
do j = 1, size(X, 2)
Xj = X(:, j, i-1) - X(:, j, i) + [L/2, L/2, L/2]
Xj = Xj - L*floor(Xj/L)
d = [L/2, L/2, L/2] - Xj
Xu(:, j, i) = Xu(:, j, i-1) + d
end do
end do
end subroutine

I think you can replace the line Xj = Xj - L*floor(Xj/L) with Xj = modulo(Xj, L), good point, I didn't realize that you can use modulo for floating point numbers.

• Can you please elucidate/elaborate upon the logic of this code ? It would be of immense help !! Sep 19 '15 at 9:28
• I tried to use your line of logic, it doesn't seem to work for me :( Could be that I was wrong, but can you please tell me the idea behind this code ??? Is this code with the condition that the origin is at the centre of the cubic box ?? Sep 20 '15 at 11:59
• The problem is in some cases, I am getting a difference of L, between folded and unfolded coordinates !! Sep 20 '15 at 12:35
• As documented in the subroutine, the box is assumed to be [0, L]^3. The logic of the code is that it checks if an atom leaves the box and appears on the other side, and instead it keeps the trajectory continuous, i.e. it unfolds the trajectories. If it is still not clear, then please ask a specific question, preferably with code that doesn't work. Sep 21 '15 at 2:10

I solve my own problem. As I said in the question I'm using PBC:

do k = 1, 3
pos(k,i) = modulo(pos(k,i),length)
end do

and there is a dynamics for the positions $pos(k,i)$ of the particles

pos(k,i)=pos(k,i)+vel(k,i)*dt+0.5*force(k,i)*dt**2

For the computation of the diffusion coefficient I used a second variable of position, let say $pos0(k,i)$. This new variable evolve like the original one:

pos0(k,i)=pos0(k,i)+vel(k,i)*dt+0.5*force(k,i)*dt**2

but to them I don't apply the PBC. Therefore the computation of the diffusion coefficient is something like:

do it=1,nstep
diff=0
do i = 1,nbody
pos2(:)=pos2(:)+pos0(:,i)
diff=diff+dot_product(pos0(:,i),pos0(:,i))
end do
diff=diff/nbody-dot_product(pos2(:),pos2(:))/nbody**2
end do
diff=diff/nstep/6

Other quantities instead are computed using the original variables $pos(k,i)$, for which the PBC holds:

do k = 1, 3
pos(k,i) = modulo(pos(k,i),length)
end do

I hope that this answer can be helpful to who is looking for something similar.