The present and future of computational mineral physics: A vision

1
The present and future of computational mineral
physics A vision

• Ronald E. Cohen Geophysical Laboratory
• Carnegie Institution of Washington

Thanks to Oguz Gulseren, Iris Inbar, Don Isaak,
Fred Marton, Gerd Steinle-Neumann Marcelo
Sepliarsky, David Singh, Michael Wu Supported
by NSF, ONR, DARPA, DOE, CIW, Keck
2
The vision

• Theoretical mineral physics as a primary tool
  used in conjunction with experiments.
• Real time tools for predicting expected
  properties
• Choose experiment
• Design experiment (best part of spectrum, P, T,
  etc. to study)
• Tools for comparing results with first-principles
  results
• If agree, good.
• If disagree outside accuracy of experiment and
  theory, new physics, or problem with experiment.
• Theoretical mineral physics as a tool to make
  predictions for Earth science.

Theory and experiment as equal partners.
3
The visionWhat is needed?

• Straightforward, trustworthy, user-friendly,
  well-documented methods and software
• Training and education in theoretical physics
• Computer power easily accessible and available
• Fast networks for exchange of data
• Easy to use, efficient, visualization and data
  analysis tools
• More accurate first-principles methods
• Multi-scale (time and space)
• Excited state

4
Why theory?

• Make predictions
• Provide guidance for planning experiments
• Help interpret experiments
• Understand material behavior

Why first-principles?

• Experiment generally provides dressed
  interactions
• More robust to study ground state potential and
  forces and unrenormalized frequencies, and derive
  materials properties from them
• Comparison with experiment shows whether
  fundamental understanding is achieved
• Better reliability of predictions outside the
  experimental data

5
Failures of theory

• Failure of theory can lead to new understanding
  of the material, the physics, or may be due to
  approximations.
• Failures are often more interesting than success.
• Theoretical methods are under continual
  development
• Increases in accuracy
• Applicability to more classes of materials
• Applicability to more complex problems

6
Present vs. Future of mineral physics

• Present
• Pure compounds
• Simpler structures
• point defects (rare)
• single crystals
• diffusion (rare)
• equations of state
• elasticity
• thermal properties (rare)
• phase transitions
• well-behaved materials
• mostly Fe-free oxides and silicates

• Future (some happening now)
• Complex solid solutions
• (real compositions)
• Complex structures
• extended defects
• grain boundaries
• plasticity
• mineral-mineral, mineral-fluid, and mineral-melt
  interactions
• Mott insulators (Fe2 etc.).

7
Calculation of physical properties from
first-principles
Exact theory is known Schrödingers equation
Local Density Approximation of DFT
Complicated many-body interactions in material of
interest (atom, molecule, crystal) at each point
e-
are like those of homogeneous electron gas with
same density as the density at that point.
Becomes more accurate with increasing P
8
Calculation of physical properties from
first-principles
Exact theory is known Schrödingers equation
Generalized Gradient Approximation (GGA)
Complicated many-body interactions in material of
interest (atom, molecule, crystal) at each point
e-
e-
are like those of homogeneous electron gas with
same density as the density at that point, and
imposed gradients at that point.
9
Calculation of physical properties from
first-principles
Exact theory is known Schrödingers equation
Weighted Density Approximation (WDA)
Complicated many-body interactions in material of
interest (atom, molecule, crystal) at each point
e-
e-
are like those of some model with the same
weighted density averaged over a sphere centered
at that point.
Becomes more accurate with increasing P
10
Calculation of physical properties from
first-principles
Exact theory is known Schrödingers equation
Quantum Monte Carlo

• Effectively solve the Schrödinger equation
  stochastically
• Ground state (variational MC, diffusion MC,
  Greens function MC, ..)
• Finite T Path integral MC

11
LDA is not always successful
Sometimes GGA helps.
Sometimes it doesnt.
Figure 1. Energies of quartz and stishovite
computed with GGA and LDA. The LDA gives
stishovite as the zero pressure structure, but
the GGA gives quartz properly as the ground state
structure for SiO2. From Hamann, 1996.
12
Density functional theory (DFT)

