Transcript PAT Directorate S&T highlights

FADFT2007 ISSP 7/24/06
First principles studies of materials under extreme
condition
Tadashi Ogitsu
Quantum Simulations Group
Lawrence Livermore National Laboratory
This work was performed under the auspices of the U.S. Dept. of Energy at the
University of California/ LLNL under contract no. W-7405-Eng-48.
Collaborators
Andrea Trave, Alfredo Correa, Jonathan DuBois, Kyle Caspersen,
Eric Schwegler, and Andrew Williamson (Physic Ventures)
Lawrence Livermore National Laboratory (theory)
Gillbert Collins, Andrew Ng, Yuan Ping
Lawrence Livermore National Laboratory (experiment)
David Prendergast
Molecular Foundry, Lawrence Berkeley Laboratory
François Gygi and Giulia Galli
University of California, Davis
Eamonn Murray and Steven Fahy
Tyndall National Institute, University College, Ilreland
David Fritz and David Reis
University of Michigan Ann Arbor (experiments)
Stanimir Bonev
Dalhousie University, Canada
Outline of my talk:
• DFT: my viewpoint
• Why we need large scale simulations?
• Phase diagram of materials under
pressure
– Temperature and pressure are
extremely high
– Equilibrium property
• Dynamical response of materials upon
ultra fast laser pulse
– Ultra fast (sub ps) time resolved
measurement
– Non-equilibrium dynamics
(electrons and ions)
– Non-adiabatic?
• Rigorous theory for the ground state,
but…
– We need approximations
(LDA/GGA, pseudopot) to apply it
to a real system
– The KS eigenvalues are not
supposed to represent electron
excitation in theory, while 104,000
papers on DFT band structure (as
of 7/11/07) are found by google
Computational cost
DFT: my viewpoint
Goal!
Tight binding
DFT
QMC
Rigorousness
– So confusing… (as of April 1989)
• Why justified for excited state?
– Huge amount of literature seem to
suggest qualitatively ok (sort of
defacto standard)
– For a certain limit, some theoretical
requirement is satisfied (eg.
Koopman’s theorem)
Good cost efficiency made
DFT popular, but need further
improvement
Why large scale simulations?
• Complex material: elemental boron (8/1/07 at 17:00)
• Finite size effects
• Canonical ensemble
– Long time scale simulations
• A simple calculation could be expensive
– Eg. 2()
• Non-equilibrium (and/or non-adiabatic) dynamics?
Phase diagram of materials under pressure:
Significance of ab-initio approach
• Phase boundaries are rich in physics
– Crossing line of Gibbs free energies of different phases
• Change in structure (static total energy)
• Potential energy surface (ion dynamics -> entropy)
• Electronic structure (direct and indirect)
• Important applications in various sub-field of physics
– Modeling of interior of planets
– Fundamental questions in condensed matter physics
– Designing a novel material
Method: melting line calculation
• Two-phase simulation method (nucleation is already introduced)
T>Tmelt
T<Tmelt
J. Mei and J. W. Davenport, Phys. Rev. B 46, 21 (1992)
A. Belonoshko, Geochim. Cosmochim. Acta 58, 4039 (1994)
J. R. Morris, C.Z. Wang, K.M. Ho, and C.T. Chan, Phys. Rev. B 49, 3109 (1994)
• Ab initio method (GP and Qbox codes by Gygi at UC Davis)
– Density Functional Theory with PBEGGA
– Planewave expansion, nonlocal pseudo potential for ions
– 432 atom cell, Ecut=45Ry, -point sampling
Ab-initio two-phase MD at P=100GPa
T=2300K: melt
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
T=2200K: solidify
LiH melting line
Ogitsu et al. PRL 91, 175502 (2003)
• LiH: simple yet its phase diagram is not
well understood
Liquid
– Ionic crystal with rocksalt structure
(B1)
• What is left?
 Tmexp = 965 K
