Transcript Ionic and Molecular Liquids - Jawaharlal Nehru Centre for

Ionic and Molecular Liquids:
Long-Range Forces and
Rigid-Body Constraints
Charusita Chakravarty
Indian Institute of Technology Delhi
1
Electrostatic Interactions
• Ionic liquids: NaCl, SiO2, NH4Cl
• mobile charge carriers which are
atomic or molecular entities
• Simple ionic melts: model ions as point
charges with Coulombic interactions.
Short-range repulsions control ionic
radii.
water
• Molecular Liquids
• Electronic charge distributions show
significant deviations from spherical
symmetry
• Can be modeled by: (a) multipole
moment expansions or (b) arrays of
partial charges
benzene
2
Multipole Expansion for Electrostatic
Potential
Represent the electronic charge
distribution of a molecule by a set
of multipoles:
Q

4r

Coulombic
 cos
2
dipolar
r
 3 cos2   1

  quadrupolar
3
2
r


3
Range of Electrostatic Interactions
Type
Range
Energy (kJ/mol)
Comment
Ion-Ion
1/r
250
Ion-dipole
1/r2
15
Dipole-Dipole
1/r3
Dipole-quadrupole
1/r4
Fixed Orientation / Linear
Quadrupole-quadrupole 1/r5
Fixed Orientation / Linear
2 Static dipoles (solid phase)
Long-range interactions: Tail correction will diverge for 1/rn interactions with
n greater than or equal to 3; therefore minimum image convention cannot
be applied
4
Partial Charge Distributions of Some Typical
Molecules
Dipole + Quadrupole
2d
(-0.8472e)
O
109.47
H
d
(+0.52e)
2d
(1.04e)
TIP4P Water
Hd

(+0.52e)
H
OO
H
109.47
d
d
(+0.4238e)
(+0.4238e)
SPC/E Water
Quadrupole
Molecular
Nitrogen
5
Array of Point Charges
Coulombic contribution to the potential energy for an array
of N charges that form a charge-neutral system:
Electrostatic potential
• Particle i interacts with
all other charges and their mirror images but not with itself
• Gaussian units
• Cannot apply minimum image convention because
sum converges very slowly
6
Poisson Equation


U (r ) and U(r )

 (r )

 (r )
Charge
Distribution
Potential
Electrostatic
Potential
Energy and
Forces


  (r )  4(r )
2
Linear differential equation:

Charge distribution :  i (r )
i

Electrostatic potential:  i (r )
i
7
Ewald Summation for Point Charges

Co P (r )
Screened charge
distribution:
Converges fast in real
space


 P (r )  G (r )

- G (r )
Point Charge
Distribution:
Converges slowly
Gaussian compensating
charge distribution:
can be analytically
evaluated in real space
8
Ewald Summation
•
Screening a point charge to convert the long-range
Coulombic interaction into a short-range interaction
•
Evaluating the real-space contribution due to the
screened charges
•
The Poisson equation in reciprocal space for the
compensating screened charge distribution
•
Evaluating the reciprocal space contribution
•
Self-interaction correction
9
Screening a point charge


  (r )  4(r )
2
Single pointchargecentredat origin :



 SP (r )  qid (r ) and  SP (r )  qi / r
Single Gaussian chargecentredat origin :

 SG (r )  qi ( /  )3 / 2 exp(r 2 )
Electrostatic pot entialdue to Gaussian charge
dist ribution :

1   2  SG (r ) 
 
  4SG (r )
2
r
r

10
Electrostatic potential due to a Gaussian charge
distribution
11
Electrostatic potentialof pointcharge screened
by a compensating Gaussian chargedensity:
 short -range  (qi / r )  (qi / r )erf (  r )  (qi / r )erfc(  r )
12
Real Space Contribution


Elect rost at ic int eract ion at rj due t o charge qi at rj :



q
SR (rj )  erfc  rij
rij
Int eract io
n energy of an arrayof chargesq j locat ed
at posit ionsrj wit h t heelect rost at ic pot ent ialat rj :
U real
U real

1
  q jSR (rj )
2 j

qi q j
1
 
erfc  rij
2 i j i rij

The value of  must be chosen so that the range of interaction of the
screened charges is small enough that a real space cutoff of rc < L/2
can be used and the minimum image convention can be applied
13
Poisson Equation in Reciprocal Space
Fourier series representation of a function in a cubic box with edge length L
and volume V under periodic boundary conditions:

 


~ 
f (r )  (1 / V ) f (k ) exp(ik  r ) where k  2 (l x , l y , l z )
l 
Fourier coefficients

 
 
f ( k )   dr f ( r ) exp( ik  r )
V


Poisson equation in real space    r   4(r )
2~
k  (k )  4~(k )
Poisson equation in reciprocal space
2
Reminder: Fourier transforms of derivatives
~
df
FT of
 kf ( k )
dX

dn f
n ~


FT of


k
f (k )
n
dX
14



4
~
~
Unit Point Charge at origin:
 S (r )  d (r ) g ( k )  2