• All of the ground state properties of an
  electronic system are determined by the charge
  and spin densities.
• DFT in an exact many-body theory, but the exact
  functional is unknown. However, exact sum-rules
  are known.
• In practice the Local Density Approximation (LDA)
  or Generalized Gradient Approximations (GGA) are
  used.
• LDA--exchange-correlation interactions at each
  point in space are like those for the homogeneous
  electron gas with the same density
• GGA--includes gradients to the density without
  destroying sum rules satisfied by LDA.

Solve the Kohn-Sham equations
13
All Electron Methods
Linearized Augmented Plane Wave (LAPW) method

• no shape approximations for charge density or
  potential
• all-electron (pseudopotential also available)
• relativistic
• LDA, GGA
• straightforward to converge
• (Rkmax, k-points)
• forces
• linear response, Berrys phase

Wei and Krakauer (PRL, 1985) Singh (Planewaves,
Pseudopotentials, and the LAPW Method, Kluwer,
1994)
Pb pseudopotentials
Pseudopotential methods

• Potential from core electrons is replaced with
  fake (pseudo) potential that reproduces correct
  wave functions outside of some cut-off radius
• Basis can then be smooth
• Plane waves. FFT methods can be used.
• We have been using a mixed-basis with local
  pseudoorbitals plus a few plane waves to give a
  more efficient basis.
• Results depend on quality of pseudopotential

6s1, 6p0.5, 6d0.5
6s1, 6p1, 5d10
14
How do we know these methods work?

• Results of many studies at zero and high
  pressures agree well with experiments
• cases that dont agree are generally understood
  as problem cases
• There are a few cases that are genuine
  prediction.
• SiO2 stishovite CaCl2 transition
• Al2O3 corundum Rh2O3-II transition

15
CaCl2 transition in SiO2
Prediction A1g Raman mode in stishovite
decreases until phase transition to CaCl2
structure, then increases. Does NOT go to zero at
transition.
Prediction C11-C12 decreases until phase
transition to CaCl2 structure, then increases.
Does go to zero at transition-gt superplasticity
Cohen, R.E., In High Pressure Research in
Mineral Physics, Ise, Japan, January, 1991,
425-432.
16
Raman frequencies in SiO2 stishovite and
CaCl2-structure
Predicted transition (Cohen, 1991) was found by
Raman (Kingma et al. , Nature 1995).
exp.
LAPW
17
Simplified methods for dynamics and finite T
Gordon-Kim models

• Based on LDA
• The charge density is modeled rather than
  computed self-consistently
• The total energy is calculated approximately (the
  kinetic energy is approximated).

Self-energy
Overlap
Electrostatic

• VIB
• O2- is not stable in the free state
• Oxygen ion changed size with its environment.
• In VIB the O2- ion is stabilized with a 2
  charged Watson sphere in the quantum mechanical
  atomic calculation
• The radius of the sphere is varied for each O ion
  and the total energy minimized.

contours at 5 me-/bohr3
18
MgO thermal properties

• Quasiharmonic lattice dynamics (Isaak et al.,
  1990, 1992)
• thermoelasticity, cross derivatives
• seismic parameter dlnVs/dlnVp
• Molecular dynamics
• melting (Cohen and Gong, 1994 Cohen and Weitz,
  1998)
• thermal equation of state (Inbar and Cohen, 1995)
• thermal conductivity (Cohen, 1998)
• diffusion (Ita and Cohen, 1997, 1998)

19
Seismic parameter dlnVs/dlnVp
Isaak, Anderson, and Cohen, GRL, 1992.
High values of d ln Vs/ d ln Vp can be
explained as a pressure effect on dK/dT or ?s.
This is also true for perovskite (Wang and
Weidner, 1996). Values of 2-3 are not
unreasonable.
20
Molecular dynamics