? TmGGA = 795 K (18% lower than exp)
X
Solid (Rocksalt)
 B1 phase stable up to 100GPa (exp)
? No B2 (CsCl) phase found
- All the other alkali hydride
exhibit B1-B2 transition <
30GPa
Quantum Monte Carlo Corrections to the
DFT Melting Temperature of LiH (Tmdft=790K, Tmexp=965K)
Solid LiH (T=TM)
Liquid LiH (T=TM)
-92.2
Total Energy (au)
-91.4
-92.6
-93.0
DFT
QMC
-91.8
DFT equilibrium
volume
-92.2
QMC equilibrium
volume
-93.4
DFT
QMC
Internal Energy
Correction
-92.6
-93.0
-93.8
23.0
23.5
24.0
Simulation Cell Volume (au)
24.5
23
24
25
26
27
Simulation Cell Volume (au)
• QMC predicts corrections to the internal energy and equilibrium volume
• These equation of state corrections are larger in the solid than the liquid
• Preliminary results predict an increase in TM from 790K to 880K (exp=965K)
28
Solid/solid phase boundary: Quasi Harmonic Approximation (QHA)
Karki, Wentcovitch, Gironcoli, and Baroni, PRB 62, 14750 (2000)
• Free energy surface of phases match at the phase boundary
• Free energy surface, G(P,T), of solid can be well described by
harmonic phonon model
F (V,T) = U(V) + ZPE(V) + FH(V,T)
G(P,T) =F(V,T)+PV; P = -dF/dV

-point phonon
(a)
This
work
LO
1080
1071
TO
606
593
(a) PRB 28 3415 (1983)

1
 (V )
2
F (V ,T )   kT ln[1 exp{  (V )/kT}]
ZPE(V )  
s,q
H
s,q
s,q
LiH: NaCl phase
s,q
LiH phase diagram
Theory:
• B1-B2 boundary determined by ab initio
QHA method
• B2-liquid boundary determined by ab
initio two-phase method
Experiments:
• Low-T B1-B2 boundary is being
explored by DAC experiments (Spring-8)
• High-T B2-liquid boundary by isentrope
experiments (LLNL)
Property of LiH fluid under pressure
• Strong correlation between Li and H
dynamics
– Velocity distributions reflect the
mass difference
– Diffusion constants of Li and H
are almost the same
• Dynamical H2 (Hn) formation observed
at high temperature
– Charge state of H2 in LiH fluid is
nearly neutral
– Ionicity of LiH fluid is weakened
upon dynamical H2 formation
Ab-initio two-phase method:
Computational cost
• Two approaches successfully mapped liquid/solid phase boundaries of
materials in ab-initio level
– Two-phase: Ogitsu et al. PRL 2003
– Potential switching: Sugino and Car PRL 1995
• Which is more cost efficient?
– Two-phase method is computationally intensive, while potential
switching method demands intensive human labor (many many MD
runs on P, T and the switching parameter space)
Example with LiH:
Each two-phase simulation was roughly 2-10 ps MD run with 432 atoms cell
In total, to map out the melting line for 0-200 GPa, the CPU cost equivalent to a half year
with a linux cluster (128 cpu) was used (2002-2003)
Note: low density costs more (nature of planewave + faster dynamics at higher pressure)
Summary on LiH phase diagram
• It has been demonstrated first time that ab-initio two-phase method is
feasible
– LDA/GGA seems to underestimate the melting temperature
– QMC corrections look promissing
• B1/B2/liquid phase boundaries of LiH have been calculated for a wide
range of P, T space
• Property of compressed LiH fluid has been studied from first-principles
– Correlated Li and H dynamics
– Dynamical Hn clustering yielding weakening of ionicity
Melting line of hydrogen
Metallic hydrogen under pressure [Wigner and Huntington (1935)]
Group
I
Group
VII
Molecular hydrogen solid
Atomic hydrogen solid
Pressure
bcc
bcc
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
r
bcc
bcc
bcc structure (?) - metal
hcp structure - insulator
bcc
bcc
Large zero-point motion => possible low-T liquid state?
[Brovman, Kagan, Kholas, JETP (1972)]
Two possible scenarios [N. Ashcroft, J Phys. Cond. Matt. (2000)]
Liquid
Liquid
Insulator
Metal
Insulator
Metal
A hypothetical scenario towards the low-T liquid
H H K ~ EH 2
liquid
liquid: H

