278 8 Differential equations
Å (solid argon),σ= 3. 405 Å,ε= 1. 0318 · 10 −^2 eV. Another common practice is puttingE 0 = 4 ε,
affecting the conversion factorsF 0 ,T 0 andt 0.
Normally distributed random numbers are obtained by performing a Box-Muller transform
on uniformly distributed numbers. Letuandvbe uniform numbers in the interval(− 1 , 1 ).
These numbers will only be accepted for the transformation ifs=u^2 +v^2 is in the interval
( 0 , 1 ). In that case, we obtain two normally distributed numbersn 1 andn 2 by multiplyingu
andvwith a constant,
n 1 =Su, n 2 =Sv, S=√
−lns
s. (8.24)
n 1 andn 2 will have standard deviations of 1, but multiplying all generated numbers with a
constant will give a distribution with that constant as the standard deviation.
a) Write a program that generates anNc×Nc×Ncunit cell face centered cubic lattice of argon
atoms. If you use an object oriented programming language, each atom and/or the entire
lattice should be objects of a class.
For easy testing of the lattice arranger and the later MD implementation, you should al-
ready visualize your atoms. VMD is a visualization program with a simple output format
and pretty graphics. It can be downloaded and run from your home area, and a description
of its output format can be found in the appendix.
b) Consider first free particles with initial independent Maxwell-Boltzmann distributed veloc-
ities. These correspond to normally distributed values with standard deviation
√
kBT/mfor
the desired temperatureT. In a system ofNatoms, all 3 Nvelocity componentsvare set
using
v=√
kBT/mξ (8.25)
whereξis a normally distributed number with mean 0 and standard deviation 1. Remove
any initial total linear momentum from the system.
Integrate the dynamical equation (N2L) using the symplectic and numerically stable veloc-
ity Verlet algorithm. For each particlei, the steps are as follows (currently settingUi= 0 ):vi(t+∆t/ 2 ) =vi(t)+
Fi(t)
2 m
∆t (8.26)
ri(t+∆t) =ri(t)+vi(t+∆t/ 2 )∆t (8.27)
Fi(t+∆t) =−∇iUi({r}(t+∆t)) (8.28)vi(t+∆t) =vi(t+∆t/ 2 )+
Fi(t+∆t)
2 m
∆t (8.29)The particles will now spread out into space. We are only interested in bulk atoms in a
material, so the next step is implementing periodic boundary conditions. Every time the
position of a particle is updated, the program must check if it has gone though one of the
sides.
c) Create a function for calculating the force between all particles. Use the Lennard-Jones
potential, which has the following form:
Ui j(ri j) = 4 ε[(
σ
ri j) 12
−
(
σ
ri j) 6 ]
(8.30)whereri j=ri−rj. Differentiate the expression analytically for finding theforce. Summing
up the potential energy between all particles can also be useful.
Use the minimum image convention when calculating the distance between particles. E.g.
the distance between atoms/atom replicaeiandjin thexdirection becomesminδ(xi−xj+δL)