k

 
~
 P (r )   qid (r  ri )
Array of point charges

 
i
 P (k )   qi exp(k  ri )
i
 



4
~
P (k )  g (k )  P (k )  2  qi exp(k  ri )
k
i
Array of Gaussian charges
Convolution of Point Charge distribution and smearing function

   dr  ( r )  ( r  r )
G (r )   qi (r  ri )
P


~  ~ i ~ 
 (k )  g (k )r (k )P (k )
Green’s
Function
FT of Point charge Array
Fourier transform
of smearing Function
~
~ ~
Reminder: f1 ( x)  f 2 ( x)  f3 ( x) then f1 ( x)  f 2 ( x) f3 ( x)
15
Fourier Part of Ewald Sum
Array of N Gaussian point charges with periodic images:
N

   2
32
G (r )   q j    exp   r  rj  nL
j 1




n
Corresponding Electrostatic potential in reciprocal space
 4
 (k )  2
k

 
 
2
exp


k

r
exp

k
4

i
N
i 1

Electrostatic potential in real space can be obtained using:

 
1
 (r )   (k ) exp(k  r )
V k 0
16
Fourier Part of Ewald Sum (contd.)
Electrostatic potential in real space

 
1
 (r )    (k ) exp(k  r )
V k 0
N 4q
  
j
  2 exp[ik  (r  rj )]exp(k 2 / 4 )
k  0 j 1 k
Reciprocal Space contribution to potential energy

U rec  (1 / 2) qi  (r )
i
N
 (1 / 2) 
4q j qi
  
exp[ik  (ri  rj )]exp(k 2 / 4 )
Vk 2
 2
1
4
2


(
k
)
exp

k
4

2
2V k  0 k
k  0 i , j 1


17
System is embedded in a medium with an inifinite dielectric constant
Correction for Self-Interaction
Must remove potential energy contribution due to a
continuous Gaussian cloudof charge qi and a point charge qi
located at the centre of the cloud.
Electrostatic potential due to Gaussian centred at origin


qi
erf  r
r
2q 
Gauss (r  0)  i
r 
Gauss (r ) 
U Self
1 N
  qiself (ri )
2 i 1
   
12
N
q
i 1
2
i
18
Coulombic Interaction expressed as an
Ewald Sum
U Coul
1

2V

 2
4
2

k
exp

k
4


2

k 0 k
   
1 N
2

SelfInteraction
q
i 1
2
i

1 N qi q j erfc  rij
 
2 i j
rij

Reciprocal
space

Real Space
Important: For molecules, the self-interaction correction must be modified
because partial charges on the same molecule will not interact
19
Accuracy of Ewald Summation
• Convergence parameters:
– Width of Gaussian in real space, 
– Real space cut-off, rc
– Cutoff in Fourier or reciprocal space, kc=2/Lnc
rc  s / 
nc  sL / 
20
Calculating Ewald Sums for NaCl
Na+
Cl-
Liquid
Crystal
21
Hands-on Exercise:
Calculating the Madelung Constant for NaCl
The electrostatic energy of a structure of 2N ions of charge +/- q is
q2
 q2
U   N
 N 
rnn
rij
where  is the Madelung constant and rnn is the distance between the
nearest neighbours.
Structure

