Molecular Dynamics: Diffusion with PBC

Question

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

Answer 1

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.

Answer 2

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.