Ab initio molecular dynamics via the Car-Parrinello method: Basic

Document Sample

```					Ab initio molecular dynamics via the Car-Parrinello method: Basic ideas, theory
and algorithms

Mark E. Tuckerman
Dept. of Chemistry
and Courant Institute of Mathematical Sciences
New York University, 100 Washington Sq. East
New York, NY 10003
1808:
“We are perhaps not far removed from the time
when we shall be able to submit the bulk of
chemical phenomena to calculation.”

Joseph Louis Gay-Lussac (1778-1850)
“The underlying physical laws necessary for the mathematical
theory of a large part of physics and the whole of chemistry are
thus completely known, and the difficulty is only that the exact
solution of these laws leads to equations much to complicated to be
soluble.”

Paul Dirac on Quantum Mechanics (1929).
“Every attempt to refer chemical questions to mathematical
doctrines must be considered, now and always, profoundly
irrational, as being contrary to the nature of the phenomena.”

August Comte, 1830
Motivation
• Car-Parrinello is a method for performing molecular dynamics with
forces obtained from electronic structure calculations performed “on
the fly” as the simulation proceeds. This is known as ab initio
molecular dynamics (AIMD).

• As a result, AIMD calculations are considerably more expensive
than force-field calculations, which only involve evaluation of simple
functions of position.

• Force fields, although useful, are, with notable exceptions, unable to
treat chemical bond breaking and forming events.

• Force fields often lack transferability to thermodynamic situations in
which they are not designed to work.

• Polarization and manybody interactions included implicitly.
From ISI Citation Report
R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985)

Total Cites = 4,812
The “Universal” Hamiltonian
M Electrons

N Nuclei               ˆ ˆ ˆ ˆ              ˆ     ˆ
H  Te  Tn  Vee  Ven  Vnn

Operator Definitions:
Electronic:                         Nuclear:
1 M 2                                1 N 1 2
Te    i
ˆ                                   Tn   
ˆ               I
2 i 1                               2 I 1 M I
M                                   N
1                                ZI ZJ
Vee  
ˆ                                Vnn  
ˆ
i j   ri  r j
ˆ ˆ                        I J   ˆ
R R  ˆ
I       J

Coupling:
M    N
ZI
Ven  
ˆ
i 1 I 1   ˆ ˆ
ri  R I
Molecular energy levels
Electron coordinates                Nuclear coordinates

Notation:        r  r1 ,..., rM                      R  R1 ,..., R N
x  r1 , sz ,1 ,..., rM , sz , M

1            0
sz ,i ,  ,   or ,  ,  
0           1

Complete energy level spectrum:

ˆ
H (x, R)  E(x, R)

ˆ ˆ                     ˆ ˆ          ˆ
Te  Tn  Vee (r)  Ven (r, R)  Vnn (R)  (x, R)  E (x, R)
ˆ
                                         
Born-Oppenheimer Approximation
à la W. H. Flygare, Molecular Structure and Dynamics

Mass disparity: M H  2000me
Quasi adiabatic separability ansatz for wave function:

(x, R)   (x, R)  (R)


Schrödinger equation separates if

 I  (R )          I  (x, R)                  . . .               

Electrons in fixed back-        ˆ ˆ ˆ           ˆ ˆ
Te  Vee (r )  Ven (r, R)   (x, R)    (R) (x, R)
ground nuclear geometry R                                
Nuclei on each electronic
ˆ       ˆ     ˆ ˆ
Tn   (R)  Vnn (R)   (R)  E  (R)
hypersurface                                        
Born-Oppenheimer (electronic) surfaces and nuclear energy levels

 n ( R)

ε2

ε1 (no bound levels)

Vibrations
Rotations
ε0

R
Classical nuclear motion on an electronic surface

Consider the ground-state electronic surface    0 (R )

Nuclear Hamiltonian:

ˆ ˆ          ˆ         ˆ
H  Tn   0 (R)  Vnn (R)

“Demote” to a classical Hamiltonian:
N
PI2
H ( P, R )              0 (R )  Vnn (R )
I 1 2 M I

Nuclear motion now given by Hamilton’s equations:

H                            H
RI                         PI  
PI                           R I
Classical nuclei
(R,P)

Quantum electrons
0 (x, R )
Hellman-Feynman Theorem
Ground-state electronic surface as expectation value:

ˆ
 0 (R)  0 (R) H (e) (R) 0 (R)                         ˆ            ˆ
H (e) (R)  Te  Vee (r)  Ven (r, R)
ˆ         ˆ

 0                 ˆ
H (e)                  0                                              0
 0 ( R )             0 ( R )             ˆ                            ˆ
H (e) (R ) 0 (R )  0 (R ) H (e) (R )
R I                R I                   R I                                             R I

 0                   ˆ
H (e)                            0                            0 
 0 ( R )                0 ( R )      0 (R )          0 ( R )  0 ( R )         
R I                  R I                             R I                           R I 

0

Because       0 ( R ) 0 ( R )  1                    0 ( R ) 0 ( R )  0
R I
Kohn-Sham density functional theory

Except for very small systems, we cannot solve for the exact    0 (x1 ,..., x M )
Density functional theory represents a compromise between accuracy and computational cost.

Wave function ansatz:
 1 (x1 )        1 (x M )
0 (x1 ,..., x M ) 
 M (x1 )        M (x M )

Single-particle orbitals:     i ( x)       i  j   ij

Electron density:
M        /2
n(r )  M              dr             drM 0 (x1 ,..., x M )          
2
 i ( x)
2
2
s1     sM                                            i 1 s z  / 2
Kohn-Sham density functional theory
Total energy functional:
             
E[{ }, R]  Ts [{ }]  EH [n]  Exc [n]  Eext [n, R]

Energy definitions:
1                                         1          n(r ) n(r)
Ts [{ }]     i  2  i                 E H [ n ]   dr dr 
2 i                                       2            r r'
n(r )
Eext [ n, R ]   Z I  dr
I                 r  RI
Ground-state energy via constrained minimization

                                  
{ }

 0 ( R)  min  E[{ }, R]   ij  i  j   ij           
             i, j                 

Kohn-Sham equations ( i are eigenvalues of ij )
 1 2                                                              
   VKS (r )  i (r )   i i (r)
 2                                              VKS (r)                    EH  Exc  Eext 
                                                                 n(r )
The Born-Oppenheimer Algorithm

Electrons

Nuclei

short time Δt with F

e.g. Verlet:
t 2
R I (t )  2R I (0)  R I (t )       FI (0)
MI
The Car-Parrinello scheme
Avoid explicit minimization with a fictitious adiabatic dynamics for electronic orbitals:

Lagrangian (note μ not a mass! It has units of energy x time2):

                   
M
1
L     i  i   M I R 2  E   , R    ij  i  j   ij
I              
i 1       2 I                          i, j

Equations of motion:
E                                                      E
 i          ij  j                             M I RI  
 i   j                                                R I
Conditions: 1) “Near” Born-Oppenheimer
M
1 N
  i  i
i 1
 I
2 I 1
M I R2                0
“Seed” the CP equations of motion with initially minimized orbitals.
Energy Conservation in Born-Oppenheimer and Car-Parrinello dynamics

CP 5 a.u
CP 10 a.u.
BO 10-6, 10 a.u.

BO 10-6, 100 a.u.
CP 10 a.u.

BO 10-5, 100 a.u.

BO 10-6, 100 a.u.
CP 10 a.u.

BO 10-4, 100 a.u.

System: 8 Silicon atoms     Marx and Hutter, Modern Methods and Algorithms of Quantum Chemistry
(NIC Series) 1, J. Grotendorst, ed. (Forschungszentrum, Jülich, 2000)
Energy conservation timing comparison
Marx and Hutter, Modern Methods and Algorithms of Quantum Chemistry
(NIC Series) 1, J. Grotendorst, ed. (Forschungszentrum, Jülich, 2000)

System: 8 Silicon atoms
Consider a simple 2 degree-of-freedom system:

p                                          pR
                                            R
m                                          mR

p  F ( R, )  Fbath ( p ,  bath, )    pR  FR ( R, )  Fbath ( pR ,  bath, R )

 bath,  G ( bath, , p , T )           bath, R  G ( bath, R , pR , TR )


m         mR                   T        TR
R
Analysis of the dynamics
Full phase-space vector       ( p , pR , , R, bath, , bath, R )   evolves according to

  iL
Liouville operator:
p                    pR                
iL         F ( R, )           FR ( R, )      iLbath, (T )  iLbath, R (TR )
m               p mR R              pR

Subdivision of Liouville operator:
p 
iLref ,           iLbath, (T )
m 
pR 
iLref , R           iLbath, R (TR )
mR R
                
iL  iLref ,     F ( R, )      FR ( R, )
p              pR

iL  iL  iL ref ,R
Analysis of dynamics (cont’d)
Evolution of phase space over a time Δt characteristic of nuclear motion:

(t )  eiLt (0)

Trotter factorization:

e   iLt
e
iL t / 2 iLref,R t iL t / 2
e           e                    O t          3

Exact Trotter theorem:

n
iL t / 2
    t
FR
        t
F

iLref ,
t       t
F

4 n p
t
FR


4 n pR        4 n pR                                         4 n pR
e                 lim e               e             e              2n
e             e             
n 

                                                                             

Evolution of momentum:
t / 2
pR (t / 2)  pR (0)                     dt FR ( R(0), (t ))
0
Analysis of dynamics (cont’d)

Time-average equated to phase-space average:                 (  1/ kBT             R  1/ kBTR )
  V ( R , )
2 t / 2                        d FR ( R, )e                             1
0 dt FR ( R, (t ))                       V ( R ,)

R 
Z ( ,  )
t                                   d e
Partition function for slow variable:
 R / 
Q(  R ,  )   dpR dR e          R pR / 2 mR
 Z ( R,  ) 
2

              

1                                  1
A( R)          ln Z ( R,  )                 ln P( R)
                             R
Annealing property: T  0,   
        pR2

Q(  R ; )   dpR dR exp    R       min V ( R, )  

        2mR                 
TR  10T
mR  5m
Model Problem:

                     1
2
V ( R, )  D0 R  a
2     2
  2   R
2

vR (0)vR (t )                     TR  10T
mR  300m

v (0)v (t )
Methods: Plane-wave basis sets (periodic box, FFTs)

1 ik r                    2 n           1
 ik (r )    e  cik,g eig r      g                  | g |2  Ecut
V      g                   L             2

1                                    1
n(r )   n (g)eig r                          | g |2  4 Ecut
V g                                  2

orbitals                     density

E                                    E
Car-Parrinello
 c   k*   ij cik,g
k
M I RI  
ci ,g                                R I
i ,g
j
N
ZI
l=0                                                 Eliminating      Vext (r )  
ˆ
core electrons               I 1   r  RI
l=1                                                                            M
Eext [n]    i Vext  i
ˆ
i 1
l=2
n(r )
  Z I  dr
I          r  RI
N           l
Vpseud    vl  r  R I
ˆ                                          lm         lm
I 1 l  0 m  l


                                                 lm
N            l
   vl  r  R I   vl  r  R I   vl  r  R I                                lm
I 1 l  0 m  l

N l 1

                            lm
N                                        l
  vl  r  R I     vl  r  R I   vl  r  R I                                   lm
I 1                        I 1 l  0 m  l

M                              N
E pseud [n,{ }]    i V pseud  i    dr n(r)vl  r  R I   ENL [{ }, R]
ˆ
i 1                           I 1
Why a real-space basis?

• Plane-waves are elegant but scale as N 2M

• Slow convergence of plane waves to the basis set limit.

• Ease of localizing orbitals.

• Ease of representing position-dependent operators.

• Exact representation of 
2

• Common choice – Gaussians

                                                         |r  R I | / 2
2    2

 i (r)                i
C G (r; R I )   G (r; R I )  N x y z e
 
, , ,I
Selecting a real-space basis (why not Gaussians?)

• Retain simplicity of plane waves.

• Systematic convergence to the basis-set limit.

• Spatially localized for possible linear-scaling.

• Position independence and orthonormality.

• No BSSE

• For flexibility of use, seek noncompact support.

• Choice: Discrete variable representations (DVRs).

J. C. Light, et al. J. Chem. Phys. 82, 1400 (1985); Edwards, Tuckerman, Friesner, Sorensen,
J. Comp. Phys. 110, 82 (1994).R. A. Friesner, Chem. Phys. Lett. 116, 39 (1985);
Bacic and Light, Ann. Rev. Phys. Chem. 40, 469 (1989); J. T. Muckerman, Chem. Phys. Lett. 173,
200 (1990); Colbert and Miller, J. Chem. Phys. 96, 1982 (1992); Light and Carrington, Adv. Chem.
Phys. 114, 263 (2000); Littlejohn and Cargo, J. Chem. Phys. 117, 27, 37, 59 (2002); Varga, et al.
Phys. Rev. Lett. 93, 176403 (2004).
Definition of a DVR
Plane-waves (at the Γ (k=0)-point) -- momentum eigenfunctions:

1
 i (r ) 
V
g
Cg ,i eig r

Discrete-variable representations (position eigenfunctions): Begin with a set of
N square-integrable orthonormal functions φi(x)
N
ui ( x)  ai  l* ( xi )l ( x)
l 1

On an appropriately chosen quadrature grid {x1,…,xN}

 ij
ui ( x j ) 
ai
Expand orbitals as:

 i (r )       i
Clmnu l ( x)um ( y )un ( z )
l ,m,n

Y. Liu, D. Yarne and MET, PRB 68, 125110 (2003); H. –S. Lee and MET, JPCA 110, 5549 (2006)
DVR convergence for a 32 water box vs. plane-waves with TM PPs

N
1
Force measure:      F       F
2
I
N   I 1

DVR basis sets allow the complete basis set limit to be reached with the ease of plane waves
Is Exc = BLYP water overstructured?
Plane-wave basis (70-85 Ry cutoff)

Grossman, et. al. JCP 120, 300 (2004)
Pseudopotentials: Hamann (1989)
85 Ry cutoff

Mantz, et. al. JPCB 110, 3540 (2006)
Pseudopotentials: Troullier-Martins
70 Ry cutoff

292 K
Morrone and Car,                                             Gaussians: TZV2P                    318 K
PRL 101, 017801 (2008)                                       VandeVondele, et. al.
Pseudopotentials:                                            JCP 122, 014515 (2005)
Troullier-Martins
70 ry cutoff
Radial distribution functions for BLYP Water
DVR                    Grid = 753, t =60 ps
Neutron
X-ray
Ensemble: NVT, 300 K, μ = 500 au

3.5

3.0                       DZVP
DZVP+BSSE-BLYP
SCP-BLYP
2.5

2.0

gOO(R)
1.5

1.0

0.5
Grossman, et. al. JCP 120, 300 (2004)
0.0
2   2.5   3   3.5    4      4.5   5   5.5   6
r(Å)
R [Å]
H. –S. Lee and MET, JPCA 110, 549 (2006)               From Akin-Ojo, et al. JCP 129, 064108 (2008)
H. –S. Lee and MET JCP 125, 154507 (2006).
H. –S. Lee and MET JCP 126, 164501 (2007).
Neutron: Soper, et. al. JCP 106, 247 (1997)                                      When basis sets are too small!
X-ray: Hura, et. al. Chem. Phys. 113, 9140 (2000)                                from C. J. Mundy (2008)
Selected References
1. R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985)
2. D. K. Remler and P. Madden, Mol. Phys. 70, 921 (1990)
3. G. Galli and M. Parrinello in Computer Simulations in Chemical Physics
(NATO ASI Series C) 397, 261 (1993)
4. M. Parrinello, Solid State Commun. 102, 107 (1997)
5. D. Marx and J. Hutter, Modern Methods and Algorithms of Quantum Chemistry
(NIC Series) 1, J. Grotendorst, ed. (Forschungszentrum, Jülich, 2000)
6. M. E. Tuckerman, J. Phys. Condens. Matter, 14, R1297 (2002)
7. F. Krajewski and M. Parrinello, Phys. Rev. B 73, 041105 (2006)
8. T. D. Kunhe, M. Krack, F. R. Mohamed and M. Parrinello,
Phys. Rev. Lett. 98, 066401 (2007)
9. H. –S. Lee and M. E. Tuckerman, J. Phys. Chem. A 110, 5549 (2006);
J. Chem. Phys. 125, 154507 (2006); J. Chem. Phys. 126, 164501 (2007).
10. E. Bohm, et. al. IBM J. Res. Devel. 52, 159 (2008)
Ab initio molecular dynamics codes:
CPMD:        http://www.cpmd.org
CP2K:        http://cp2k.berlios.de
VASP:       http://cms.mpi.univie.ac.at/vasp
PINY_MD:     http://www.nyu.edu/PINY_MD/PINY.html
OpenAtom:    http://charm.cs.uiuc.edu/OpenAtom
NWChem:     http://www.emsl.pnl.gov/docs/nwchem/nwchem.html
SIESTA:     http://www.lrz-muenchen.de/services/software/chemie/siesta

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 31 posted: 11/18/2011 language: English pages: 35
How are you planning on using Docstoc?