liquid: H2
d H 2 ~ d H 2 H 2
H2
solid
solid

H
H
Measured melting line of hydrogen:
Gregoryanz et al. PRL 90 175701 (2003)
?
Exp can reach the P, T range of interest
Exp could not locate the melting point above 44GPa
Ab initio melting curve supports low-T liquid scenario
Experiments:
Gregoryanz et al. PRL (2003).
Datchi et al. PRB (2000).
Diatschenko et al PRB (1985).
Theory:
Bonev, Schwegler, Ogitsu and Galli, Nature, 2004.
Non-molecular fluid
Metallic super fluid at around 400GPa?
Babaev, Subda and Ashcroft, Nature 2004
Babaev and Ashcroft, PRL 2005
Molecular
fluid
Solid
Reminder: experiment can reach to the P, T range
Could we suggest how to detect melting?
• Change in distribution comes from the
tail of MLWF spread
– Net overlap is changing at high P
• Stronger asymmetry observed in liquid
MLWFs at high P
– Suggest enhanced IR activity in
liquid
solid
liquid
MLWF spread distribution at Tm
Summary on the melting line of hydrogen
• Maximum in melting line of
hydrogen is found by ab-initio
two-phase method
• The negative slope is explained by
weakening of effective inter
molecular potential. Dissociation
of molecule is not necessary
• IR activity measurement might be
able to detect the high pressure
melting curve (given that the
condition is experimentally
accessible)
Why higher pressure phase has not been well
understood? Limit in computational approach
• Does LDA/GGA work?
– ~200GPa might be OK (Pickard and Needs, Nature Physics Jul, 2007)
– No well established reference system to compare with
• Quantum effect of proton
– DFT/path-integral (maybe DMC/path-integral) is feasible, however, within
adiabatic approx.
– Full (elec & ion) path-integral: lowest temperature record is about 5000K
• Non-adiabatic electron-phonon coupling
– Crucial if metallic
Breakthrough in computational approach needed
What is limiting high pressure experiments?
• To reach high P, T itself is challenging (diamonds break)
• Small sample
• Probe signal needs to go through diamond/gasket
– S/N ratio problem
By Russel Hemley at Carnegie Institution
• Direct structural measurement (X-ray, neutron) cannot
reach too high pressure
– X-ray cannot determine the orientation of H2 (X-ray
scatter off electrons)
• Most reliable experimental techniques, Raman/IR,
provide only indirect information to the structure
– Hidden challenge for theory: How do we know the
structure? [Pickard and Needs, Nature Phys 2007]
Dynamical response of materials upon ultra fast laser
pulse
• Advance in the pump and probe
experiments made sub pico second
time resolution possible with
– Ultra-fast Electron Diffraction
– Dielectric function measurement
– Raman/IR
– X-ray
• Time evolution of phase transition,
chemical reaction (breaking/making a
bond) can be directly measured!
• Big challenge for theory since
– Non-equilibrium
– Adiabatic approximation might be
breaking down
Time evolution of electron
diffraction of Al
At t = 0, the laser pulse (70 mJ/cm2) is
induced
Siwik et al. Science 302, 1382 (2003)
The Jupiter Laser Facility at LLNL
Probe
Pump
Transmitted Probe
Reflected Probe
Schematics of experiment
Pump laser pulse
Epulse  2.9 106 J / kg,150fs FWHM,   400nm(3.1eV)

