CMD® Studies and Special Lecture:
Path Integral Simulations
Motoyuki Shiga(Japan Atomic Energy Agency)

CMD Workshop 2020
September 4, 2020, onlineSelf-Introduction:
Development and Applications of Ab Initio Simulations
Ab initio Molecular Dynamics

Path Integral Simulations

dynamics, thermodynamics

nuclear quantum effects

Nature Commun. 307, 393 (2014)

JCP. 115, 9149 (2001)
JACS 127, 11906 (2005)

Electronic Structure Theory
molecular structures, electronic properties
spin phase transition of TTTA

Angew. Chem.
48, 4190 (2009)

Science 307, 393 (2005)

QM/MM

protonated water clusters

large molecules,condensed phases
JCP 139, 033120 (2013)mechanochemicalring opening reactions

grain boundary
embrittlement of Ni

Rare Event Samplingchemical reactions
J. Phys. Chem. Lett. 9 6207 (2018)

NonadiabaticDynamics

Onsager-Machlup
path sampling

photochemistry
JACS 129, 6405 (2007)
aqueous solutions

ultrafast dynamics
of phenol blue

reaction path search algorithemSelf-Introduction:
Software Development
PIMD - an open-source software for parallel molecular simulations
Language: MPI Fortran 90 (300,000 lines)
License: Apache 2.0 (for anyone, free of charge)
URL: https://ccse.jaea.go.jp/software/PIMD/index.en.html
Manual: https://ccse.jaea.go.jp/software/PIMD/doc/manual.pdf (235 pages)
Author: M. Shiga (JAEA)
Contributors:
K. Ishimura (IMS), Y. Kawashima (RIKEN),
H. Kimizuka (Osaka Univ.), Y. Nagai (JAEA), T. Morishita (AIST),
Y. Ota (RIST), S. Ruiz-Barragan (JAEA), B. Thomsen (JAEA)
Your contribution is also welcome!

Features of PIMD
Simulation methods•Statics: Geometry optimization, phonon calculations•Reaction paths: string method, steepest descent, gentlest ascent method•

Classical mechanics: molecular dynamics (MD), hybrid Monte Carlo•

Classical statistics: replica exchange method•

Quantum statistics: path integral MD•

Semiclassical dynamics: centroid and ring polymer MD•

Nonadiabatic dynamics: surface hopping, mean field dynamics•

Free energies: metadyamics, constrained MD, mean force dynamics

Statistical ensembles• constant energy (NVE), constant temperature (NVT),
constant pressure (NPT), constant stress (NtT)

Features of PIMD
Potentials• Ab initio quantum chemistry: SMASH (HF, B3LYP, MP2)
• Classical force fields: AMBER, CHARMM, OPLS, EAM• Interface with other codes (ab initio):
CP2K, GAMESS, GAUSSIAN, MOLPRO, NTCHEM, ORCA,QUANTUM ESPRESSO,TURBOMOLE, VASP
• Interface with other codes (semiempirical): DFTB+, MOPAC• Multiscale methods: ONIOM, QM/MM• Machine learning potentials: AENET• A subroutine intended for user-defined potentials.

Boundary conditions• free boundary, periodic boundary (cubic cell, parallel-piped cell)

Features of PIMDFeatures of PIMD

Quantum Simulations
Quantum systems/phenomenaMethods

Basics of Molecular Dynamics

Electronic structures

First principles/ab initio methods

Quantum statistics, Semiclassical dynamics

Path integral simulations

Vibrational energy levels

Vibrational SCF

Quantum wave packet dynamics

DVR method, MCTDH method

Ground state of strongly correlated systems

Variational/Diffusion Monte CarloIntroduction
Nonadiabatic Dynamics
Surface hopping, mean field dynamics
to Molecular Dynamics

Ab initio path integral simulations
Electronic structure calculations – Born-Oppenheimer approximation:
Schroedinger Eq. for electron-nucleus system with electrostatic interaction

𝐑: nuclei … 8ixed𝐫, 𝐩: electrons% Ψ𝐑 𝐫 = 𝐸𝐑 Ψ𝐑 𝐫𝐻𝐑 𝐫", 𝐩Ψ𝐑 : electronic ground state
Ab initio molecular dynamics – classical trajectory:𝑀𝐑̈ = −∇𝑉 𝐑

where𝑉 𝐑 = 𝐸𝐑

Nuclear quantum effects (beyond AIMD):?, 𝐏? 𝜑" 𝐑 = 𝐸" 𝜑" 𝐑
𝐻 𝐑

where?#𝐏?, 𝐏? =?
𝐻 𝐑+𝑉 𝐑2𝑀

