Relativistic quantum chemical methods for large systems

1
Relativistic quantum chemical methods for large
systems

• Lucas Visscher
• department of Theoretical Chemistry
• Vrije Universiteit Amsterdam

2
Relativististic Quantum Mechanics

• 1905 STR
• Einstein E mc2
• 1926 QM
• Schrödinger equation
• 1928 RQM
• Dirac equation
• 1949 QED
• Tomonaga, Schwinger Feynman

3
The DIRAC programme
4
DIRAC developers

• Faegri / Helgaker(Oslo)
• Basis sets / Integrals
• Fleig (Düsseldorf)
• Configuration Interaction
• Jensen (Odense)
• MCSCF
• Kaldor (Tel Aviv)
• Coupled cluster theory
• Norman (Linköping)
• Molecular properties
• Saue (Strasbourg)
• Density Functional Theory
• Visscher (Amsterdam)
• Relativity electron correlation

Collaborations 2000-2003
5
General framework

• Fix the nuclei T 0 K
• Use Dirac equation without QED-corrections
• Compute electron wave function and energy
• Only the valence electrons or all
• 1, 2 or 4-components wave function
• High quality wave functions are expensive
• Challenge find out how accurate the results are
• Compute to the accuracy that is needed
• Pay attention to computational science aspects
  code optimization,portability, parallelization,
  user interfaces etc.

6
Accuracy of electronic structure calculations
Computational scaling yNx
7
Quaternion Modified Dirac equation L. Visscher
T. Saue, J. Chem. Phys. 113 (2000) 3996.
8
MO-integrals in quaternion form L. Visscher, J.
Comp. Chem. 23 (2002) 759.
9
Nuclear Quadrupole Couplings

• The coupling between the nuclear quadrupole
  moment Q and the electric field gradient at the
  nucleus q gives an energy splitting that can be
  measured with high precision in diatomic
  molecules
• Quantum chemistry gives q and can be used to
  obtain accurate values of Q or to predict and
  rationalize NQR or NMR observations

10
Iodine
11
How to reach this accuracy

• Basis sets
• Energy optimize large eventempered uncontracted
  sets
• Check basis set convergence of EFG in two
  molecules
• Choose near-complete HF basis set and smaller set
  for calculation of correlation contribution
• Computational method
• DC-HF with large sets to find near HF-limit
  molecular EFG
• DC-MP2 and/or SFDC-CCSD(T) for core correlation
  without or with minor truncation of virtual
  orbital space
• DC-CCSD(T) for valence correlation with more
  rigorous (typically at 20 au) truncation of
  virtual space
• Latest developments
• New libraries of relativistic all-electron basis
  sets at dz and tz level will facilitate
  application of 4-component methods
• Further development of efficient integral-direct
  (approximation) methods decrease time spend
  evaluating integrals involving the small
  component part of the wave function

12
The small component wavefunction

• The large component wave function resembles the
  non-relativistic wave function
• Exact relation between large and small component
  wave functions

• Small component wave function is related to the
  first derivative of large component wave function
• Prefactor damps singularity in the vicinity of
  nuclei
• The small component wave function is an
  embarrassingly local quantity !

13
Towards linear scaling

• Observation Major bottleneck lies in processing
  of (SSLL) and (SSSS) electron repulsion
  integrals
• Simple Coulombic Correction Neglect all
  (SSSS) integrals
• Accurate for most practical purposes
• Method requires an a posteriori correction based
  on neglected electronic charge
• May be inadequate for sensitive properties that
  probe the wave function around the nuclei
• (SSLL) type integrals strongly dominate
  calculation time
• 1-center approximation Neglect/approximate
  multi-center (SSXX) integrals
• Balance nuclear attraction and electron repulsion
• No a posteriori corrections necessary
• Work associated with (SSLL) type integrals is
  also reduced
• Implementation in progress

14
Computational scaling
15
Computational scaling
16
Implementation details

• Parallelization
• Most steps are dominated by AO-integral
  evaluation and are parallelized by distributing
  over the integral calculation tasks (master-slave
  algorithm)
• The CC algorithms are formulated in MO-basis and
  are parallelized by distribution of the
  transformed integrals (fixed distribution, no
  master necessary)
• Difficult aspects
• HF DFT the algorithms scale well but need
  substantial memory (due to the large Fock
  matrices on each node)
• 4-index trade-off between CPU and memory
  efficiency is difficult for high-angular momentum
  function shells
• CC communication of intermediate quantities is
  needed (synchronization steps) and slows down
  calculation on Beowulf-type architectures

17
Insertion of Pd into C-H in methane

• Part of larger study by Bickelhaupt and
  co-workers into oxidative insertion in specific
  bonds to develop selective catalysts
• Questions
• Is the reaction profile computed by DFT reliable
  ?
• Initial LANL2DZ-CCSD(T) calculations by Miquel
  Solà gave no clear picture larger basis sets
  give negative reaction barrier. Problem in ECP or
  in CC method ?