t
T*
Probe laser pulse
R* 50nm thick free standing
gold foil
Broad band =400~800 nm
1. t=0: electrons are excited by 3.1eV photons
2. t>0: Transmission and Reflection (T*, R*) gives 2()
• Electronic states evolve (Auger, el-el and el-ph scattering)
• Atomic configuration evolves (energy dissipation from electrons)
Time evolution of 2() of 50nm Au film triggered by a
laser pulse [Ping et al, PRL 96, 25503 (2006)]
Epulse  2.9 106 J / kg,150fs FWHM,   400nm(3.1eV)
• Fine time resolution, simple and
reliable technology
• Interpretation of results is challenging
due to missing information

– Electronic states
– Atomic configurations
1.
2.
3.
For 1.2-4 ps, 2(), does not change (quasi steady state)
The inter-band transition peak at 2.5eV is present in the quasi steady state
The peak is enhanced from ambient condition

Parallel pair of bands (ll1) contribute on a peak in
2()
• Inter-band transition no-momentum
change
– Kubo-Greenwood formalism
 ki | p | kj 
( fi  fj ) (Ei  Ej   )
E

E
k, i, j
i
j
 2 ( ) ~ 
Ef
• Intra-band transition require change in
momentum
– Electron-phonon coupling
– The transportation function
(2F()) to DC conductivity and
the Drude form
7
6
5
4
3
2
1
0
1
1.5
2
2.5
3
Photon ene rgy (eV)
3.5
4
Current formalisms for 2() does not describe low
and high energy regimes seamlessly
• Inter-band transition no-momentum
change
– Kubo-Greenwood formalism
• Intra-band transition require change in
momentum
Ef
Super-cell
+ Kubo-Greenwood
– Electron-phonon
coupling
No-inelastic
el-ph scattering
– The transportation
function counted
(2F()) to DC conductivity and
the Drude form
• In a disordered system, elastic scattering becomes dominant, therefore,
Kubo-Greenwood formula is good enough
• The quasi-steady state of warm dense gold: ordered or disordered?
- The inter-band transition peak suggests presence of long range order
Procedure of ab-initio 2() calculation
Two temperature model:
Underlying assumptions:
• Electrons are in thermal equilibrium
• Heating of ions is slow (el-ph coupling of gold is small)
• Ions are also in thermal equilibrium
Kubo-Greenwood with
elevated Tel
Tel(t=0) + el-ph
Tion(t)
Ab-initio MD at T
2()
Note: TD-DFT-MD plus non-adiabatic correction might provide the direct answer
Comparison of exp and ab-initio 2()
No enhancement of inter-band peak observed in ab-initio 2()
• Missing el-ph coupling (eg. intra-band transition)?
• Thermalized electrons (Fermi distribution) incorrect?
• Inter-band peak implies long range order of lattice?
Note: single 2 calculation generate 1TB data
Summary
• Ab-initio 2() does not agree with experimental measurement
– No inter-band peak above 2eV
• There are many assumptions to be re-examined
– Electron distribution function
– Application of Kubo-Greenwood formalism to small  (Drude) regime
– Electron-phonon coupling constant upon excited electrons
How fast do electrons thermalize?
E=120J/cm2
• There seems to be a general consensus
on electronic thermalization time scale
of several hundred femto second
• Only one quantitative experimental
measurement on gold found [PRB 46,
13592 (1992)]
– Residual in high energy is not
explained
– Energy density is very small
compared to Ping’s experiments
E=300J/cm2
– Residual seems to grow as a
function of input energy
Thermalization time scale as a function of
input energy should be re-examined
Concluding remark
• Physics under extreme condition provide exciting and
challenging problems to computational physics community
– Significance of computational approach in highpressure physics has been and will be growing
– Ultra-fast pump and probe experimental technique
provide exciting new physics that challenges theory.
Novel computational approaches will be needed
 Ab-initio MD beyond BO approximation
 Seamless transport calculation formalism (elastic
and in-elastic el-ph scattering)