Properties of minerals at high pressure

First principles computation of thermo-elastic properties of crystals

Transcript

1. Ab initio computational methods for the prediction of stability, structure and properties of minerals at high pressure and high temperature conditions Mauro Prencipe Dip. Scienze Mineralogiche e Petrologiche – Università di Torino Ferrara 13-15 Settembre 2010. 89° Congresso SIMP: L'evoluzione del sistema Terra: dagli atomi ai vulcani La Terra profonda attraverso esperimenti e modelli teorici
2. SIMP 2010 - Ferrara 13-15 Settembre 2010 Typical pressures > 130 GPa Temperatures > 3000 K Chemistry: Mg, Si, O, Fe, Ca, Al Minerals: MgSiO3 , MgO... The Environment
3. SIMP 2010 - Ferrara 13-15 Settembre 2010 Some Motivations for Theoretical Simulations •Validation of models for the Earth’s lower mantle; •Validation of experimental results concerning structure and properties of mineral phases, at high pressure and temperature conditions; •Prediction of new mineral phases and phase transitions; •Interpretation of properties and behaviour of solids, to changes of P/T conditions, in terms of fundamental physics.
4. SIMP 2010 - Ferrara 13-15 Settembre 2010 Methods of Quantum Chemistry: a Very Brief History •Discovery of the inadequacy of classical physics to explain experimental phenomena at the atomic level; formulation of quantum physics and exact application to simple models; •Development of approximated methods of solution for multielectronic systems, to be applied to simple molecules (Valence Bond method: VB; Hartree-Fock: HF; perturbative post Hartree-Fock: MP2). Density Functional Theory: DFT. Molecules Development of highly sophisticated and accurate methods of solution: CI, MCSCF, GVB, MPn, CC, MBPT (DFT) Solids (crystals, or large molecular clusters as models of crystals) Mainly DFT; HF; hybrid HF/DFT techniques; MP2 Time 1900 Today Calculations by hands... Computational Resources Mainframes; cluster of processors
5. SIMP 2010 - Ferrara 13-15 Settembre 2010 The Fundamental Equation ) , ( ) , ( R x E R x H    Hamiltonian Operator Wave Function, depending upon the electronic spatial and spin coordinates (x) and nuclear (R) Energy of the system ) , ( s r x  r: space coordinates; s: spin coordinates of all the electrons The (time independent) Schrödinger equation:
6. SIMP 2010 - Ferrara 13-15 Settembre 2010                  k i i j i h k h hk k h ij ik k n i i k N k k r e Z Z r e r e Z m M H , , , 2 2 2 1 2 2 2 1 2 2 1 2 1 2 1 2   Kinetic energy of the nuclei Kinetic energy of the electrons Electrons-nuclei potential energy Interelectronic potential energy Internuclear potential energy The Hamiltonian Operator By following the rules imposed by the Correspondence Principle, for the translation of classical dynamic variables (observables) into the correspondent quantum expressions:
7. SIMP 2010 - Ferrara 13-15 Settembre 2010           k i i j i ij ik k n i i el r e r e Z m R H , , 2 2 1 2 2 2 1 2 ) (  ) , ( ) ( ) , ( ) ( R x R E R x R H el el    The electrons-nuclei potential contains the parametric dependence of the electronic Hamiltonian, the electronic energy and wave function ψ, upon the nuclear coordinates R. A First (usually very good) Approximation: Born-Oppenheimer Factorization of the total wave function Ψ in the product of an electronic wave function ψ, which parametrically depends on the nuclear coordinates, and a nuclear wave function φ: ) ( ) , ( ) , ( R R x R x      ψ is the solution of the electronic wave equation
8. SIMP 2010 - Ferrara 13-15 Settembre 2010 A First (usually very good) Approximation: Born-Oppenheimer     h k h hk k h el r e Z Z R E R E , 2 2 1 ) ( ) ( ) (R E NN The total energy, at a fixed nuclear configuration R, is given by the sum of the electronic energy and of the internuclear potential ENN (R). The E(R) function defines the Born-Oppenheimer surface, whose minima (R0 ) identify stable or metastable phases (at least at T=0K). ) ( ) ( ) ( ) ( 2 2 2 R E R R E R E NN el k k                The nuclear motion, in the effective electronic potential Eel (R), could be determined by solving a Schrödinger equation. However, it is generally solved within the harmonic approximation, in a small neighbourhood of R0
9. SIMP 2010 - Ferrara 13-15 Settembre 2010 Static Pressure: No Thermal and Zero Point Effects Static pressures (Pst ) at given cell volumes (V) can easily be calculated as first derivatives, with respect to V, of the static energy curve Est (V) V V E V P st st d ) ( d ) (   Beryl (Be6 Al4 Si12 O36 ) B3LYP Hamiltonian Pst (GPa) Est (hartree) V (Å3) Negative static pressure field The static bulk modulus Kst is defined as 0 | 0 d ) ( d V V st st V V P V K          where V0 is the volume at the static equilibrium (Pst =0). In the case of beryl, Kst =180 GPa
10. SIMP 2010 - Ferrara 13-15 Settembre 2010 Thermal and Zero Point Pressure: Normal Mode Vibrations j i N j i j i N q q q q E E q q q E R E        3 , 2 3 2 1 2 1 ) 0 ( ) , , , ( ) (  The effective potential which rules the nuclear motion, can be expanded in a Taylor’s series of some mass-weighted coordinates q; by truncating the series at the second order in the q’s (harmonic approximation), and by taking into account that at the equilibrium (E/qj )=0, we have:          i i i i i i i Q h Q Q H i ) ( 2 1 2 1 2 2 2   2 2 2 2 1 2 1 ) ( Q Q Q h        ) ( ) ( , , , i i m i m i i m i Q Q h     i i i i m h m      2 1 , ) 2 1 ( ,    By representing the Hamiltonian (H) in the space of the q’s, and diagonalizing, the wave function (), describing the nuclear motion, is factorized in the product of 1D-functions (normal modes)   i i i m N m Q Q Q i ) ( ) , , ( , 3 1   
11. SIMP 2010 - Ferrara 13-15 Settembre 2010 Thermal and Zero Point Pressure: the Quasi-Harmonic approximation It is assumed that the frequency of vibration of each normal mode does not depend upon the temperature of the crystal, but does depend upon the volume. For each normal mode, a Grüneisen parameter () can be defined, as: V V i i i        A whole vibrational spectrum is evaluated for a number of values of the unit cell volume in given range, and the Grüneisen’s parameters are then estimated. Beryl, B3LYP Hamiltonian. Distribution of Grüneisen’s parameters
12. SIMP 2010 - Ferrara 13-15 Settembre 2010                j j j j j j j T st T n V h V h V E T V P      ) , ( 2 ) , ( 0 Total Pressure Static Pressure (Pst ) Zero Point Pressure (Ppz ) Thermal Pressure (Pth ) Thermal and Zero Point Pressure: the Quasi-Harmonic approximation                   j kT V T j Z kT T V F V F T V P T V TS T V U T V F / ) ( e , lnZ ) , ( ) , ( ) , ( ) , ( ) , (  Given the classical relation between F and P, and the statistical definition of F in terms of the partition function Z: Derivation of F(V,T), with respect to V Number of phonons of the jth oscillator, at the temperature T
13. SIMP 2010 - Ferrara 13-15 Settembre 2010 Thermal Pressure: the Phonon Gas 1 1 ) , (   kT h e T n   Distribution of the average number of phonons, at the temperature T, according to the Bose-Einstein statistics: Contributions to the pressure of the phonon gas (thermal pressure) from the phonons of different frequencies, at three different temperatures. Case of beryl at a total pressure of 0GPa (B3LYP calculation)
14. SIMP 2010 - Ferrara 13-15 Settembre 2010 T (K) 0 100 200 300 400 500 P st , P th (GPa) -1.3 -1.2 0.0 Ppz Zero Point pressure (right P axis) Thermal pressure (left P axis) Static pressure (left P axis) 1.25 1.30 GPa Negative thermal expansion Positive thermal expansion Total Pressure and Thermal Expansion Case of Beryl (B3LYP) Thermal expansion against a zero external pressure. To keep the total pressure (P) of the phonon gas P=Pst +Pzp +Pth at the constant value of 0GPa, being Pzp almost constant in the investigated T range, Pst = -Pth At T lower than about 280K, Pst must increase with T; this is realized through a contraction of the unit cell volume: negative thermal expansion. The reverse is true at higher temperatures.
15. SIMP 2010 - Ferrara 13-15 Settembre 2010 Thermal Expansion Case of beryl (B3LYP calculation) Thermal expansion coefficient (), as a function of temperature, at three different pressures. As the total pressure increases, the frequencies of the modes having negative ’s decrease (low-frequency region of the vibrational spectrum) and therefore, for each of such mode, the correspondent number of phonons increases. This leads to enhanced negative contributions of the phonon gas (thermal) pressure to the total pressure and therefore, at high pressures, thermal expansion is negative even at relatively high temperatures.
16. SIMP 2010 - Ferrara 13-15 Settembre 2010 Compressibility Once the total pressure (P) has been obtained as a function of the cell volume (V) and of the temperature (T), the (in)compressibility can be calculated as a function of both P and T Beryl (WC1LYP calculation) P=0 Gpa, T=300 K B3LYP WC1LYP EXP V0 697.5 692.1 676.8 K0 168 180 179 K’ 3.5 4.0 3.7 P 0 3 6 9 (dKT /dT)P -0.020 -0.022 -0.028 -0.048 Slope at different pressures (Gpa/K)
17. SIMP 2010 - Ferrara 13-15 Settembre 2010 Current Works Applications of the described algorithm: • Perovskyte/post-perovskite at the lower mantle’s (D’’ layer) conditions (Torino; poster); • (Mg,Fe)O at the lower mantle’s conditions (Torino, Genova); • Carbonates: aragonite (HP/HT) and calcite (Torino); disorder in dolomite (Perugia; poster). Near future: • Test of new hybrid HF/DFT schemes to get more accurate results (Torino; Genova; Perugia); • intrinsic anharmonicity (beyond the quasi-harmonic approximation): ad-hoc calculations on aragonite and beryl (Torino); • Solids solutions. Far future: • General perturbative/variational approach to the intrinsic anharmonicity problem; • Post Hartree-Fock perturbative methods (MP2) in solids state calculations (CRYSCOR project; Torino)
18. SIMP 2010 - Ferrara 13-15 Settembre 2010 Conclusions… Gottfried Wilhelm von Leibniz (1646-1716) Calculemus!