• Computational details
• Single point calculations at ZORA-BLYP geometries
• Hamiltonian-method (SF)DC-CCSD(T)
• Aug-cc-pVXZ sets for methane, Faegri basis
  correlating functions for Pd
• Include core-valence correlation (Pd 4s, 4p)
  explicitly
• Add counterpoise-correction
• Symmetry C2v and Cs

18
Hartree-Fock
B3LYP
BLYP
19
Reaction profile no counterpoise
10.0
5.0
0.0
-5.0
Relative energy (kcal/mol)
-10.0
-15.0
-20.0
-25.0
Reactants
Associative complex
Transition state
Products
Pd (TZ nof) CH4 (aug-cc-pVDZ)
Pd (TZ 1 f) CH4 (aug-cc-pVDZ)
Pd (TZ 1 f) CH4 (aug-cc-pVTZ)
Pd (TZ 4 f) CH4 (aug-cc-pVTZ)
Pd (TZ 4 f p) CH4 (aug-cc-pVTZ)
Pd (TZ 4 f p g) CH4 (aug-cc-pVTZ)
Pd (TZ nof) CH4 (aug-cc-pVDZ) BLYP
20
Reaction profile with counterpoise
20.0
15.0
10.0
5.0
Relative energy (kcal/mol)
0.0
-5.0
-10.0
-15.0
Reactants
Associative complex
Transition state
Products
Pd (TZ nof) CH4 (aug-cc-pVDZ)
Pd (TZ 1 f) CH4 (aug-cc-pVDZ)
Pd (TZ 1 f) CH4 (aug-cc-pVTZ)
Pd (TZ 4 f) CH4 (aug-cc-pVTZ)
Pd (TZ 4 f osani p) CH4 (aug-cc-pVTZ)
Pd (TZ 4 f osani pg) CH4 (aug-cc-pVTZ)
Pd (TZ nof) CH4 (aug-cc-pVDZ) BLYP
21
Some observations

• DFT gives indeed a reasonable description of the
  activation energy
• Spin-orbit effects are negligible (max 0.4
  kcal/mol)
• BSSE errors are not to be neglected !
• Relativistic CC can be used as standard tool, if
  youre patient enough.
• Bottleneck is the N7 scaling of the conventional
  CC algorithm that is magnified by a relativistic
  prefactor

22
How to treat large systems ?

• Expansion in predefined (atomic) set of 4-spinors
• Can be implemented via matrix transforms
• Should save on CPU due to the integral-direct
  techniques
• Will reduce memory requirements
• Can be extended to frozen core or 4-component
  ECPs
• Can be combined with e. g. Dyalls NESC or
  Barysz-Sadlej-Snijders infinite-order
  Douglas-Kroll schemes
• Density fitting and other linear scaling
  techniques
• Should take care of SS-type Coulomb interactions
• Can also be used in the correlated methods
• Further implementation of the spinfree algorithm
• Reduces time spent in CI and CC to
  non-relativistic level
• Makes perturbative approaches possible (molecular
  mean-field approximation)
• Embedding techniques to treat solvent effects

23
Plans for the future Combine DFT and ab initio
methods

• Embedding of wave functions by means of DFT
• Couples the efficiency of DFT to the accuracy of
  (relativistic) ab initio methods
• Should give improved description of active sites
  in large molecules (Ac in extraction ligands, Me
  centers in proteins)
• Interaction between DFT and WF region via a
  potential and a nonadditive kinetic energy
  functional
• General formalism allows for easy inclusion of
  various wave function based techniques

Molecular Mechanics
Ab Initio
Density Functional Theory
24
Actinide Chemistry

• Small molecules
• Use 1- or 2-component ZORA-DFT and 4-component ab
  initio methods.
• Use relativistic CCSD(T) calculations to validate
  the minima found by the DFT method and for
  establishing the relative stabilities of
  different structures.
• Goal is to calibrate the intrinsic accuracy of
  the DFT method in these systems in cases where
  direct comparison with experiment is difficult.
• Larger model systems
• Use ZORA-DFT in a QM/MM scheme to study more
  realistic models of (e.g.) ligands that are used
  in actinide separation schemes.
• Start by studying solvation and complexation of
  the uranyl ion.
• Key questions
• How to model solvation effects ?
• How many explicit solvent molecules are needed ?
• Can they be treated in a classical manner (QM/MM)
  ?
• How far can we get with ab initio methods ?
• Can we handle realistic systems (efficiency) ?
• How well does DFT perform for these systems ?
• Are the functionals adequate ?

25
The CUO molecule

• Change of ground state due to interaction with
  the heavier noble gasses
• DFT calculations by Zhou, Andrews, Li Bursten
  show that frequencies of CUO in argon correspond
  to a triplet ground state whereas the singlet
  ground state is found in neon matrices (see e.g.
  Zhou et al. JACS 121 (1999) 9712, Andrews et al.,
  JACS 125 (2003) 3126).
