Molecular Dynamics (1)
How to start up?
Initial configuration:
Criteria for choosing initial configuration,quick relaxation to the
equilibrium state.
- Random initial configuration:
Advantage,liquid-like,easy to generate.
Drawback,substantial overlaps at high densities.
Use,Monte-Carlo for systems without hard-core interactions,Not
suitable for hard-core interactions and MD at high densities,
- Crystal initial configuration:
e.g.,fcc crystal
Advantage,avoiding overlaps,easy to generate.
Drawbacks,not liquid-like
Use,suitable for all type of MC and MD.
Initial velocities:
Maxwell-Boltzmann distribution:
)2ex p (2)(
2
kTkTf vmmv ixiiix
Uniform distribution:
Each velocity components drawn from a uniform distribution in
a range (-vmax,vmax).
Constraint:
0
1

N
i iivm
P
Exercise
1) Write a program for generating an initial configuration with
the particles places on a fcc lattice.
2) Write a program for generating initial velocities from an
uniform distribution in the range (-vmax,vmax) (vmax=5) with the
constraint that the total initial velocity is zero.
Equilibration
Monitoring the equilibration of system with several observables,
e.g.,internal energy,pressure etc..
A
t
The production run should not start before the system reaches
the targeted equilibrium state!
Molecular Dynamics for Hard Spheres
Read initial configuration and initial velocities
Identify the next collision
Move all particles forward until collision occurs
Calculate the velocity change of the colliding pair
Calculate observables or write out the current configuration
End of program
loop
Identify the next collision
Collision condition,bij = rij.vij < 0
i jrij
vij
going awayi jrij
vij
approaching
Collision time:
| rij + vij tij | = s => vijtij2 + 2bijtij +rij2 - s2 = 0
If bij2 - vij2(rij2 - s2)? 0 (condition for collision taking
place effectively),the collision time is given by
tij = {- bij - [bij2 - vij2(rij2 - s2)]1/2}/vij2
Collision-Time List and Partner list
Initializing CT list - COLTIM(N):
From the initial condition,
Calculate for each particle its collision time with all the other
particles,
Find out the minimum of the collision times of particle I and put
the minimum collision time into COLTIM(I) (if particle I will
not collide with any other particle,put a very large value,e.g.,
TIMBIG,into COLTIM(I)).
Repeat this for all the particles.
Partner list - PARTNR(N):
PARTNR(I) contains the identity (number of the particle) of the
collision partner of particle I.
Move all particles forward until collision occurs
From all the collision times stored in COLTIM,Find out the minimum
one,(tij)min.
Advance the position of all the particles with this time interval,i.e.,
do i=1,N
RX(i)=RX(i)+VX(i)*(tij)min
RY(i)=RY(i)+VY(i)*(tij)min
RZ(i)=RZ(i)+VZ(i)*(tij)min
enddo
Calculate the velocity change of the colliding pair
Dynamics of elastic collision:
i j
rij vij(before collision)
vij(after collision)
Energy conservation requires,
|vij(after collision)| = |vij(before collision)|
Since the momentum change takes place along rij and no momentum
change in the direction perpendicular to rij,the collision is specular
with respect to rij,This leads to
vi(after) = vi(before) + dv
vi(after) = vi(before) + dv
dv = [-(bij/s2) rij]collision
Update COLTIM and PARTNR
There is no need to recalculate totally COLTIM and PARTNR!
It is obvious that it is necessary to update for the two particles
involved in the collision (say,I and J),i.e.,COLTIM(I),
COLTIM(J),PARTNR(I) and PARTNR(J).
Moreover,we need also to update for all the particles which have
I or J as collision partner!
Now,one can repeat all the previous steps to form a loop.