Determination of Forces from a Potential in Molecular Dynamics

0.54 MB·

by  Bernard Monasse

Determination of Forces from a Potential
4
in Molecular Dynamics
1
0
(note)
2

n
a
J Bernard Monasse Frédéric Boussinot

6 Mines-ParisTech, Cemef Mines-ParisTech, Cemef

] [email protected] [email protected]
h
c January 2014
e
m
-
t Abstract
a
t InMolecularDynamics(MD),theforcesappliedtoatomsderivefrom
s
. potentialswhichdescribetheenergyofbonds,valenceangles,torsionan-
t
a gles,andLennard-Jonesinteractionsofwhichmoleculesaremade. These
m definitions are classic; on the contrary, their implementation in a MD
- system which respects the local equilibrium of mechanical conditions is
d usually not described. The precise derivation of the forces from the po-
n
tentialandtheproofthattheirapplicationpreservesenergyistheobject
o
ofthisnote. Thisworkispartofthebuildingofamulti-scaleMDsystem,
c
presently under development.
[

1
Keywords. Molecular Dynamics ; Force-Fields ; Potentials ; Forces.
v
1
8
1 Introduction
1
1
. Numericalsimulationatatomicscalepredictssystemstatesandpropertiesfrom
1
a limited number of physical principles, using a numerical resolution method
0
4 implemented with computers. In Molecular Dynamics (MD) [2] systems are
1 organicmolecules,metallicatoms,orions. Weconcentrateonorganicmolecules,
: but our approach could as well apply to other kinds of systems. The goal is to
v
i determine the temporal evolution of the geometry and energy of atoms.
X
AtthebasisofMDistheclassical(newtonian)physics,withthefundamental
r equation:
a
→− →−
F =ma (1)
→− →−
where F istheforceappliedtoaparticleofmassmand a isitsacceleration
(second derivative of the variation of the position, according to time).
1
A force-field is composed of several components, called potentials (of bonds,
valenceangles,dihedralangles,vanderWaalscontributions,electrostaticcontri-
butions,etc.) andisdefinedbytheanalyticalformofeachofthesecomponents,
andbytheparameterscaracterizingthem. Thebasiccomponentsusedtomodel
molecules are the following:
• atoms, with 6 degrees of freedom (position and velocity);
• bonds, which link two atoms belonging to the same molecule; a bond
between two atoms a,b tends to maintain constant the distance ab.
• valenceangles, whicharetheangleformedbytwoadjacentbondsbaetbc
in a same molecule; a valence angle tends to maintain constant the angle
a(cid:99)bc. A valence angle is thus concerned by the positions of three atoms.
• torsionangles(alsocalleddihedralangles)aredefinedbyfouratomsa,b,c,d
consecutivelylinkedinthesamemolecule: aislinkedtob,btoc,andcto
d; a torsion angle tends to priviledge particular angles between the planes
abc and bcd. These particular angles are the equilibrium positions of the
torsionpotential(minimalenergies). Inmostcases, theyareTrans(angle
of 180◦), Gauche (60◦) or Gauche’ (−60◦).
• van der Waals interactions apply between two atoms which either belong
to two different molecules, or are not linked by a chain of less than three
(or sometimes, four) bonds, if they belong to the same molecule. They
are pair potentials.
All these potentials depend on the nature of the concerned atoms and are
parametrized differently in specific force-fields. Molecular models can also con-
sider electrostatic interactions (Coulomb’s law) which are pair potentials, as
van der Waals potentials are; their implementation is close to van der Waals
potentials, with a different dependence to distance.
Intra-molecularforces(bonds,valenceangles,torsionangles)aswellasinter-
molecularforces(vanderWaals)areconservative: theworkbetweentwopoints
does not depend on the path followed by the force between these two points.
Thus, forces can be defined as derivatives of scalar fields. From now on, we
consider that potentials are scalar fields and we have:
→− →−
F(r)=−∇U(r) (2)
→−
where r denotes the coordinates of the point on which the force F(r) applies,
and U is the potential from which the force derives.
TheworkpresentedhereispartofaMDsystempresentlyunderdevelopment[1];
the defined forces and their implementations have been tested with it.
2
Structure of the text
Bonds are considered in Section 2, valence angles in Section 3, torsion angles
in Section 4, and finally Lennard-Jones potentials in Section 5. Section 6 sum-
marises the force definitions, and Section 7 concludes the text. A summary of
the notations used in the paper is given in the Annex.
2 Bonds
Abondmodelsasharingofelectronsbetweentwoatomswhichproducesaforce
betweenthem. Thisforceisthederivativeofthebondpotentialdefinedbetween
thetwoatoms. Fig. 2.1showsa(attractive)forceproducedbetweentwolinked
atoms a and b.
Figure 2.1: Attractive bond between two linked atoms
A harmonic bond potential is a scalar field U which defines the potential
energy of two atoms placed at distance r as:
U(r)=k(r−r )2 (3)
0
where k is the strength of the bond and r is the equilibrium distance (the
0
distance at which the force between the two atoms is null). We thus have:
∂U(r)
=2k(r−r ) (4)
∂r 0
The partial derivative of U according to the position r of a is:
a
∂U(r) ∂U(r) ∂r
= . . (5)
∂r ∂r ∂r
a a
But:
∂r
=1 (6)
∂r
a
We thus have:
∂U(r)
=2k(r−r ) (7)
∂r 0
a
3
→− →−
Letaandbbetwoatoms,and u =norm(ba)bethenormalizationofvector
→−
ba. The force produced on atom a is:
→− ∂U(r) →− →−
f =− .u =−2k(r−r ).u (8)
a ∂r 0
a
and the one on b is the opposite, according to the action/reaction principle:
→− →−
f =−f (9)
b a
→−
Therefore, if r >r , the force on a is a vector whose direction is opposite to u
0
andtendstobringaandbcloser(attractiveforce),whileittendstobringthem
apart (repulsive force) when r <r .
0
→− →−
According to the definition of f and f , the sum of the forces applied to a
a b
and b is null (i.e. equilibrium of forces):
→− →−
f +f =0 (10)
a b
Note that no torque is produced as the two forces are colinear.
3 Valence Angles
Valence angles tend to maintain at a fixed value the angle between three atoms
a, b and c such that a is linked to b and b to c, as shown on Fig. 3.1.
Figure 3.1: Valence angle
The forces applied to the three atoms all belong to the plane abc defined by
the points a, b, c.
A harmonic valence potential is a scalar field U which defines the potential
energy of an atom configuration forming a valence angle θ by:
U(θ)=k(θ−θ )2 (11)
0
wherek isthestrengthofthevalenceangleandθ istheequilibriumangle(the
0
oneforwhichenergyisnull). ThepartialderivativeofU accordingtotheangle
θ is thus:
∂U(θ)
=2k(θ−θ ) (12)
∂θ 0
4
The partial derivative of U according to the position r of a is:
a
∂U(θ) ∂U(θ) ∂θ
= . (13)
∂r ∂θ ∂r
a a
that is:
∂U(θ) ∂θ
=2k(θ−θ ). (14)
∂r 0 ∂r
a a
As a describes a circle with radius |ab|, centered on b, we have1:
∂θ 1
= (15)
∂r |ab|
a
→− →−
Let p be the normalized vector in the plane abc, orthogonal to ba :
a
→− →− →− →−
p =norm(ba×(ba×bc)) (16)
a
The force applied on a is then:
→− ∂U(θ) →− →−
f =− .p =−2k(θ−θ )/|ab|.p (17)
a ∂r a 0 a
a
In the same way, the force applied on c is:
→− →−
f =−2k(θ−θ )/|bc|.p (18)
c 0 c
→− →−
where p is the normalized vector in plane abc, orthogonal to cb :
c
→− →− →− →−
p =norm(cb×(ba×bc)) (19)
c
The sum of the forces should be null:
→− →− →−
f +f +f =0 (20)
a b c
Thus, the force applied to b is:
→− →− →−
f =−f −f (21)
b a c
3.1 Torques
→−
Let us now consider torques (moment of forces). The torque exerted by f on
a
→− →− →− →− →− →− →−
b is ba×f and the torque exerted by f on b is bc × f . As ba and f are
a c c a
orthogonal, one has2:
→− →− →− →−
|ba×f |=|ba||f |=|ba||−2k(θ−θ )/|ab||=|2k(θ−θ )| (22)
a a 0 0
1 The length of an arc of circle is equal to the product of the radius by the angle (in
radians)correspondingtothearcofcircle.
2 Ifu⊥v then|u×v|=|u||v|.
5
For the same reasons:
→− →−
|bc×f |=|2k(θ−θ )| (23)
c 0
→− →− →− →−
Thus, the two vectors ba×f and bc ×f have the same length. As they are
a c
by construction in opposite directions, the sum of the two torques on b is null:
→− →− →− →−
ba×f +bc×f =0 (24)
a c
As a consequence, no rotation around b can result from the application of the
→− →−
two forces f and f .
a c
4 Torsion Angles
A torsion angle θ defined by four atoms a,b,c,d is shown on Fig. 4.1. In the
OPLS [3] force-field, as in many other force-fields, potentials of torsion angles
have a “triple-cosine” form. This means that the potential U of a torsion angle
θ is defined by3:
U(θ)=0.5[A (1+cos(θ))+A (1−cos(2θ))+A (1+cos(3θ))+A ] (25)
1 2 3 4
Figure 4.1: Torsion angle θ
Thepartialderivativeofthetorsionanglepotentialaccordingtotheposition
r of a is:
a
∂U(θ) ∂U(θ) ∂θ
= . (26)
∂r ∂θ ∂r
a a
The partial derivative of the potential according to the angle θ is:
∂U(θ)
= 0.5(−A sin(θ)+2A sin(2θ)−3A sin(3θ)) (27)
∂θ 1 2 3
= −0.5(A sin(θ)−2A sin(2θ)+3A sin(3θ)) (28)
1 2 3
3 Inthefollowing,wewillalwaysomitthelastparameterA4.
6
4.1 Forces on a and d
Let us call θ1 the angle a(cid:99)bc. Atom a turns around direction bc, on a circle of
radius |ab|sin(θ ). The partial derivative of θ according to the position of a is:
1
∂θ 1
= (29)
∂r |ab|sin(θ )
a 1
We thus have:
∂U(θ) −0.5
= (A sin(θ)−2A sin(2θ)+3A sin(3θ)) (30)
∂r |ab|sin(θ ) 1 2 3
a 1
Similarly, for atom d, noting θ2 the angle b(cid:99)cd :
∂U(θ) −0.5
= (A sin(θ)−2A sin(2θ)+3A sin(3θ)) (31)
∂r |cd|sin(θ ) 1 2 3
d 2
→− →−
Letp thenormalizedvectororthogonaltotheplaneabc,andp thenormalized
1 2
→− →−
vector orthogonal to the plane bcd (the angle between p and p is θ):
1 2
→− →− →−
p =norm(ba×bc) (32)
1
→− →− →−
p =norm(cd×cb) (33)
2
The force applied on a is:
→− 0.5 →−
f = (A sin(θ)−2A sin(2θ)+3A sin(3θ)).p (34)
a |ab|sin(θ ) 1 2 3 1
1
In the same way, the force applied on d is:
→− 0.5 →−
f = (A sin(θ)−2A sin(2θ)+3A sin(3θ)).p (35)
d |cd|sin(θ ) 1 2 3 2
2
4.2 Forces on b and c
→− →−
We now have to determine the forces f and f to be applied on b and c. The
b c
equilibrium conditions imply two constraints: (A) the sum of the forces has to
be null:
→− →− →− →−
f +f +f +f =0 (36)
a b c d
and (B) the sum of torques also has to be null4. Calling o the center of bond
bc, this means:
−→ →− −→ →− →− →− →− →−
oa×f +od×f +ob×f +oc×f =0 (37)
a d b c
From (37) it results:
→− →− →− →− →− →− →− →− →− →−
(ob+ba)×f +(oc+cd)×f +ob×f +oc×f =0 (38)
a d b c
−→ −→ −→ −→
4 Itisnotpossibletosimplydefine fb =−fa and fc =−fd,asthesumoftorqueswould
benon-null,thusleadingtoanincreaseofpotentialenergy.
7
and:
→− →− →− →− →− →− →− →− →− →−
(−oc+ba)×f +(oc+cd)×f −oc×f +oc×f =0 (39)
a d b c
which implies:
→− →− →− →− →− →− →− →− →−
oc×(−f +f −f +f )+ba×f +cd×f =0 (40)
a d b c a d
From (36) it results:
→− →− →− →− →− →−
−f +f −f +f =2(f +f ) (41)
a d b c d c
Substituting (41) in (40), one gets:
→− →− →− →− →− →− →−
oc×(2(f +f ))+ba×f +cd×f =0 (42)
d c a d
thus:
→− →− →− →− →− →− →− →−
2oc×f +2oc×f +ba×f +cd×f =0 (43)
d c a d
which implies:
→− →− →− →− →− →− →− →−
2oc×f =−2oc×f −cd×f −ba×f (44)
c d d a
→−
and we finally get the condition that the torque from f should verify in order
c
(37) to be true:
→− →− →− →− →− →− →− →−
oc×f =−(oc×f +0.5cd×f +0.5ba×f ) (45)
c d d a
Let us state:
→− →− →− →− →− →− →−
t =−(oc×f +0.5cd×f +0.5ba×f ) (46)
c d d a
→− →− →− →−
Equation oc × x = t has an infinity of solutions in x, all having the same
c
→−
component perpendicular to oc. We thus simply choose as solution the force
→−
perpendicular to oc defined by:
→−f =(1/|oc|2)→−t ×→−oc (47)
c c
Equation (45) is verified because:
→−oc×→−f =(1/|oc|2)→−oc×(→−t ×→−oc) (48)
c c
thus5 :
→−oc×→−f =(1/|oc|2)|oc|2→−t =→−t (49)
c c c
→−
Thevalueof f isfinallydeducedfromequation(36)statingtheequilibrium
b
of forces:
→− →− →− →−
f =−f −f −f (50)
b a c d
→− →− →− →−
Wehavethusdeterminedfourforcesf ,f ,f ,f whosesumisnull(36)and
a b c d
whose sum of torques is also null (37).
5 ifu⊥v,thenu×(v×u)=|u|2v.
8
5 Lennard-Jones Potentials
A Lennard-Jones (LJ) potential U(r) between two atoms placed at distance r
is defined by6:
σ 12 σ 6
U(r)=4(cid:15)[( ) −( ) ] (51)
r r
Inthisdefinition,parameterσ isthedistanceatwhichthepotentialisnull,and
parameter (cid:15) is the minimum of the potential (corresponding to the maximum
of the attractive energy).
Stating A=σ12 and B =σ6, Eq. (51) becomes:
A B
U(r)=4(cid:15)( − ) (52)
r12 r6
The partial derivative of U according to distance is thus:
∂U(r) A B
= 4(cid:15)(−12 +6 ) (53)
∂r r13 r7
A B
= 24(cid:15)(−2 + ) (54)
r13 r7
24(cid:15) A B
= (−2 + ) (55)
r r12 r6
24(cid:15) σ 12 σ 6
= − (2( ) −( ) ) (56)
r r r
Let a and b be two atoms. The force on a is:
→− 24(cid:15) σ 12 σ 6 →−
f = (2( ) −( ) ).u (57)
a r r r
→− →−
where u is the normalization of ba. From the action/reaction principle, one
deduces that the force on b should be the opposite of the force on a:
→− →−
f =−f (58)
b a
→− →−
According to the definition of f and f , the sum of the forces applied to a
a b
and b is null:
→− →−
f +f =0 (59)
a b
Note that, as for bonds, no torque is produced because the two forces are col-
inear.
6 Resume
Theforcesdefinedintheprevioussectionsaresummedupinthefollowingtable:
6 thisisthe6-12form;otherformsofLJpotentialsexist.
9
→− →−
Bond ab 8 f =−2k(r−r ).u
a 0
→− →−
9 f =−f
b a
→− →−
Valence abc 17 f =−2k(θ−θ )/|ab|.p
a 0 a
→− →− →−
21 f =−f −f
b a c
→− →−
18 f =−2k(θ−θ )/|bc|.p
c 0 c
Torsion abcd 34 →−f = 0.5 (A sin(θ)−2A sin(2θ)+3A sin(3θ)).→−p
→−a |a→−b|sin(θ→−1) 1→− 2 3 1
50 f =−f −f −f
b a c d
→− →− →−
47 f =(1/|oc|2).t ×oc
c c
35 →−f = 0.5 (A sin(θ)−2A sin(2θ)+3A sin(3θ)).→−p
LJ ab 57 →−fd = |2c4d(cid:15)|s(i2n((σθ2))12−1(σ)6).→−u 2 3 2
→−a r→− r r
58 f =−f
b a
Bond InEq. 8,k isthebondstrengthconstant,r isthedistancebetweenatoms
aandb,andr istheequilibriumdistance,forwhichenergyisnull. Vector
→− 0 →− →−
u is defined by u =norm(ba).
Valence In17and18, k istheanglestrengthconstant,θ istheanglea(cid:99)bc,andθ0 is
→− →−
theequilibriumangle,forwhichenergyisnull. In17,p isdefinedbyp =
→− →− →− →− →− a →− →− →−a
norm(ba×(ba×bc)). In18, p isdefinedby p =norm(cb×(ba×bc)).
c c
Torsion In 34 and 35, θ is the torsion angle, θ1 is the angle a(cid:99)bc, θ2 is the angle
b(cid:99)cdandA1, A2 andA3 arethepara→−meterswhichdefi→−nethe“three→−-cosin→−e”
form of the torsion angle. Vector p is defined by p = norm(ba×bc)
1 1
→− →− →− →−
and p =norm(cd×cb). In 47, o is the middle of bc and t is defined by
2 c
→− →− →− →− →− →− →−
t =−(oc×f +0.5cd×f +0.5ba×f ).
c d d a
LJ In 57, σ is the distance at which the potential is null and (cid:15) is the depth of
→− →−
thepotential(minimumofenergy). Asforbonds,onehas u =norm(ba).
In each case (bond, valence, torsion, LJ interaction), the sum of the forces
that are applied to atoms is always null (Eq. (10), (20), (36), (59)). Moreover,
no torque is induced by application of these forces: no torque is produced by
bonds and LJ interactions, as the produced forces are colinear; we have verified
inSec. 3.1thatnotorqueisproducedbyvalenceangles(24); fortorsionangles,
we have chosen the forces in such a way that the sum of the forces and the
global sum of torques are always null (37). This means that no energy is ever
added by the application of the forces during the simulation process.
7 Conclusion
We have precisely defined the forces that apply on atoms in MD simulations.
The definitions are given in a purely vectorial formalism (no use of a specific
10