• Interesting for further methodological study
• Do all functionals give the same small energy
  difference ?
• Do spin-orbit interactions change the picture ?
• What do the ab initio methods predict ?

26
The CUO molecule

• Do other functionals also give this small energy
  difference between singlet and triplet ?

27
The CUO molecule

• Do spin-orbit interactions change the picture ?

28
The CUO molecule

• What do ab initio methods predict ?
• SO and electron correlation are important

29
The CUO molecule

• Preliminary conclusions
• Choice of functional is not very critical
• SO effects stabilize the triplet state by 8 (DFT)
  to 3 (ab initio) kcal/mol
• Origin of near-degeneracy may be more complicated
  than suggested by a scalar relativistic treatment
• Whats next ?
• Check consistency of the ab initio description
• Use MCSCF to optimize orbitals for triplet state
• Check core correlation and basis set convergence
• Check distance dependence
• Ab initio structures instead of DFT ?
• Interaction with noble gasses for SO-DFT and CCSD

30
Uranylfluoride in water

• Case study to validate the QM/MM approach
• EXAFS experiment in aqueous solution by Vallet et
  al., (Inorg. Chem., 2001, 40, 3516) indicates
  hepta- coordinated uranium

Can the solvent molecules be described at the MM
level of theory ?
31
Computational Section
QM/MM Method IMOMM (simple mechanical coupling)
QM region Hamiltonian-Method ZORA-DFT
(ADF) Functional BPW91 Basis Set TZ2P on all
atoms MM region Force-field parameters Amber
(ADF-implementation) Electrostatic interactions
with U with multipole-derived charges from the QM
calculation
32
Results and Discussion
A. Structure of the UO2F4(H2O)2- complex in the
gas and liquid phases
B. Electronic structure of the UO2F4(H2O)2-
complex in the gas-phase (HEXA vs HEPTA)

• C. Electronic structure of the UO2F4(H2O)2-
• complex in the liquid-phase (HEPTA vs HEXA)

33
Structure of the UO2F4(H2O)2- complex in the
gas phase
34
Liquid
QM/MM Partitioning Focus Computational
efficiency and accuracy
0 H2O
3 H2O
11 H2O
35
Influence of the partitioning on the structure
indicates the separation between the QM and MM
regions
36
Most efficient partitioning for structure

solvent molecules
MM
QM
37
Building-up the solvent
First shell
Second shell
38
Electronic structure of the UO2F4(H2O)2-
Complex in the gas-phase
Why do we not find a seven-fold coordination ?
Fragment Decomposition Scheme
EAB EA EB DEINT
DEINT DEel DEpauli DEoi
T.Ziegler and A.Rauk, Inorg. Chem., 1979, 18,
1558
39
Fragment A UO2F42- Fragment B H2O
DEint
Pauli
Or.int.
Total
El.int.
Left bar HEXA Right bar HEPTA
40
Electrostatic interaction
41
Electronic structure of the UO2F4(H2O)2-
Complex in the liquid-phase
Minimal description of first shell 3 H2O
Gap (HEXA-HEPTA) 25.9 (QM/MM) kcal/mol
23.3 (QM)
42
First shell 11 H2O
Gap (HEXA-HEPTA) 7.0 (QM/MM) kcal/mol
10.7 (QM)
43
First and partial second shell 16 H2O
Gap (HEXA-HEPTA) 5.3 (QM/MM) kcal/mol
3.8 (QM)
44
Electronic structure of the UO2F4(H2O)2-
Complex in the liquid-phase
Energy trend of the gap DEHEXA-EHEPTA
45
Fragment Decomposition Scheme
29
16
11
3
DEHEXA-EHEPTA
46
Decomposition of the bonding energy
47
Charge transfer to solvent
48
Uranyl embedding conclusions and perspectives
B3LYP and BPW91 give similar results
QM/MM gives good structures but single point QM
calculation recommended for quantitative results
HEXA and HEPTA structure energy gap is inverted
by favourable interaction with solvent molecules
going from gas to liquid phase
Next step run MD calculations to study entropic
effects
49
Acknowledgements

• Coworkers Collaborations
• Iodine NQM calculations
• Joost van Stralen (VUA)
• CUO uranyl embedding
• Ivan Infante (VUA)
• Pd insertion
• Theodoor de Jong Matthias Bickelhaupt (VUA)
• Development of DIRAC
• Trond Saue (Strasbourg), Hans-Joergen Jensen
  (Odense),
• Uzi Kaldor (Tel Aviv), Timo Fleig (Düsseldorff),
  Knut Faegri (Oslo),
• Trygve Helgaker (Oslo), Patrick Norman
  (Linköping)
• Financial other support
• Netherlands Organization for Scientific Research
• COST-D23 (Dirac) and COST-D26 (NQMs)
• National Computing Facilities (computer time)
• Pacific Northwest National Laboratories (computer
  time)