Introduction

This is a molecular dynamics simulator I programmed in Fortran in order to understand MD in more detail. Periodic boundary conditions, energy minimizer, thermostat and barostat, as well as a correct physical scale and forcefields for the noble gasses are already implemented and working well. Future steps might include parallelization of the calculation of particle interactions and maybe introducing electrostatics and more complex molecules like water.

Github

Theory

The code can be found on Github in the link above. I tried to comment everything as good as possibile to the extent, that I think most of it can be easily understood. I will still explain the main integration loop and features a bit in the following text.

Prerequisites

Most of the algorithm is implemented in the main file. The constants file contains all necessary constants and parameters needed. The forcefield file contains four different force returning functions.

Iteration Loop

! Integration loop
integrator: do iter = 0, maxiter

The first thing to be done each iteration is fulfilling the periodic boundary conditions (pbc). The first step is checking for particles, which are not in the boundaries. If they are not, they will be “teleported” across the boundary by adding or subtracting the current box size in that dimension.

  ! "Particle-teleporter" to fulfill periodic boundary condition
  do particle = 1, nparticles ! Iterate over all particles
    do dim = 1, 3 ! Iterate over 3d space

      if (r(dim, particle, 0) > boxsize(dim)) then
        r(dim, particle, :) = r(dim, particle, :) - boxsize(dim)
      else if (r(dim, particle, 0) < 0.0_dp) then
        r(dim, particle, :) = r(dim, particle, :) + boxsize(dim)
      end if

    end do
  end do

The next step is calculating the interactions of the particles with each other. For this, we iterate over each particle, which iterates over all other particles expect itself again. For each interaction, we need to calculate the real-space distance in all three coordinates called d0. Now, to fulfill the pbc, we need to calculate interactions across the boundaries. To do that, we iterate over all three dimensions, each from -1 to +1. This means, for the calculation, not only the real-space cell is used, but also all 26 nearest neighbors (containing the image particles). For each of the image particle interaction, the distance is modified by the respective additional distance given by our current boxsize.

After the distance in thee coordinates is calculated, the total distance l is calculated. This is needed for the cutoff decision as well as the calculation of the force, which is then scaled by the distance in the three spatial coordinates. The force on each particle which is exert by the other particles is added up.

  ! Calculation of particle interactions
  f = 0
  do particle1 = 1, nparticles   ! Interaction for each particle
    do particle2 = 1, nparticles ! ..with each other particle
      if (particle1 /= particle2) then  

        d0 = r(:, particle2, 0) - r(:, particle1, 0)
        do dim1 = -1,1
          do dim2 = -1,1
            do dim3 = -1,1

              ! Define new position of image and for dim1=dim2=dim3=1 of real particle interaction
              d = d0
              d(1) = d(1) + (dim1) * boxsize(1)
              d(2) = d(2) + (dim2) * boxsize(2)
              d(3) = d(3) + (dim3) * boxsize(3)

              l = SQRT(d(1) ** 2 + d(2) ** 2 + d(3) ** 2)

              if (l < cutoff) then
                f(:,particle1,1) = f(:,particle1,1) + fitted_lj_f(d, l, ff_parameters(:,atom_type))
              end if

            end do
          end do
        end do

      end if
    end do
  end do

After that, the integration step can be made using the Verlet algorithm. The velocity is also calculated.

! Do Verlet integration step
  r(:,:,1) = 2 * r(:,:,0) - r(:,:,-1) + (f(:,:,1) / mass(atom_type)) * (dt ** 2)

! Calculate velocities by numerical derivation
  v = (r(:,:,1) - r(:,:,-1)) / (2 * dt) ! v(t+dt)

After the integration step, the whole positional array is shifted upwards and thus making it ready for the next iteration. This is needed, because r in this case is used in the following way: r(xyz-dimension,particlenumber,[t-dt,t,t+dt,t+2dt].

! Move array to t = t+dt for next integration step
  r = cshift(r,  1, DIM=3)

Now, the thermostat follows. This is easily done by implementing the Berendsen weak coupling thermostat.

! Thermostat
  ! Calculate E and instantaneous temperature
  E = (mass(atom_type) / 2) * sum(v ** 2)
  T_inst = ((2 * E) / (3 * kB * nparticles))
  ! Weak coupling thermostat
  v = v * SQRT(1 + (dt / tau_T) * (T0 / T_inst - 1))
  ! Apply velocity to r(t-dt) for algorithm
  r(:,:,-1) = r(:,:,0) - (v(:,:) * dt)

Additionally, the Berendsen barostat is implemented.

! Barostat
  ! Calculate instantaneous volume and pressure
  V_inst = boxsize(1) * boxsize(2) * boxsize(3)
  P_inst = ((2 * E) / (3 * V_inst)) + (1 / (3 * V_inst)) * (sum(r(:,:,-1) * f(:,:,1)))
  ! Scaling factor
  mu = (1 - (dt / tau_P) * (P0 - P_inst)) ** (1./3.)
  ! Scale coordinates and boxsize by mu
  r = r * mu
  boxsize = boxsize * mu
end do integrator

With this, the elementary steps of the whole simulation are explained.