𝜑! : rovibrational state(However, NQE will be solved here in a different manner, via path integrals)

Nuclear Quantum Effects
ExperimentHale and Querry, Appl. Opt. 12, 555 (1973).

IR absorption spectrum of water

15000

10000

500000kT ≒ 200 cm-1 （300 K）

OH stretching

low frequency orhigh temperature

!w
>1
kT
intra-molecularvibration
HOH bending

1000
2000
3000
-1
wavenumber cm

4000

Most vibrational modes have quantum naturehigh frequency or
low temperatureNuclear Quantum Effects
Classical averages

TheoryA(R)exp(- bV )dRòA =ò exp(- bV )dR
mass-independent

Classical energy3N
H =+ V
2b

Quantum averages
Tr Aˆ exp - bHˆA =Tr exp - bHˆ

(( ())))(

mass-dependent

Quantum energyPˆ 2
H =+ Vˆ
2M(principle of equipartition)

mass-independent

Classical free energy

Quantum free energyò exp(- bV )dR
F - F = - b log
ò exp(- bV )dR-1

B

mass-dependent

BAA

mass-independent

FB − FA = −β log−1( ( ))Tr (exp (−β Ĥ ))Tr exp −β Ĥ A
B

mass-dependent

Isotope effects on thermal properties are mainly from nuclear quantum effects

Review of Quantum Mechanics
Quantum mechanics: Variables are expressed as operators
Positionr

Momentum p
Hamiltonian operator

Position operatorr̂

¶¶r
222ˆp!
¶Hˆ =+ V (rˆ ) = + V (rˆ )22m2m ¶r

ˆ = -i!
Momentum operator p

When A is a function of positionA(rˆ) r = A(r ) r
When A is a function of momentumA( pˆ ) p = A( p ) pReview of Quantum Mechanics
Closure relation for orthonormal basis functions:å f f =1iii

r-representation

ò r r dr = 1
p-representation

1p p dp = 1ò2p!
p-basis functions are plane waves

æ ipr öp r = expç ÷,è ! øæ ipr ör p = expç÷,è ! ør r¢ = d (r - r¢)

Review of Quantum Mechanics
Time-dependent Schrödinger equation∂i! Ψ = Ĥ Ψ
∂tThe formal solution isæ iHˆ t ö
÷ Y(0)Y(t ) = expçç ÷!
èøTime evolution operator

The matrix element in r-representation

æ iHˆ t ö÷rG(r , r ¢; t ) = r ¢ expçç ÷!
èø“Propagator” or “Green function”

Review of Quantum Mechanics
The formal solution

æ iHˆ t ö
÷ Y(0)Y(t ) = expçç ÷!
èøUsing the closure relation

æ iHˆ t ö÷ r r Y(0) drr ¢ Y(t ) = ò r ¢ expçç ÷!
èøThereforeY(r¢, t ) = ò G(r, r¢; t )Y(r,0)dr
G has the full information: Y is obtained as a time-space function
G(r’’,r’;t) is equal to Y(r’,t) when Y(r,0) = d(r-r’’)Path Integral Theory
For a Hamiltonian operator of single particle2ˆpHˆ = Tˆ + Vˆ =+ V (rˆ )2m

The propagator can be expressed as (Feynman, 1965)

æ¢()G r , r; t = lim çP ®¥P2m ö÷ ò dr2 ò dr3 ... ò drP exp(iS [r ] / ! )è 2pi!Dt ø

where we defined

and the actiont
Dt = ,P

r1 = r ,rP +1 = r ¢é m æ rs +1 - rs ö 2 V (rs ) + V (rs +1 ) ùS [r ] = å ê ç÷ ú Δt
2s =1 êúû
ë 2 è Dt øPPath Integral Theory

æ iHˆ t ö
÷ r
G(r , r ¢; t ) = r ¢ expçç ÷
è ! øPæ iHˆ t ö ææ iHˆ Dt ö öæ iHˆ Dt ö
æ iHˆ Dt ö
÷ = ç expç ÷ ÷ = expç ÷ ´ "´ expç ÷expçç ÷ çç÷÷ç÷ç÷!
!
!
!
èø èèøøèøèøP

æ iHˆ Dt ö÷ r
G (r , r ¢; t ) = lim r ¢ expçç P ®¥! ÷øè

Dt =tPò r r dr = 1
æ iHˆ Dt öæ iHˆ Dt ö
æ iHˆ Dt öç÷ç÷÷ rG(r , r ¢; t ) = lim ò dr2 ò dr3 ... ò drP r ¢ expç rP ´ rP expç rP -1 ´ × × × ´ r2 expçç ÷÷P ®¥! ø! ø! ÷øèèè

Path Integral Theory
æ iHˆ Dt ö
æ iVˆDt ö æ iTˆDt ö æ iVˆDt ö÷ r » r ¢ expç ÷expç ÷expç ÷ r
G(r , r ¢; Dt ) = r ¢ expçç ÷ç÷ç÷ç! ø! ø è ! ø è! ÷øèèæ 1 ö1= ç÷ ò dp p p
è 2p! øæ ipr ör p = expç÷è ! øæ ipr öp r = expç ÷
è ! øæ iVˆDt öæ iTˆDt ö æ iVˆDt ö
æ 1 ö
ç
÷
÷expç ÷ r=ç
p p expçç ÷ò dp r ¢ expç ÷÷ç÷è 2p! ø
è 2! ø
è ! ø è 2! ø
æ iDt æ p 2 ö öæ 1 öæ iV (r ¢)Dt ö ¢æ iV (r )Dt öçç÷÷ ÷ p r expç =ç÷ ò dp expç ÷ r p expçç ÷÷
2! ø!
2m2!
è 2p! øèèøèøøè

æ iDt æ p 2 ö ö
æ 1 ö æ i(V (r ¢) + V (r ))Dt öæ iöçç÷÷ ÷ expç - p(r ¢ - r )÷=ç÷expç ÷ ò dp expçç ÷2!
è 2p! ø èøè !
øè ! è 2m ø ø12

æ i æ m æ r ¢ - r ö 2 V (r ¢) + V (r ) ö öæ m ö÷Dt ÷=ç÷ expç - çç ç÷ ÷ ÷
ç ! 2 è Dt ø2è 2π!iDt øèø øè
12æ m öæ iö=ç÷ expç - ´ Lagrangian ´ Dt ÷è 2π!iDt øè !
ø2nd order Trotter expansion

Path Integral Theory
G (r , r ¢; t ) µ lim ò dr2 ò dr3 ... ò drP exp(iS (r1 , r2 ,", rP +1 ) / ! )P ®¥Sum of all possible paths

exp(iS (r1, r2 ,", rP+1 ) / !)r’= rP+1 at time texp(iS (r1, r2 ,", rP+1 ) / !)r = r1 at tme 0

Imaginary Time Path Integral Theory
Quantum Statistical Mechanics
For Canonical Ensemble
Partition function

Statistical Average(A(rˆ ) = Z -1Tr A(rˆ )exp - bHˆ( )= ∫ r exp (−β Ĥ ) r dr
Z = Tr exp −β Ĥ()( )= Z ò A(r ) r exp(- bHˆ ) r dr
= Z -1 ò r A(rˆ )exp - bHˆ r dr-1)= ∫ G r,r;-iβ ! dr

= Z -1 ò A(r )G (r , r ;-ib! )drt → -iβ !

()Propagator in imaginary time G r , r;-ib! =()r exp - bHˆ r

Imaginary Time Path Integral Theory(G r,rʹ;-iβ !

) ∝ lim ∫ dr ∫ dr ... ∫ dr exp (−βV (r ,r ,!,r ))
2P→∞

3

P

Sum of all possible paths
exp(- b Veff (r1 , r2 , ! , rP ))r=r1 at time 0

r’=r1 at time -iβ !

eff

12P

Path Integral Simulations
Imaginary time path integral theory - quantum statistical mechanics
Application to molecular systems： Chandler, Wolynes, J. Chem. Phys. 74, 4078 (1981).
Density of two particles (nuclei)“Ring polymer connected by beads”

particle 1

particle 1

particle 2

1

equivalent
for infinite
number of beads

particle 214

4

2
2

3

3

Uncertainty principle： spread of polymer

Path integral molecular dynamics： MD sampling using forces on each bead
Ab initio PIMD uses Born-Oppenheimer forces： Quantum many electron-nuclei systems
Marx, Parrinello, J. Chem. Phys. 104, 4077, (1996). CPMD
Shiga, Tachikawa, Miura, J. Chem. Phys. 115, 9149 (2001)
BOMD
Shiga, “Path Integral Simulations”, Reference Module in Chemistry, Elsevier (2018).

Path Integral Simulations
Partition function of single-particle system(Z = Tr exp - βHˆ) where2pˆ
Hˆ = Tˆ + Vˆ =+ V (rˆ )2M

The path integral representation
P2æ
æ MP ö1 P
æ P M 2
öö2Pç
÷
÷
()()
Z = lim ççdRdR...
dRexpbw
RR+VRçå
å
1ò2PPss -1s ÷ ÷ with w P º2÷òòçP ®¥ 2πb !
P s =1!bè s =1 2øøèøèis isomorph to classical partition function of interacting P -particle system（ring polymer)

Q = N ò dR1 ò dR2 ... ò dRP exp(- bVeff (R1 ,!, RP ))PM 21 P
2w P (Rs - Rs -1 ) + å V (Rs )where the effective potential is Veff (R1 , ! , RP ) = å2P s =1s =1Classical MD of NP particle system effective potential Veff
→ Classical statistics of NP particle system→ Quantum statistics of N particle system

Path Integral vs Wavefunction Calculations
Rovibrational states

Distribution ρ x = Z -1

F2Ｅ2()

Energy(H = Z-1 å Ei exp(- βEi )i)i

F1Ｅ1Ｅ０()2∑ Φi x exp −β Ei

Quantum ensemble

F0
equivalent！
Path integral MDsampling⎡⎤
Distribution ρ x = Z Tr ⎣exp −β H δ r − x ⎦()EnergyClassical MDsampling-1() (H = Z-1Tr [exp(- βH ) H ]Classical ensemble)

Bead Convergence• Inversely proportional to temperaturewater• Proportional to the highest frequency in the system• Independent of system size
• 16-32 beads at the ambient temperature for water
300 K

8P=1
P=2
P=4
P=8
P=16
P=32
P=64

6
4
2
00.8

1OH distance [Å]

1.2

1632

energydistribution
distribution

10

number of beads

Shinoda, Shiga, Phys. Rev. E 71, 041204 (2005).

Applications of Path Integral Simulations
FeH

Guanine-Cytosine
Water & Ice
PRE 71, 041204 (2005).

N2H4
Floppy molecules
JPCA 113, 1985 (2009).

H-bonds in biomoleules
CPL 625, 174 (2015).

NH4+

BeH2

Dihydrogen bonds
JCP 125, 204310 (2006).

Hydrogen in metals
JPCC 116, 23113 (2012).

Ag+ in He
Helium droplet
PCCP 19, 13798 (2017).

Proton Transfer in Water

Marx, ChemPhysChem, 7, 1848 (2006)

Nature of Protons in Water Clusters
H5O2+

H3O2-

Protonated water dimer

Deprotonated water dimer

Optimized

PIMD

Optimized

PIMD

Symmetric

FloppyUnsymmetric

FloppyNo fixed shape due to large quantum fluctuations(zero-point vibrations / tunnel effects)
Tachikawa, Shiga, J. Amer. Chem. Soc. 127, 11908 (2005).
Suzuki, Tachikawa, Shiga, J. Chem. Phys. 138, 184307 (2013).

Nature of Protons in Water Clusters
2D distributions ROO, DROO

D5O2+
D [Å]

D [Å]

D = ROH* ー RH*O

H5O2+–ROO = 2.422 Å

2.418 ÅROO [Å]ROO [Å]

D3O2–

Geometric isotope effects:
Competing effects ofzero-point vibration and tunneling

D [Å]

D [Å]

H3O2––ROO = 2.498 Å2.504 Å

ROO [Å]ROO [Å]Tachikawa, Shiga, J. Amer. Chem. Soc. 127, 11908 (2005).

Nature of Protons in Water Clusters
100 K

3–O – H*・・・O

D [Å]

H*

H–O –H ・・・O
H

ROO [Å]

D [Å]

600 K

Shielding coeff [ppm]H

NMR chemical shift

2H*（古典）H*

10-10H（古典）Classical H **Quantum H
Classical H
Quantum H

H

200
400
temperature[K] [K]
temperature

600

Proton distribution→electrons on protons→chemical shift

ROO [Å]Shiga, Suzuki, Tachikawa,J. Chem. Phys. 132, 114104 (2010).

Nature of Protons in Water Clusters
sulfuric acid clusters

ーcoordination numberproton distribution

existence ratioー #waterーー#waterー

R(Os-H)
Sugawara, Yoshikawa, Takayanagi, Shiga, Tachikawa. J. Phys. Chem. B. 115, 11486 (2009).

Liquid water at ambient condition
H2O/D2O/T2O

H2 O

D2 O

T2 O

(古典論)

(量⼦論)

(量⼦論)

(量⼦論)

Machida, Kato, Shiga, J. Chem. Phys. 148, 102324 (2018).

Liquid water at ambient condition
First-principles BOMD [BLYP-D2 functional, PAW, (H2O)64, periodic boundary]sampling
15 psclassical MD H2O

radial distributions

ClassicalQuantumExptl

HHOHOO
quantum PIMD H2OMachida, Kato, Shiga, J. Chem. Phys. 148, 102324 (2018).

ÅLiquid water at ambient condition
Evaporation energy [kcal/mol]Classical MD
Quantum PIMD
Exptl

H2O

D2O

T2O

11.9
11.2
10.0

11.9
11.5
10.4

11.9
11.8

Effects of H-bonding H O
2qHO..O

ROH

qOH..O

ROO

ROH

ROH
RClassicalOH

qHO..O

ROO

H-bond parameter correlations:ClassicalROO

cf. D2O m.p. 3.82 oC, b.p. 101.43 oC

H2O

ROH

ROH

Machida, Kato, Shiga, J. Chem. Phys. 148, 102324 (2018).

Hydrogen in Pd
diffusion paths
T siteO site

Pd
fcc lattice

Kimizuka, Ogata, Shiga, Phys. Rev. B 97, 014102 (2018)
Kimizuka, Ogata, Shiga, Phys. Rev. B 100, 024104 (2019).

Hydrogen in Pd
diffusion paths

Quantum Transition State Theory：

T siteGillan, Phys. Rev. Lett. 58, 563 (1987).

O sitefree energy calculations from ring polymer()F q ≡ − ∫ dq() (f cent r δ r cent − q(δ r cent − q)

)

PIPI

f cent : force on the ring polymer centroidsPd
fcc lattice

Kimizuka, Ogata, Shiga, Phys. Rev. B 97, 014102 (2018)
Kimizuka, Ogata, Shiga, Phys. Rev. B 100, 024104 (2019).

Hydrogen in Pd
diffusion paths
O sitePd
fcc latticeH
Quantum free energy [eV]T sitequantum classicalreaction coordinate [nm]Kimizuka, Ogata, Shiga, Phys. Rev. B 97, 014102 (2018)

Hydrogen in Pd
diffusion paths
O sitePdT
Quantum free energy [eV]T sitefcc lattice
reaction coordinate [nm]Kimizuka, Ogata, Shiga, Phys. Rev. B 97, 014102 (2018)

Hydrogen in Pd
With lattice expansion of 2.4%

T siteO sitePd

Quantum free energy [eV]diffusion paths
Hfcc lattice
reaction coordinate [nm]cf: Nanoparticle expt. by Kofu, Yamamuro, et al., Phys. Rev. B 94, 064303 (2016)
Kimizuka, Ogata, Shiga, Phys. Rev. B 97, 014102 (2018)

Hydrogen in Pd
temperature [K]diffusion paths
T siteO sitePd
fcc latticediffusion coefficient [m2/s]

H

Kimizuka, Ogata, Shiga, Phys. Rev. B 100, 024104 (2019).

HTQuantumHydrogen in Pd
temperature [K]diffusion paths
O sitePd
fcc latticeHTdiffusion coefficient [m2/s]

T sitecf: low temperature expt by Fukutani’s group
Kimizuka, Ogata, Shiga, Phys. Rev. B 100, 024104 (2019).

Hydrogen adsorption on Pt(111) surface
Full coverage

H distribution

Pt-H distribution
T: atop
F: fccT

Classical

F

Quantum

2/3 coverage
T: atop
F: fcc
H: hcpClassical

T
F,H

Quantum

Yan, Yamamoto, Shiga, Sugino, Phys. Rev. B 101, 165414 (2020).

Low Temperature Helium Droplet
He300 cluster at 0.5 K
Quantum Simulation

36 Å

Classical Simulation

fictitious solid

Suzuki, Miyazaki, Takayanagi, Shiga, Phys. Chem. Chem. Phys. 20, 26489 (2018).

Low Temperature Helium Droplet
Charge transfer

Ionization dynamics of He300
Nonadiabatic RPMD, initially 0.5 K, 0.6 ps

hneーSuzuki, Miyazaki, Takayanagi, Shiga, Phys. Chem. Chem. Phys. 20, 26489 (2018).

Acknowledgement
Water:• Prof. M. Tachikawa (Yokohama City Univ)• Dr. M. Machida (JAEA), Dr. K. Kato (Kyushu Univ)
Helium droplets:• Prof. T. Takayanagi (Saitama Univ)
Hydrogen diffusion in Pd:• Prof. H. Kimizuka (Nagoya University)• Prof. S. Ogata (Osaka University)
Hydrogen adsopption on Pt(111):• Prof. O. Sugino, Dr. L. Yan, Dr. Y. Yamamoto (Univ of Tokyo)
PIMD:• Dr. B. Thomsen, Dr. Y. Nagai (JAEA)
JSPS KAKENHI•18H05519, 18H01693, 18K05208, 16K05675

Hydrogenomics