Validity of Molecular Dynamics Simulations for Soft Matter

0.13 MB·English

by  Sangrak Kim

Validity of Molecular Dynamics Simulations for Soft Matter
Sangrak Kim
Department of Physics,
Kyonggi University
154-42 Gwangyosanro, Youngtong-ku,
Suwon 440-760, Korea
(Dated: January 15, 2013)
In this work, we analytically examine the validity of molecular dynamics for a soft potential
3 systembyconsideringasimpleone-dimensionalsystemwithapiecewisecontinuouslinearrepulsive
1 potentialwallhavingaconstantslopea. Wederiveanexplicitanalyticalexpressionforaninevitable
0 energy change ∆E duetothediscrete process, which is dependenton two parameters: 1) α, which
2 is a fraction of time step τ immediately after the collision with the potential wall, and 2) µ≡ aτ,
n where p0 is themomentum immediately before thecollision. The whole space of parameters α apn0d
a µcanbedividedintoaninfinitenumberofregions,whereeachregioncreatesapositiveornegative
J energy change ∆E. On the boundaries of these regions, energy does not change, i.e,∆E =0. The
envelope of |∆E| vs. µ shows a power law behavior |∆E|∝µβ, with the exponent β ≈0.95. This
3
1 implies that the round-off error in energy introduced by the discreteness is nearly proportional to
thediscrete time step τ.

]
h PACSnumbers: 02.70.Ns,02.70.Bf,02.60.-x
p Keywords: MolecularDynamics,Finite-differenceMethods,NumericalMethods
-
p
m
Hard, or completely impenetrable, core particles have If the total energy of the system is different from the
o an excluded volume effect. To simulate a hard-core sys- initial energy, the system may show very different ther-
c tem,weuseevent-drivenmethods[1]-[2]todeterminethe modynamicanddynamicbehaviors,whichnolongercor-
.
s timeatwhichthehard-coreparticlescollide. Atthetime respond to the original objective. Therefore, it is im-
c of collision, the longitudinal component of the velocities portant to know how the total energy changes with the
i
s is exchanged, but the transverse component of the ve- system control parameters in the simulations.
y
locities remains the same after the collision. Therefore, In this article, we derive an analytical expression for
h
their total energy is completely conserved, and there is energy change ∆E in a very simple system. The system
p
absolutely no energy drift. is composed of a simple particle colliding with the soft
[
In reality, atoms usually have a soft repulsive core. potential wall in one-dimension. The soft potential wall

1 This can be modeled by exponential or power functional is modeled as a piecewise linear potential, as shown in
v
forms [3]-[4]. For example, noble gases such as argons Fig. 1:
8
can be described by the Lennard-Jones potential
6
7 V (r)=4ǫ[(r)−12−(r)−6], (1) 0, if q ≥0,
2 LJ σ σ V(q)= (3)
(cid:26)aq, otherwise,
.
1 where ǫ and σ characterize the energy and length scales,
0 respectively. where q is a coordinate of the particle and a < 0 is a
3 In molecular dynamics simulations [2], we usually cal- constantthatcharacterizestheslopeofthepotential. For
1
culatethepositionsandvelocitiesoftheparticlesbysolv- simplicity, we take the mass of the particle as m=1. In
:
v ing the Newtonian equations of motion, region I, where q ≥ 0, the particle moves freely with
Xi d2r~ the Hamiltonian H(p,q) = p2, which is quite trivial. In
r mi dt2i =F~i(r~1,··· ,r~N),i=1,··· ,N, (2) region II, where q < 0, it m2oves with the Hamiltonian
a H(p,q)= p2 +aq.
for all the particles in the system. The system governed 2
The velocityVerletalgorithm,asymplectic algorithm,
by the Newtonian equations of motion has some inter-
is used to solve Eq. (2). For region I, it is written as
esting properties such as time reversibility, energy con-
servation, and so on [5]. Molecular dynamics is a digital
computingimplementationofthiscontinuousdifferential qn = q0+nτpn, (4)
equation into its corresponding discrete difference equa- (cid:26) pn = p0,
tion. Thus, keeping energy constant is a key measure in
and for region II, it is written as
assessing the validity of molecular dynamics.
An inevitable loss of accuracy is caused by represent- qn+1 = qn+τpn− 12aτ2, (5)
ing a derivative by its finite-difference approximation, (cid:26) pn+1 = pn−aτ.
termeda round-offerrorindigitalcomputing. Manydif-
ferent molecular dynamics algorithms [6]-[8] have been Atypicaldiscretetrajectoryfora=0.1,τ =0.5,start-
proposed to solve Eq. (2). ing at q0 =5.1,p0 =−1.0, is shown in Fig. 2. Note that
2
0.9
II I
0.8
0.7
0.6
0.5
q)
V( 0.4
0.3
0.2
0.1
0
−0.1
−10 −5 0 5 10
q
FIG. 1: (Color on-line) Potential model for a soft potential FIG.3: (Coloron-line)Divisionsofαandµspaceuptonc =
wall with slope a = 0.1. If a particle has an energy E0 = 10. These divisions can continue up to nc → ∞. The cross
0.5,continuousdynamicspredictsthattheparticlewillreflect pointswithabscissa aregivenbytheformulaµ= nc2−1,nc =
at q = −5.0, but in molecular dynamics, the particle may 2,··· ,∞. All points on the vertical lines and diagonal lines,
penetrate furtherinto thesoft wall. except α = 0 show ∆E = 0. Upper right-angled triangles
show ∆E >0, denoted with the+ symbol shown in thebox.
Lower right-angled triangles show ∆E <0, denoted with the
− symbol shown in the box. The numbers shown in the box
are thenumbernc.
5
4
3
2 1
n+n
1 n c 0.8
q 0
0.6
−1
0.4
−2
0.2
−3
E
∆
−4 0
−5
−0.2
0 5 10 15 20 25 30 µ=0.4
t
−0.4 µ=0.42
µ=0.44
FIG. 2: (Color on-line) A typical particle trajectory for a = −0.6 µ=0.46
0.1,τ =0.5startingat q0 =5.1with p0=−1.0. Thenumber µ=0.48
−0.8
nisthenumberofstepsimmediatelybeforethecollisionwith 0 0.2 0.4 0.6 0.8 1
α
the potential wall when the particle enters from region I to
II. The number n+nc is the number of steps immediately
afterthecollision withthepotentialwallwhentheparticleis FIG. 4: (Color on-line) Energy change ∆E vs. parameter α
leaving region II. This corresponds to α= 0.9. Note that in for differentvaluesof µ. Thiscorresponds toseeing a change
regionI,thetrajectoryisastraightline,denotedwithastar, along the vertical lines on the regions ’4 -’ and ’5 +’ in Fig.
butinregionII,itisaparabola,denotedwithafilledsphere. 3.
inregionI,thetrajectoryis astraightline,denotedwith the wall and then finally, leaves the region II. Let n de-
a star symbol, whereas in region II, the trajectory is a notethenumberofstepsimmediatelybeforethecollision
parabola, denoted with a filled sphere symbol. We now andn+n denote thenumberofstepsimmediately after
c
examine the beginning and end of the collision with the thecollision. Then,n wouldbethenumberofstepsthat
c
wall. As the particle enters a line q =0 from region I to occur during the collision. As shownin Fig. 2, when the
region II, it experiences a constant repulsive force from particle enters the potential wall, it has a different value
3
Parameter α is in the range of 0 < α ≤ 1, and µ is in
10 α=0.1
therangeof0<µ≤1sincen isanintegergreaterthan
α=0.5 c
5 α=0.9 1. From Eq. (8), for a fixed value of µ, the number of
E
∆ collisionsnc changesby1overthewholerangeofαfrom
0
0 to 1, since n (α = 1)−n (α =0)=1 for any value of
c c
−50 0.5 1 1.5 2 µ. If ξ =µ(nc−1)=2, then ∆E =0 for any value of α.
µ Thus, vertical lines µ = 2 for n = 2,··· ,∞ divide
nc−1 c
theparameterspaceintoaninfinitenumberofsubspaces,
100 asshowninFig. 3. Straightlinesα=(n −1)(1n µ−1)
c 2 c
∆ E for nc =2,··· ,∞ divide the subspaces into two regions:
upper trianglesshowapositive changein∆E, andlower
10−5 trianglesshowanegativechangein∆E. Atpointsonthe

10−4 10−3 10−2 10−1 100 boundaries, ∆E = 0. In Fig. 3, the regions are denoted
µ
with different colors and a pair of symbols including a
number n and a symbol, ’+’ or ’-’, up to n =10.
c c
FIG. 5: (Color on-line) Energy change ∆E vs. parameter
µ along the horizontal lines in Fig. 3 for different values of
parameterα. Upperpanel: ∆E vs. µ. Lowerpanel: log|∆E| Energy changes ∆E as a function of α for different
vs. logµ. The straight line shown in the lower panel has a values of µ are illustrated in Fig. 4. These correspond
slope of 0.95. to seeing a change of ∆E along the vertical lines on the
regions’4 -’,and’5+’inFig. 4. Forthe caseofµ=0.4,
every point except α=0 gives ∆E =0. For other cases,
∆E =0 occurs when crossing the diagonal lines.
of qn+1, depending on the previous position qn. We pa-
rameterize this as α, the fraction of time during which
the particle moves in region II from qn to qn+1. Thus, Fig. 5 shows the results of ∆E as a function of µ for
α = 0 corresponds to qn+1 = 0, and α = 1 corresponds different values of α. The upper panel shows ∆E vs. µ
to q = 0. The parameter α can take a value in (0, 1]. and the lower panel shows log|∆E| vs. logµ. We find
n
The case shown in Fig. 2 corresponds to α = 0.9. In that the envelope of log|∆E| is linear with a slope of
thesimulations,weactuallyencountermanydifferentre- 0.95. This exponent governsthe energy fluctuation from
alizations of α; however, unfortunately, we do not know the round-off error of a discrete time step τ.
the distribution of α a priori.
Wederiveanalyticalexpressionsforthenumberofcol- Although we havedemonstratedthe validity of molec-
lisions nc and the relative energy change during the col- ulardynamicsfromaverysimplesystemforsoftmatter,
lision, defined by theresultscanbegeneralizedtoawiderangeofsystems.
Furthermore, if we recognize many realizations of colli-
∆E ≡ pn+nc2−pn2. (6) sions, then energy fluctuations, or energy drifts in the
p 2
n system are coming from distributions of just two vari-
ables: 1) α, which is related to the phase of the collision
The number of collisions n is calculated from the in-
c
withthe potentialwallwhenthe particleentersitand2)
equality
p , which is different from collision to collision.
n
n (n −1)
(α−1)τp +n τp − c c aτ2 ≥0, (7)
n c n
2 In summary,we considereda very simple system com-
so that n can be obtained from posed of a particle entering a linear repulsive potential
c
wall with a constant slope a. By applying the velocity
1+ 1µ+ (1+ 1µ)2+2(α−1)µ Verlet algorithm, the incoming momentum pn and the
2 2
nc =Ceil( q ), (8) outgoingmomentumpn+nc couldbederived. Further,we
µ also derived the collision number n and energy change
c
whereµ≡ aτ andthefunctionCeil(x)denotesthemini- ∆E. These two values depended on only two parame-
pn ters, α and µ, with ranges of 0 < α ≤ 1 and 0 < µ ≤ 1,
mumintegerlargerthanx. Thequantityaτ corresponds
respectively. The parameter space created by α and µ
to the magnitude of momentum generated by the force
could be divided into regions by the values n and signs
field in a time step τ, and thus, parameter µ is the ratio c
of ∆E. An energy change of 0, i.e., ∆E = 0, was only
of this momentum and the momentum of the incoming
seen on the boundaries of the regions. The envelope of
particle. The relativeenergychange∆E is calculatedby
log|∆E| followed a power law with an exponent of 0.95.
∆E =ξ(ξ−2), (9) Roughly speaking, |∆E|≈µ0.95. Finally, any two parti-
cles with a soft repulsive potential core could be treated
whereξ ≡µ(n −1). Notethatn and∆E arefunctions inthe same wayif we introducedthe reducedcoordinate
c c
of only two dimensionless parameters, α and µ. system. This needs to be examined further.
4
Acknowledgement Korea (NRF) funded by the Ministry of Education Sci-
ence and Technology (2012-0002969).
ThisresearchwassupportedbyBasicScienceResearch
Program through the National Research Foundation of
[1] B. Alder and T. Wainwright, J. Chem. Phys. 27, [5] H. Goldstein,C. PooleandJ. Safko,ClassicalMechanics
1208(1957). (Perarson, New York, 2002); J. Lowerstein, Essentials
[2] J. Haile, Molecular Dynamics Simulation (John Wiley & of Hamiltonian Dynamics (Cambridge University Press,
Sons, New York,1992). Cambridge, 2012).
[3] D. Frenkel and B. Smit, Understanding Molecular Sim- [6] L. Verlet,Phys. Rev.159, 98(1967).
ulation (Academic Press, New York,2002). [7] D. Beeman, J. Comp. Phys.130, 98(1976).
[4] M. Griebel, S. Knapek and G. Zumbusch, Numeri- [8] M. Allenand D. Tildesley, Computer Simulation of Liq-
cal Simulation in Molecular Dynamics (Springer Verlag, uids (Oxford University Press, Oxford, 1987).
Berlin, 2007).