Transcript geneste

Path Integral Molecular Dynamics
Grégory GENESTE
CEA, DAM, DIF, F-91297 Arpajon, France
14/04/2011
Workshop Abinit 2011
1
Coworkers:
Marc Torrent (parallelization over images)
François Bottin
Jessica Hermet (proton transport in oxides)
Hichem Dammak, labo SPMS, Ecole Centrale Paris
Marc Hayoun, LSI, Ecole Polytechnique
14/04/2011
Workshop Abinit 2011
2
Outline
Introduction: why studying quantum effects associated
to atomic motions ?
I - Path Integral (PI) methods
II - Technical points related to Path Integral Molecular
Dynamics (PIMD)
III – Implementation of PIMD in ABINIT, keywords, etc
IV – Langevin PIMD, choice of random number generator,
friction coefficient
V – Conclusion: in progress and future work
14/04/2011
Workshop Abinit 2011
3
Why studying quantum effects of atomic motions ?
Quantum effects of atomic motions are important
Below the Debye temperature
- CRUCIAL in the case of LIGHT elements, up to
very high temperatures:
H2, Debye T ~ 6000 K (vibration)
He: liquid down to zero K (under 1 atm) because of
quantum effects
- IMPORTANT for « ordinary solids » with
Debye temperature ~ RT or above
14/04/2011
Workshop Abinit 2011
4
Path Integral Methods
One way to account for quantum effects =
Path Integral (PI) methods
Based on the discretized version of Feynman's path
Integral applied to the density operator.
Path integral for time evolution operator:
t
x'
∣ x ' 〉= ∫ Dx s e
− i H t /ℏ
〈 x∣e
i ds L x s , dx s
ℏ ∫0
ds
x'
= ∫ Dx s e
x
x
… and for the canonical density operator:
ℏ
− 1 ∫ ds[ m dx
ℏ0
2 ds
∣ x ' 〉= ∫ Dx s e
− H
〈 x∣e
x'
i
A x s
ℏ
t= − i ℏ
2
V x s ]
x
14/04/2011
Workshop Abinit 2011
5
Path Integral Methods
Canonical partition function:
ℏ
− 1 ∫ ds [ m dx
ℏ0
2 ds
x
Z = ∫ dx ∫ Dx s e
x
x 0=x
ℏ
2
V x s ]
− 1 ∫ ds[ m dx
ℏ0
2 ds
= ∮ Dx s e
ℏ=x
2
V x s ]
ℏ
x 0=x
Discretized version:
Z = lim P
∞[
2
m P k B T 3NP/ 2
]
∫
2
h
R0
...∫
RP −1
e−
V eff R 0 ... RP− 1
d R0 ... d R P− 1
with
R s= r 1 s ... r N s
3N− vector
V eff R 0 ... R P− 1 = V eff r 1 0 ... r N 0 ... r 1 P− 1 ... r N P− 1
P− 1
N
V eff R 0 ... R P− 1 = ∑ [ ∑
s= 0
14/04/2011
i= 1
1
k P,
2
ri s − ri s
Workshop Abinit 2011
1 2
1
V r 1 s ...r N s ]
P
6
Path Integral Methods
Spring constant
2 2
m P mPkBT
= 2 2=
ℏ
ℏ2
k P,
P = « Trotter » number
ri P = ri 0
with cyclic boundary condition:
Z = lim P
2 m P k B T 3NP/ 2
]
∫
∞[
2
h
R0
...∫
RP −1
e−
V eff R 0 ... RP− 1
d R0 ... d R P− 1
The right member is the CLASSICAL partition function (for fixed Trotter number P)
of the following CLASSICAL system (of NxP particles):
(up to a multiplicative constant)
14/04/2011
Workshop Abinit 2011
7
Path Integral Methods
Example (P=4, N=5)
s=0
s=1
s=3
i=5
i=4
s=2
i=3
i=2
i=1
Harmonic interaction, with
k P,
2 2
m P mPkBT
= 2 2=
ℏ
ℏ2
« Real potential » divided by P
s = « imaginary time »
Ensemble of « quasi-particles » for fixed s = « imaginary time slice »
14/04/2011
Workshop Abinit 2011
8
Path Integral Methods
Each quantum system has, for fixed P, a classical equivalent
= « CLASSICAL ISOMORPHISM »
In the limit of infinite Trotter numbers, the physical EQUILIBRIUM properties
of the quantum system are the same as those of the classical equivalent.
One gets directly the observables in the CANONICAL ensemble.
Can be extended easily to the isothermal-isobaric (NPT) ensemble.
Ex: Barrat, Loubeyre, Klein, J. Chem. Phys. 90, 5644 (1989)
Thus, the EQUILIBRIUM properties of the classical equivalent system can be
reached by any CLASSICAL simulation technique
Monte Carlo Method
= PIMC
14/04/2011
Molecular Dynamics Method
=PIMD
Workshop Abinit 2011
9
Path Integral Methods
Path Integral is an elegant way to do Quantum Mechanics …
… without wave function !!!
NB: at low temperature, the Trotter number should be chosen
very high
=> one recovers the complexity of QM !
14/04/2011
Workshop Abinit 2011
10
Path Integral Methods
Limitations
2 m P k B T 3NP/ 2
Z = lim P ∞ [
]
∫
h2
Is obtained from
Z=∫ R R, R; d R
where
R,R' ;
R0
...∫
RP −1
e−
V eff R 0 ... RP− 1
d R0 ... d R P− 1
is the density matrix in the position representation:
R,R' ;
= 〈 R∣e− H∣R ' 〉
Such a formulation implicitely assumes DISCERNABLE particles.
In the case of INDISTINGUISHABLE particles, Z writes:
1) for BOSONS:
1
B
Z =
∑n
N ! p∈S
∫R
R , p R;
Permutation signature
2) for FERMIONS:
F
Z =
14/04/2011
dR
1
N ! p∈∑
Sn
p
∫R
R , p R;
Workshop Abinit 2011
dR
11
Path Integral Methods
IN PRACTISE: Primitive Approximation (PA)
For fixed P:
Z≈ [
2
m P k B T 3NP/ 2
]
∫
h2
Z = lim P
∞
R0
...∫
RP− 1
e−
V eff R 0 ... R P− 1
d R0 ... d R P− 1
ZP
2 m P k B T 3NP/ 2
Z P= [
]
∫
h2
R0
... ∫
R P− 1
e−
V eff R0 ... R P− 1
d R 0 ... d R P− 1
Calculations are performed for fixed P
Observables are computed using ZP, i.e. in the PA.
Quantities must be CONVERGED with P
<A> computed from ZP = Primitive Estimator of A
14/04/2011
Workshop Abinit 2011
12
Path Integral Methods
Example: Energy Primitive Estimator:
Using
2 m P k B T 3NP/ 2
Z P= [
]
∫
2
h
R0
... ∫
R P− 1
e−
∂ ln Z P
P= −
∂
E
V eff R0 ... R P− 1
d R 0 ... d R P− 1
one gets:
E
3
=
NP k B T −
P
2
P− 1 N
∑∑
s= 0 i= 1
1
k P,
2
s
r i − ri
s 1 2
Kinetic Energy Estimator
14/04/2011
Workshop Abinit 2011
1
P
P− 1
∑V
r 1s ... r Ns
s= 0
Potential Energy Estimator
13
Technical points related to PIMD
Z = lim P
2 m P k B T 3NP/ 2
]
∫
∞[
2
h
P− 1
N
V eff R 0 ... R P− 1 = ∑ [ ∑
s= 0
k P,
i= 1
R0
...∫
1
k P,
2
RP −1
e−
ri s − ri s
V eff R 0 ... RP− 1
1 2
d R0 ... d R P− 1
1
V r 1 s ...r N s ]
P
2 2
m P mPkBT
= 2 2=
ℏ
ℏ2
This term grows to infinity
as P tends to infinity
tends to ZERO as
P tends to infinity!
In the limit of large P, the classical system is mainly
HARMONIC
Standard MD algorithms such as Nose-Hoover thermostat
fail to correctly sample the phase space => non ergodic trajectory
14/04/2011
Workshop Abinit 2011
14
Technical points related to PIMD
How to recover ergodicity in PIMD ?
Marx, Parrinello, J. Chem. Phys. 104, 4080 (1996)
Tuckerman, Marx, Klein, Parrinello, J. Chem. Phys. 104, 5579 (1996)
1) Thermostat chains (Martyna, Tobias, Klein, J. Chem. Phys. 101, 4177 (1994))
+ coordinate transformations
(staging or normal mode)
=> the corresponding scheme has been extended to the NPT ensemble
Equations of PIMD in the NPT ensemble:
Martyna, Hugghes, Tuckerman, J. Chem. Phys. 110, 3275 (1999)
2) Alternative way: Langevin dynamics
Uses random numbers
14/04/2011
Workshop Abinit 2011
15
Technical points related to PIMD
LANGEVIN DYNAMICS
NVT: « Langevin thermostat »
with
s
f i =−
d2 r is
d ri s
s
m
=fi − m
dt
dt 2
1
s
s
∇ V r r 1 ... r N − k P ,
P
i
s
2 ri − ri
Ri s
s 1
− ri
s− 1
s
d ri
− m = friction force
dt
Ris = random « Langevin » force (white noise)
Langevin force: Quigley, Probert, Comput. Phys. Comm. 169, 322 (2005)
Gaussian with variance
2 m kB T
t
Time step
14/04/2011
Workshop Abinit 2011
16
Technical points related to PIMD
Example: (3D) harmonic oscillator
3
ℏ = 51.44 meV
2
PxT = 1600
Energy (meV)
P=1
(3 because 3D
oscillator)
P=320
(PIMD)
(exact)
Temperature (K)
14/04/2011
Workshop Abinit 2011
17
Technical points related to PIMD
How to implement the Langevin method in the NPT ensemble ?
This has been formulated by Quigley and Probert
J. Chem. Phys. 120, 11432 (2004); Comput. Phys. Comm. 169, 322 (2005)
d p is
= f is −
dt
d pi s
dt
d r i s pi s
=
dt
m
pG s
r
Wg i
dh p G h
=
dt W g
dpG
= V X− Pexp Id
dt
pG s
1 Tr [ pG ] s
Ri −
p −
pi
Wg i
Nf
Wg
s
1
Nf
∑
i,s
pi s 2
Id−
m
G
pG R p
in which X is an « internal pressure » defined by
V X = ∑ pi s,
i,s
'
14/04/2011
=
∂
∂h
pi s,
m
r is, f i s, −
' hT
Combines Langevin dynamics and
Martyna et al barostat
Workshop Abinit 2011
18
Implementation in ABINIT
- ABINIT already contains a structure adapted to PIMD
- this structure is devoted to algorithms using different
IMAGES of the system
- these images are propagated at the same time, according
to the chosen algorithm
- well adapted to NEB, string method, etc
- the two subroutines contained in ABINIT that control these
algorithms:
gstateimg.F90
predictimg.F90
First, let us have a look at these two existing routines:
14/04/2011
Workshop Abinit 2011
19
Implementation in ABINIT
Subroutine gstateimg.F90, general structure :
Loop on itimimage (propagation of images)
do itimimage=1,ntimimage
1) loop on the images (that evolve) to compute energy,
forces...
do idynimage=1,ndynimage
Call gstate : computed total energy, forces, stress
for each image
end do
2) Evolution of images
Call predictimg
end do
14/04/2011
Workshop Abinit 2011
20
Implementation in ABINIT
Subroutine predictimg
=> chooses the algorithm to propagate the images
(according to the value of imgmov)
Other
Algorithms
(string, etc)
select case(imgmov)
case(0)
call predict_copy(acell_timimage,itimimage,list_dynimage,natom,ndynimage,nimage,ntimimage,&
& rprim_timimage,vel_timimage,xred_timimage)
case(1)
&
Here comes
PIMD !
call predict_steepest(acell_timimage,fxcartfactor,itimimage,list_dynimage,natom,ndynimage,nimage,ntimimage,&
results_gs_timimage,rprim_timimage,vel_timimage,xred_timimage)
case(2)
call predict_string(acell_timimage,fxcartfactor,iatfix,itimimage,list_dynimage,natom,ndynimage,nimage,ntimimage,&
& results_gs_timimage,rprim_timimage,vel_timimage,xred_timimage)
case(9, 13)
! Path Integral Molecular Dynamics
call predict_pimd(acell_timimage,fxcartfactor,itimimage,list_dynimage,natom,ndynimage,nimage,ntimimage,&
& results_gs_timimage,rprim_timimage,vel_timimage,xred_timimage,mdparam)
case default
end select
14/04/2011
Workshop Abinit 2011
21
Implementation in ABINIT
PIMD implemented in: SUBROUTINE predict_pimd.F90
Two algorithms:
- Langevin dynamics
- Nose Hoover chains + coordinate transformation
In two ensembles:
- NVT
- NPT
Rq: Langevin NPT uses Quigley-Probert algorithm
14/04/2011
Workshop Abinit 2011
22
Implementation in ABINIT
Associated keywords:
Choice of algorithm:
imgmov=9 (Langevin) or 13 (Nose-Hoover chains)
Choice of statistical ensemble:
optcell = 0 if NVT, 2 if NPT
nimage = TROTTER number
ntimimage = number of time steps
mditemp = initial temperature
mdftemp = thermostat temperature
dtion = time step
Langevin thermostat:
vis= friction coeff in the case of Langevin dynamics
14/04/2011
Workshop Abinit 2011
23
Langevin PIMD
Langevin dynamics is a powerfool way to
recover ergodicity, …
… but special care should be paid to:
1) Choice of random number generator
2) Choice of friction coefficient
14/04/2011
Workshop Abinit 2011
24
Langevin PIMD
1) Random number generator:
Langevin dynamics requires random numbers
The quality of the random number generator is crucial !
Possible choices:
(1) in ABINIT: subroutine uniformrandom.f90
(2) intrinsic function of fortran = random_number
(3) random number generator: ZBQ
(R.Chandler, generator from G. Marsaglia and A. Zaman,
Annals of Appl. Probability 1 (1991), 462-480)
Test of algorithm on a long trajectory => (2) and (3)
give better result than (1)
15/04/2011
14/04/2011
Workshop Abinit 2011
28
25
Langevin PIMD
Best way to test MD algorithm is to perform
A LONG trajectory on a BIG system, but
… difficult in ab initio !
=> remove (temporary) the ab initio part of ABINIT !!!!!
and replace it by a phenomenological potential
=> 100 000 time steps on 1728 particles, T=300 K
Temperature
(K)
(1) Use of subroutine uniformrandom:
14/04/2011
Workshop Abinit 2011
26
Langevin PIMD
(2) Use of intrinsic fortran routine random_number
Temperature
(K)
(3) Use of random number generator ZBQ
Temperature
(K)
14/04/2011
Workshop Abinit 2011
27
Langevin PIMD
totoQ
Changing the
Friction coeff:
random_number
ZBQ
uniformrandom
14/04/2011
Workshop Abinit 2011
28
Langevin PIMD
2) Friction coefficient
- Langevin thermostat requires special attention to the choice
of the friction coefficient
s
d2 r is
d
r
s
i
m
=
f
−
m
i
dt
dt 2
Ri s
Random (Langevin) force randomly
chosen at each step according
to a normal law with variance
2 m kB T
t
- 1/γ roughly corresponds to the time needed to equilibrate the system
The higher γ, the faster your system is equilibrated...
but
- γ (homogeneous to a frequency) should be chosen smaller than
the characteristic phonon frequencies of the system
=> a compromise must be found !
14/04/2011
Workshop Abinit 2011
29
Langevin PIMD
Example: PIMD with ABINIT, NVT Langevin
Box of 20 H2 molecules, P=20, 100 processors (// band-FFT)
Several tests with Vis = 2.10-5, 1.10-4, 2.10-4, 5.10-4
(given in inverse atomic time units)
mditemp= 2000 K; mdftemp = 1000 K, dtion=20
1) vis = 2.10-5 = 0.9 THz
14/04/2011
Workshop Abinit 2011
30
Langevin PIMD
2) vis = 1.10-4 = 4.6 THz
3) vis = 2.10-4 = 9.3 THz
14/04/2011
Workshop Abinit 2011
31
Langevin PIMD
4) vis = 5.10-4 = 23.3 THz
Rq: stretching vibration of H2 molecule = 4160 cm-1 = 125 THz
14/04/2011
Workshop Abinit 2011
32
Langevin PIMD
Be careful to the TIME STEP !
- usually, systems interesting for PIMD contain LIGHT atoms
(H, He, etc) => high phonon frequencies
- The DEFAULT value of DTION is TOO high in such cases
Default value: dtion=100 (~ 2.10-15 s)
14/04/2011
Workshop Abinit 2011
33
CONCLUSION
- PIMD works in NVT and NPT in the case of Langevin
Dynamics, in // (k-point, band, fft)
- Nose-Hoover chains implementation in progress
(NVT and NPT)
- with P=1, one recovers classical dynamics => PIMD
also involves classical MD (=> new algorighms in ABINIT)
- systems of interest:
H2, LiH
proton transport in oxides
14/04/2011
Workshop Abinit 2011
34
CONCLUSION
PARALLELIZATION:
- Loop over images : highly parallelizable !
- currently IMPLEMENTED by M. TORRENT
(for all image algorithms)
=> will provide a 4th level of // in ABINIT
14/04/2011
Workshop Abinit 2011
35