• Integrate Newtons law Fma forward in time.
• Obtain P,V,T for equation of state (NEV and NPT
  ensembles)
• Classical
• high T, above Debye temperature

21
MgO thermal expansivity
Inbar and Cohen 1995 (VIB-MD)
Duffy and Ahrens 1993
Isaak et al 1990. (PIB-QLD)
Inbar and Cohen 1995 (VIB-MD)
22
Diffusion constants from MD and thermodynamic
integration
Zf number of ways of forming vacancy 1 Zm number
of equivalent diffusion paths 12 l jump
distance MD lattice a/?2 ? attempt
frequency FFT ?Gf free energy of vacancy
formation FE integration ?Gm free energy of
vacancy migration FE integration kB Boltzmanns
constant W solubility factor 2 T temperature
23
Diffusion constants from MD and thermodynamic
integration
Milman et al., PRL 70, 2928, 1993.

• Adiabatic Switching
• Sample H in an MD simulation as a progress
  variable is switched
• For ?Gf
• switch the perfect crystal to a harmonic Einstein
  crystal
• switch the defective crystal to a harmonic
  Einstein crystal
• the difference in free energies gives the
  formation free energy
• For ?Gm
• perform an MD simulation as one atom is pushed
  towards the vacancy

24
Diffusion constants from MD and thermodynamic
integration
formation
migration
25
Intrinsic diffusivity in MgO
Ita and Cohen, 1997, 1998
26
MgO diffusivity--test of homologous temperature
law
homologous temperature relation (fit all)
homologous temperature relation (fit g at P0)
Ita and Cohen, 1997
27
MgO diffusioncomparison with recent experiments
Experiments Van Orman, Fei, Hauri, and
Wang,Science, 2003.
Prediction Ita and Cohen, 1997 V3.1 cm3/mole
28
Energy versus rotation in MgSiO3 perovskite
Mg2, Si3.4,O1.8-
29
Pyrolite model (VIBexp.) compared with seismology
30
Elasticity versus P and T Method
Calculate Helmholtz free energy F(V,T,e). Get
elastic constants from strain energy
density ?Fstrain1/2cijkleijekl F(V,T,e)E0(V
,e)Eel(V,T,e)-TSel(V,T,e)Fvib(V,T,e) E0
Static energy (density functional
theory). Eel-TSel Electronic free energy,
Fermi Dirac distribution. Fvib Mean
field approach (particle in a
cell). Supercell of 52 atoms.
Displace wanderer along special direction.
APS 2003 - S34.01
31
Cell model calculations
For the exact partition function, one must
integrate over all atoms of all atoms in the
crystal.
a 3N-d integral!
In the particle in a cell model, instead a 3-d
integral is evaluated (assumes no correlations in
atomic motions).
a tremendous reduction in effort!
32
Thermoelasticity of Fe computed at inner core
density (LAPWMBPP_PIC)
Thermoelasticity of Fe computed at inner core
density compared with moduli obtained for the
inner core from free oscillation data give an
estimate of the temperature of the inner core.
G. Steinle-Neumann, L. Stixrude, R. E. Cohen, and
Oguz Gulseren, Nature 413, 57-60 (2001).
33
Calculated and experimental Hugoniot for hcp Fe
Wasserman, E, L. Stixrude and R. E. Cohen, Phys.
Rev. B 53, 8296-8309 (1996).
34
Ta Hugoniot
Rankin-Hugoniot equation
vary T at constant V until equation is satisfied.
35
Hugoniot temperatures for Fe
Yoo et al. (PRL, 1993)
PIC model (Wasserman et al.,1996)
Brown and McQueen model (JGR, 1986)
36
Ta sound velocities along the Hugoniot
premelting
Exp. (Brown and Shaner, 1984)
melting
Gulseren and Cohen, 2001
37
Thermoelasticity of Tantalum
400 GPa
Gülseren and Cohen, 2002 (experimental work
Walker et al., 1980)
0 GPa
38
Thermoelasticity of Tantalum (PIC)
?P100 GPa
Gülseren and Cohen, 2002 (experiment Walker et
al., 1980)
39
Finite temperatures and complex solutions

