## 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.

## Leave a Reply