Based on the attached pseudo codes, write a full MD code.

Argon parameters can be used for the LJ parameters in your code.

Calculate Radial distribution function from the output of your code.

Algorithm $5($Calculation of the Forces$)$

subroutine force$($f$,$en$)$

en$=0$

do i$=1,$npart

f$($i$)=0$

enddo

do i$=1,$npart$-1$

do j$=$i$+1,$npart

xr$\}=\backslash$textrm$\{$x$\}($i$)-\backslash$textrm$\{$x$\}($j

xr$=$xr$-$box$*$nint $($xr$/$box$)$

r$2=$xr$***2$

if $($r$2.$lt$.$rc$2)$ then

r$2$i$=1/$r$2$

r$6$i$=$r$2$i$**3$

ff$=48$r$6$i$*($r$6$i$-0\mathrm{.}5)$

f$($i$)=$f$($i$)+$ff$*$xr

f$($j$)=$f$($j$)-$ff$*$xr

en$=$en$+4($r$6$i$-1)-$ecut

endif

enddo

enddo

return

end

determine the force

and energy

set forces to zero

loop over all pairs

periodic boundary conditions

test cutoff

Lennard$-$Jones potential

update force

update energy

Comments to this algorithm:

For efficiency reasons the factors $4$ and $48$ are usually taken out of the force

loop and taken into account at the end of the calculation for the energy.

The term ecut is the value of the potential at $r=r_c$; for the Lennard$-$Jones

potential, we have

ecut $=4(\frac{1}{r_c^{12}}-\frac{1}{r_c^6})$

Algorithm $4($Initialization of a Molecular Dynamics Program$)$

subroutine init

sumv$=0$

sumv$2=0$

do i$=1,$npart

x$($i$)=$lattice$\_$pos$($i$)$

v$($i$)=(\backslash$operatorname$\{$ranf$\}()-0\mathrm{.}5)$

sumv$=$sumv$+$v$($i$)$

sumv$2=$sumv$2+$v$($i$)**2$

enddo

sumv$=$sumv$/$npart

sumv$2=$sumv$2/$npart

fs$=$sqrt $(3*$temp$/$sumv$2)$

do i$=1,$npart

v$($i$)=($v$($i$)-$sumv$)*$fs

xm$($i$)=$x$($i$)-$v$($i$)*$dt

enddo

return

end

initialization of MD program

place the particles on a lattice

give random velocities

velocity center of mass

kinetic energy

velocity center of mass

mean$-$squared velocity

scale factor of the velocities

set desired kinetic energy and set

velocity center of mass to zero

position previous time step

Comments to this algorithm:

Function lattice$\_$pos gives the coordinates of lattice position i and

ranf$()$ gives a uniformly distributed random number. We do not use a

Maxwell$-$Boltzmann distribution for the velocities; on equilibration it will be$-$

come a Maxwell$-$Boltzmann distribution.

In computing the number of degrees of freedom, we assume a three$-$di$-$

mensio

Algorithm $6($Integrating the Equations of Motion$)$

subroutine integrate$($f$,$en$)$

sumv$=0$

sumv$2=0$

do i$=1,$ npart

xx$=2$f$($i$)$

vi$=($xx$-$xm$($i$))/(2*$delt$)$

sumv$=$sumv$+$vi

sumv$2=$sumv$2+$vi$**2$

xm$($i$)=$x$($i$)$

x$($i$)=$xx

enddo

temp$=$sumv$2/(3*$npart$)$

etot $=($en$+0\mathrm{.}5*$sumv$2)/$npart

return

end

integrate equations of motion

MD loop

Verlet algorithm $(4\mathrm{.}2\mathrm{.}3)$

velocity $(4\mathrm{.}2\mathrm{.}4)$

velocity center of mass

total kinetic energy

update positions previous time

update positions current time

instantaneous temperature

total energy per particle

Comments to this algorithm:

The total energy etot should remain approximately constant during the sim$-$

ulation. A drift of this quantity may signal programming errors. It therefore

is important to monitor this quantity. Similarly, the velocity of the center of

mass sumv should remain zero.

In this subroutine we use the Verlet algorithm $(4\mathrm{.}2\mathrm{.}3)$ to integrate the equa$-$

tions of motion. The velocities are calculated using equation $(4\mathrm{.}2\mathrm{.}4).$nal

Algorithm $7($The Radial Distribution Function$)$

subroutine gr$($switch$)$

if $($switch$.$eq$\mathrm{.}0)$ then

ngr$=0$

delg$=$box$/(2*$nhis$)$ bin size

do i$=0,$nhis

g$($i$)=0$

enddo

else if $($switch$.$eq$\mathrm{.}1)$ then

ngr$=$ngr$+1$

do i$=1,$ npart$-1$

do j$=$i$+1,$npart

xr$=$x$($i$)-$x$($j$)$

xr$=$xr$-$box$*$nint $($xr$/$box$)$

r$=$sqrt$($xr$**2)$

if $($r$.$lt$.$box$/2)$ then

ig$=$int $($r$/$delg$)$

g$($ig$)=$g$($ig$)+2$

endif

enddo

enddo

else if $($switch$.$eq$\mathrm{.}2)$ then

do i$=1,$nhis

r$=$delg$*($i$+0\mathrm{.}5)$

vb$=(($i$+1)3)*$delg$**3$

nid$=(4/3)$vb$*$rho

g$($i$)=$g$($i$)/($ngrnid$)$

enddo

endif

return

end

radial distribution function

switch $=0$ initialization,

$=1$ sample, and $=2$ results

initialization

nhis total number of bins

sample

loop over all pairs

periodic boundary conditions

only within half the box length

contribution for particle i and $j$

determine $g(r)$

distance $r$

volume between bin $i+1$ and $i$

number of ideal gas part. in vb

normalize $g(r)$

Comments to this algorithm:

For efficiency reasons the sampling part of this algorithm is usually combined

with the force calculation $($for example, Algorithm $5).$

The factor $=3\mathrm{.}14159$dots. system $($in fact, we approximate $N_f$ by $3$ N $).$