• Effective Hamiltonians (Rabe, Vanderbilt,
  Waghmare, Garcia, Krakauer)
• Include only softest modes
• Taylor series expanded around high symmetry
  (cubic) structure in strain, soft-mode
  displacements and soft-mode-strain coupling
• Simulate system at finite T using Monte Carlo or
  Molecular Dynamics
• Parameters are fit to first-principles (LDA)
  calculations

• Potential models fit to first-principles results
• (Sepliarsky, Stachiotti, Phillpot, Strieffer,
  Migoni, Cohen)
• Perform selected self-consistent computations
• Fit results to shell (or other) model
• Use model potential in MD

40
Fitting results for PbTiO3 Total energies
- Input data correspond to LDA-LAPW
calculations - Soft mode displacement and cell
distortions correspond to theoretical values
Soft mode displacement along 001 direction
41
Fitting results for PbTiO3
It is necessary to include explicitly stress in
the fitting
In LDA, the pressure at each point was estimated
through the equation of state
42
PbTiO3 Molecular dynamics simulations

• System 6 x 6 x 6 primitive cells
• 3-D periodic boundary conditions
• Nose-Hoover (N-stress-T) algorithm
• 16 ps

Zero Pressure
Cell evolution
Polarization evolution
This is mostly LDA error!
43
PbTiO3 M.D. under negative pressure
Pressure -5.0 GPa
At T 0 K - V 62.56 A3 - c/a 1.083
Cell evolution
Polarization evolution
Again, this is LDA error!
44
PMN model simulations
PbMg1/3Nb2/3O3 - 12 ordered along 111

• System 6 x 6 x 6 primitive cells
• 3-D periodic boundary conditions
• Hoover (N-v-T) algorithm
• 8 ps

T 0
45
PMN random configuration
At T 0 - Zero total polarization - Individual
cells has a net polarization
Temperature behaviour System of 1080 atoms
46
Different PMN configurations
polar non-polar non-polar non-polar non-polar
47
Why WDA?

• LDA provides very accurate predictions of
  ferroelectric and piezoelectric properties of,
  such as phonons, polarizations, crystal
  structure, elasticity.
• However these properties are extraordinarily
  sensitive to volume. LDA predicts smaller volume,
  while GGA overestimates it.
• Even at experimental volume, LDA gives incorrect
  strain (e.g. PbTiO3), and overestimates the
  dielectric constant e0.

48
Experiment V427.27 Bohr3 c/a1.063
LDA At expt. V c/a1.112 Fully relaxed
V408 (-4.6) c/a1.051
GGA At expt. V c/a1.068 Fully relaxed
V473 (11) c/a1.22
49
PbTiO3 LDA, GGA, new WDA
LDA and GGA are very poor for fully
relaxed PbTiO3. LDA gives reasonable and GGA
gives excellent results for properties at the
experimental volume. The new WDA gives a good
volume as well.
50
Summary

• Theoretical Mineral Physics is set to join
  experiments as an equal partner.
• Much development is needed, in training,
  education, codes, computational power, and
  fundamental theory.
• At the present time, theory can already
  contribute to our understanding of Earth
  materials, and provide provisional data for the
  Earth.

51
(No Transcript)
52
(No Transcript)
53
Ability to design materials computationally
Nature   402, 60 - 63 (04 Nov 1999)
Also
Nature   413, 54 - 57 (06 Sep 2001)
54
Correlated systems (LDAU and DMFT)
Nature   410, 793 - 795 (12 Apr 2001)
55
(No Transcript)
56
(No Transcript)
57
(No Transcript)
58
(No Transcript)
59
Summary

• Theoretical Mineral Physics is set to join
  experiments as an equal partner.
• Much development is needed, in training,
  education, codes, computational power, and
  fundamental theory.
• At the present time, theory can already
  contribute to our understanding of Earth
  materials, and provide provisional data for the
  Earth.