Sodium Chloride (NaCl)
1.747565
Cesium Chloride (CsCl)
1.762675
Zinc blende (cubic ZnS)
1.6381
22
Rotational and Vibrational Modes of
Water
Rotational
Asymmetric Strech
3756cm-1
Intermolecular
vibrations
n (cm-1)
Librations
800
OO stretch
200
OOO bend
60
http://www1.lsbu.ac.uk/water/vibrat.htm
Bend
1595cm-1
http://chsfpc5.chem.ncsu.edu/~franzen/CH795N/lecture/XIV/XIV.html
Symmetric Stretch
3657cm-1
Constants
n (cm-1)
A
40.1
B
20.9
C
13.4
23
Molecules: Multiple Time-scales
• Bonded interactions are much stronger than non-bonded
interactions
• Intramolecular vibrations have frequencies that are
typically an order magnitude greater than those of
intermolecular vibrations
• MD/MC: time step will be dictated by fastest vibrational
mode
• Fast, intramolecular vibrational modes do not explore
much of configuration space- rapid, essentially harmonic,
small amplitude motion about equilibrium geometry
• Require efficient sampling of orientational and
intermolecular vibrational motions
24
Simulation Methods for Molecules
• Freeze out all or some intramolecular modes:
– Serve to define vibrationally averaged molecular
structure and are completely decoupled from
intermolecular vibrations, librations or rotations
– Rigid-Body Rotations:
• Characterize the mass distribution of the molecule by its
moment of inertia tensor
• MC: Use orientational moves
• MD: Propagate rigid-body equations of motion
• Will not work if there are low-frequency vibrational modes
– Apply Geometrical Constraints
• MD: SHAKE
• MD: RATTLE
• Multiple Time-Scale Algorithms
25
Rotations of Rigid-Bodies
Space-fixed (SF) and Body-fixed (BF)
axes (Goldstein, Classical Mechanics)
N2
H2
CH4
Moments of Inertia of Molecules:
Ia < Ib < Ic
Linear:
Spherical polar angles: ,
Non-linear:
Euler angles: ,,y
(Atkins, Physical Chemistry)
NH3
SF6
lCl5
H2O
Ix  Iy26 Iz
Monte Carlo Orientational Moves for
Linear Molecules
Orientation of a linear molecule is specified by a unit vector
u . To change it by a small amount:
1.
Generate a unit vector v with a random orientation.
See algorithm to generate random vector on unit
sphere
2.
Multiply vector v with a scale factor g, which
determines acceptance probability of trial orientational
move
3.
Create a sum vector: t = u + gv
4.
Normalise t to obtain trial orientation
27
Euler Angles (,,y
Euler’s rotation theorem: Any rotation of a rigid-body may be
described by a set of three angles
•Rotation, A: Initial orientation of body-fixed axes (X,Y,Z) to final orientation (X’’,Y’’’, Z’)
http://stackoverflow.com/questions
28
http://mathworld.wolfram.com/EulerAngles.html
Euler angles (,,y(contd.)
•
Rotate the X-, Y-, and Z-axes about the Z-axis by   <  <  ,
resulting in the X'-, Y'-, and Z-axes.
•
Rotate the X'-, Y'-, and Z-axes about the X'-axis by  0<<,
resulting in the X'-, Y''-, and Z'-axes.
•
Rotate the X'-, Y''-, and Z'-axes about the Z'-axis by y   <y< ,
resulting in the X''-, Y'''-, and Z'-axes.
Rotation A = BCD, therefore new coordinates are:
29
Monte Carlo Orientational Moves for
Non-linear Molecules
• Specify the orientation of a rigid body by a quaternion of
unit norm Q that may be thought of as a unit vector in
four-dimensional space.
30
Applying Geometrical Constraints
Lagrangian Equations of Motion
  L   L 

  

t  q   q 
Cartesian coordinates
1 2

L( x, x)  mx  U ( x)
2
  L   L 



t  x   x 
 U 
mx  

 x 
Kinetic Energy
Potential Energy
L  K (q )  U (q)
Geometrical Constraints
  (q)  0
 
d   
dqi
i qi
 
 0  q     0
t
Define constraint equations and require that system moves tangential to
the constraint plane.
31
Introduce a new Lagrangian that contains all the constraints:
L'  L     q 

The equations of motion of the constrained system are:
  L'   L' 

  

t  q   q 
 
U
i  
mi q
  
qi
qi

 Fi   Gi

Constraint force acting
along coordinate qi
32
Bond Length Constraint for Diatomic
Molecule
m1
m2
 2
Bond constraint   r1  r2  d 2  0
Constraint forces:
•
•
•
•
lie along bond direction
are equal in magnitude
opposite in direction
do no work
Verlet algorithm:
 
 

  2(r1  r2 ) or G1   (r1  r2 )
r1
 
 

  2(r1  r2 ) or G2   (r1  r2 )
r2



t 2
t 2
ri (t  t )  2ri (t )  ri (t  t ) 
Fi (t ) 
Gi (t )
mi
mi
New position
with constraint

t 2
 ri ' (t  t ) 
Gi (t )
mi
Unconstrained
position
Constraint forces
on atom i
33
Diatomic molecule (contd)
Three unknowns




2
r1 (t  t )  r1 ' (t  t )   (t mi )(r1 (t )  r2 (t ))




2
r2 (t  t )  r2 ' (t  t )   (t mi )(r1 (t )  r2 (t ))


2
r1 (t  t )  r2 (t  t )  d 2
In the case of a diatomic molecule, can obtain quadratic
equation in . However, cannot be conveniently generalised
for larger molecules with more constraints.
34
Multiple Constraints
Verlet algorithm for N atoms in the presence of l constraints:
Taylor expansion of the constraints with
respect to unconstrained positions.
Impose the requirement that after one time
step, the constraints must be satisfied.
i k (t  t ) using
unconstrained positions
Constraint force acting
on i-th atom due to k-th
constraint
Substitute from Taylor
35
expansion above
Multiple Constraints (contd.)
Need to find the unknown Lagrange multipliers:
Solution by Matrix Inversion
Since the Taylor expansion was truncated at first-order
iterative scheme will be required to obtain convergence.
36
SHAKE Algorithm
• Iterative scheme is applied to each constraint in turn:
RATTLE: Similar approach within a velocity Verlet scheme
37
References
• D. Frenkel and B. Smit, Understanding
Molecular Simulations: From Algorithms to
Applications
• A. R. Leach, Molecular Modelling: Principles and
Applications
• D. C. Rapaport, The Art of Molecular Dynamics
Simulation (Details of how to implement
algorithms for molecular systems)
• M. P. Allen and D. J. Tildesley, Computer
Simulation of Liquids (SHAKE, RATTLE, Ewald
subroutines)
38