n - CINaM

Density Functional based Tight Binding
developpements et applications
aux agrégats moléculaire
+
Mathias Rapacioli
Laboratoire de Chimie et Physique Quantiques,
Université Toulouse III
• DFT + TB = DFTB
• Traitement des agrégats moléculaires
 neutres
 chargés
• Algorithme et implémentation
Density Functional Theory
- Théorèmes Hohenberg et Kohn
E = E[n]
+ Principe variationel
- Formalisme de Kohn and Sham
Systeme d’électrons non intéragissant avec la bonne densité :
E[n] = ∑
i
∆
1
ϕ i − + Vext +
2
2
∫
'
Zα Zβ
n(r ')
1
ϕ i + E xc[n(r)]+ ∑
r − r'
2 αβ Rα − Rβ
En pratique : fonctionelles approchées pour la contribution d’échange corrélation à
l’énergie.
Density Functional Theory
- Théorèmes Hohenberg et Kohn
E = E[n]
+ Principe variationel
- Formalisme de Kohn and Sham
Systeme d’électrons non intéragissant avec la bonne densité :
E[n] = ∑
i
∆
1
ϕ i − + Vext +
2
2
∫
'
Zα Zβ
n(r ')
1
ϕ i + E xc[n(r)]+ ∑
r − r'
2 αβ Rα − Rβ
DFTB = DFT + base de valence minimale LCAO
+ développement de Taylor de l’energie en n=n°+δn
+ on néglige les termes à trois centres
(Seifert et al., IJQC, 1996, Porezag et al., PR B, 1995, Elstner et al., PRB, 1998)
n = n0 + δ n
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
DFTB
Foulkes and Haydock, PRB, 1989
n = n0 + δ n
DFTB
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
∂ E[n]
∂ E[n]
E = E [ n°(r)] − ∫
n°(r) + ∫
(n°(r) + δ n(r))
∂ n(r) n°
∂ n(r) n°
Foulkes and Haydock, PRB, 1989
n = n0 + δ n
DFTB
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
∂ E[n]
∂ E[n]
E = E [ n°(r)] − ∫
n°(r) + ∫
n(r)
∂ n(r) n°
∂ n(r) n°
Foulkes and Haydock, PRB, 1989
DFTB
n = n0 + δ n
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
Foulkes and Haydock, PRB, 1989
∂ E[n]
∂ E[n]
E = E [ n°(r)] − ∫
n°(r) + ∫
n(r)
∂ n(r) n°
∂ n(r) n°
E rep (n°)
HKS(n°)
DFTB
n = n0 + δ n
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
Foulkes and Haydock, PRB, 1989
∂ E[n]
∂ E[n]
E = E [ n°(r)] − ∫
n°(r) + ∫
n(r)
∂ n(r) n°
∂ n(r) n°
E rep (n°)
E=
E rep (n0 )
HKS(n°)
+
∑ϕ
iocc
i
H (n°) ϕ i
DFTB
n = n0 + δ n
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
E=
Foulkes and Haydock, PRB, 1989
+
E rep (n0 )
∑ϕ
iocc
LCAO + principe variationel :
H (n°) Ci = εi S Ci ; ϕ i = ∑ Ciµ φ µ
µ
avec:
∇2
H = φi −
+Veff [n°] φ j
2
0
ij
H (RAB ) = φi∈A
0
ij
∇2
−
+Veff [n°] φ j∈B
2
i
H (n°) ϕ i
DFTB
n = n0 + δ n
E = E [ n°(r)] +
E=
Seifert et al., IJQC, 1996,
Porezag et al., PR B, 1995
∂E[n]
∫ ∂n(r) δn(r)
n°
+
E rep (n0 )
∑ϕ
i
H (n°) ϕ i
iocc
LCAO + principe variationel :
H (n°) Ci = εi S Ci ; ϕ i = ∑ Ciµ φ µ
µ
avec:
∇2
H = φi −
+Veff [n°] φ j
2
0
ij
Néglige les termes à 3 centres :
H ij0 (RAB ) = φi∈A −
∇
+Veff [n°A + n°B ] φ j∈B
2
2
Si n° est densité des atomes dissociés :
Paramétrisation possible par paire d’atome
des éléments de matrice
DFTB
n = n0 + δ n
E = E [ n°(r)] +
E=
Seifert et al., IJQC, 1996,
Porezag et al., PR B, 1995
∂E[n]
∫ ∂n(r) δn(r)
n°
E rep (n0 )
∑ϕ
+
i
H (n°) ϕ i
iocc
Avec un raisonnement similaire (tight binding) E rep (n0 ) peut
s’écrire comme une somme de contributions de paires
E=
∑E
A,B
rep
(RAB )
+
∑ϕ
iocc
i
H (n°) ϕ i
DFTB
n = n0 + δ n
Elstner et al., PRB, 1998
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
E=
SCC
∑E
rep
+
+
(RAB )
A,B
∑ϕ
iocc
LCAO + variational principle :
H (n°) Ci = εi S Ci ; ϕ i = ∑ Ciµ φ µ
µ
i
H (n°) ϕ i
1
*
2
∂ 2 E[n(r)]
∫ ∂n(r)∂n(r') δn(r)δn(r')
n°
DFTB
n = n0 + δ n
SCC
Elstner et al., PRB, 1998
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
+
1
*
2
∂ 2 E[n(r)]
∫ ∂n(r)∂n(r') δn(r)δn(r')
n°
δn(r) = ∑ qat F at00 (r)
at
q = charges de Mulliken
E=
∑E
+
rep (RAB )
A,B
∑ϕ
i
H (n°) ϕ i
+
iocc
Natoms
0.5 *
∑Γ
0
(n
AB
AB )qA qB
A,B
LCAO + variational principle :
H (n°) Ci = εi S Ci ; ϕ i = ∑ Ciµ φ µ
ΓAA = Paramètre de Hubbard
µ
avec
ΓAB
1
∂ 2 E xc[n°]
+
=
R ∂ n(rA )∂ n(rB ) n°
DFTB
n = n0 + δ n
Elstner et al., PRB, 1998
∂E[n]
E = E [ n°(r)] + ∫
δn(r)
∂n(r) n°
E=
∑E
SCC
rep (RAB )
A,B
LCAO + variational principle :
H (n°) Ci = εi SCi
+
+
∑ϕ
iocc
i
H (n°) ϕ i
1
*
2
+
∂ 2 E[n(r)]
∫ ∂n(r)∂n(r') δn(r)δn(r')
n°
Natoms
0.5*
∑
A,B
0
Γ AB (nAB
)qA qB
DFTB
DFTB = tight binding paramétré sur DFT
DFTB SCC (2nd ordre) = tight binding paramétré sur DFT + autocohérence
Paramétrisation non empirique ou sur un système en particulier
relativement transferable, réactivité chimique
Temps de calcul = diagonalisation de l’équation séculaire (base minimale) ☓
nbre cycles autocohérents (typiquement 10)
De nombreuses revues sur les différentes applications de la DFTB :
M. Elstner, Theor. Chem. Acc. 116, 316 (2006).
M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim,
S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
Seifert and Joswig, wiley Inter. Reviews: Comp. Molec. Science, 2, 456 2012
Quelques extensions de la DFTB
• Dépendance angulaire des fluctuations de densité
• DFTB third order (Gauss et al., JCTC, 2011, 4, 931)
• Spin polarized DFTB (Köhler et al., Chem. Phys. 309, 23, 2005) = Unrestricted DFTB
• La TDDFTB : Approximations de la DFTB appliquées à la TDDFT
- Approche de la réponse linéaire (Niehaus et al. PRB 63, 085108, 2001)
- Dynamique dans les états excités (Niehaus et al. EPJD, 2005, Mitric et
́ al., JPC A, 113,
127005, 2009)
• GW : approche DFTB pour la fonction de Green (Niehaus et al. PRA, 71, 022508 2005)
• Transport : Pecchia & Di Carlo, Rep. Prog. Phys. 67, 1497 (2004); Pecchiaet al., Nano Lett. 4, 2109
(2004), Di Carlo et al., Physica B 314, 86, 2002)
• DFT + TB = DFTB
• Traitement des agrégats moléculaires
 neutres
 chargés
• Algorithme et implémentation
Agrégats moléculaires neutres
Interaction intermoléculaires
Repulsion de Pauli
En DFTB :
Erep +
Base valence non orthogonale
Interactions de dispersion
Coulomb entre multipoles des
fragments
Absentes
Terme second ordre
(1/R * charges Mulliken)
Dimère de benzene neutre
Instable !!!
Absence de dispersion : Un problème hérité de la DFT
DFT : champ de recherche très actif en DFT
- nouvelles fonctionelles (MO6, Zhao & Thrular, Theo.Chem Acc. 2008)
- couplage short / long range separation DFT/MP2 ( Savin et al., Recent devel. DFT,
1996)
- empirical correction (Piacenza & Grimme, JACS, 2005; Goursot et al. , JCTC, 2007)
DFTB : correction empirique Natoms
(Elstner et al., JCP,
2001, Zhechkov, JCTC, 2005)
ij
E Disp =
∑
ij
− fdamp (Rij )
C6
Rij6
Plusieurs choix possible pour la fonction de coupure
Dimère de benzene neutre
Prise en compte de la polarisation des liaisons :
Charges de Mulliken remplacée par charges CM3 (Li et al. JPC 1998)
CM 3
k
q
=q
Mull
k
+
Natoms
∑
k'
Dkk' Bkk'
D = parametre empirique
B = Mayer’s Bond order
-> Amélioration du dipole moléculaires (Kalinowski et al. JPC A, 2004, Simon et al., PCCP,
12)
09)
energies
-> Les Binding
charges
CM3 sont plus proches des charges DFT Electrostatic Potential Fitted
et corrigent le potentiel coulombien à grande distance (Rapacioli et al. JCP, 130, 244304,
DFTB1 (kcal/mol)
2.97
2.41
1.47
1.45
CCSDT2 (kcal/mol)
2.84
2.62
1.53
1.53
Lee et al., JPC A,
2007
• DFT + TB = DFTB
• Traitement des agrégats moléculaires
 neutres
 chargés
• Algorithme et implémentation
Agrégats moléculaires cations
Interaction intermoléculaires neutre +
En DFTB :
Polarisation
Partiellement prise en
compte (terme second ordre)
Résonance de charge
entre différentes unités
Erreur d’auto-interaction
+
Polarisation
En DFTB : bonne description de la polarisabilité moléculaire
mauvaise description de polarisabilité atomique
Molécule d’hydrogene en présence d’une charge
QA
Hydrogen molecule
-δ
+δ
Atome d’hydrogène en présence d’une charge charge
QA
Pas de polarisation pour cet atome d’hydrogène
Polarisation atomique
Introduction empirique d’une contribution de polarisation :
E pola = −∑
c
E pola = −∑
c
α c 2
Ec
2

 

R ac R bc
damp
damp
α
f
(R
)
f
(R
)
∑ c
ac
bc
3
3
R
R
ac
bc
c

 2
αc
Rab
pol
qb f damp (Rab ) 3 = ∑ qaqbγ abpol
γ
∑
ab = −
avec
2 b
Rab
ab
E DFTB = E rep (n0 ) +
∑
iocc
-
ϕi H (n°) ϕi +
1
pol
(γ ab + γ ab
)qaqb
∑
2 ab
La polarisation atomique s’introduit “très facilement” comme un terme de second ordre
- γ
pol
ab

( Ra, Rb, Rc, Rd,...) doit être calculé “en vol”, c’est un terme à n-corps
Iftner et al., J. Chem. Phys., 140, 034301 (2014)
Dontot et al., PhD, 2014
Comportement à grande distance
Dissociation (Na-O)+
Comportement en 1/R^4 à grande distance
Oxyene polarisé par Na+
Paramètre Na-O
Tarrat et al. in prep
Iftner et al., J. Chem. Phys., 140, 034301 (2014)
Dontot et al., phD (2014)
Résonance de charge
Problème : En DFT (donc en DFTB) les effets de délocalisation de charge à grande
distance sont mal décrits (problème d’auto-interaction des fonctionnelle d’échange
corrélation)
+
E(C6H6+) + E(C6H6)
Density Functional Theory
l’auto-interaction
- Théorèmes Hohenberg et Kohn
E = E[n]
+ Principe variationel
- Formalisme de Kohn and Sham
Systeme d’électrons non intéragissant avec la bonne densité :
E[n] = ∑
i
∆
1
ϕ i − + Vext +
2
2
∫
'
Zα Zβ
n(r ')
1
ϕ i + E xc[n(r)]+ ∑
r − r'
2 αβ Rα − Rβ
Auto-interaction en principe corrigée par la fonctionnelle
d’exchange-correlation
L’approche CI-DFTB
CI-DFT développée par Wu & Van Voorhis (PRA, 72, 024502, 2005) pour des molécules
Ψ = α Ψ 1+ β Ψ 2 + γ Ψ 3+ η Ψ 4 + ….
Où ψi = configuration « charge localisée sur la molécule i »
+
Pour un dimère
+
+
=α
+β
Etape 1 : Les formes localisées obtenues avec DFTB + contrainte de localisation de
charge
( H[n] + λi Pi ) ϕ k = εkSϕ k
|Ψ1 & |Ψ 2 = dét. de Slater construits à partir des OM occupées DFTB contrainte
Etape 2 : Petite interaction de configuration (=model valence bond)
E 12: = Ψ1 | H |Ψ 2
 E E  α 
 1 S  α 
1

 E 12
÷  =ε 
E 2  β 
 S12
12
12
÷  
1  β 
avec
S12 = Ψ1 |Ψ 2
Rapacioli et al., JCTC, 2011, 7, 44
Dimère de benzène
+
- Bonne asymptote de dissociation
- EBinding = 18.9 kcal/mol ≈ Eexp = 17.0 kcal/mol (Rusyniak et al., JPC A, 107, 7656, 2003)
17.6 kcal/mol (Meot-Ner et al., JACS, 100, 5466,1997)
20.6 kcal/mol (Hiraoka et al., JCP, 95, 8413, 1991)
E ab initio : 18.34 kcal/mol (EOM-IP-CCSD, Pieniazek, JCP 2007)
Rapacioli et al., JCTC, 2011, 7, 44
Dimère d’eau
+
+
DFT : La structure la plus stable dépend clairement du choix de la fonctionnelle
Sources : 1 Lee and Kim, JCTC, 5, 976, 2009
2 Rapacioli et al., Phys. Stat. Sol. B 249, No. 2, 245–258 (2012)
5 Cheng et al., JPCA, 113, 13779, 2009
Potentiels d’ionisation agrégats pyrène
1- Structures agrégats neutres : replica exchange (parallel tempering) Monte Carlo
Potentiels d’ionisation agrégats pyrène
1- Structures agrégats neutres : replica exchange (parallel tempering) Monte Carlo
Pour l’hexamer, il existe 2 structures dégénérées (<0.01 eV)
Structure favorable pour l’ion
Structure défavorable pour l’ion
Potentiels d’ionisation agrégats pyrène
1- Structures agrégats neutres : replica exchange (parallel tempering) Monte Carlo
2- Potentiels d’ionisations verticaux
Pour l’hexamer, il existe 2 structures dégénérées (<0.01 eV)
Potentiels d’ionisation agrégats pyrène
1- Structures agrégats neutres : replica exchange (parallel tempering) Monte Carlo
2- Potentiels d’ionisations verticaux
Pour l’hexamer, il existe 2 structures dégénérées (<0.01 eV)
Données expériementales au dispositif
SOLEIL,
Projet GASPARIM (C. Joblin, P. Brechignac …)
Etats électroniques excités
Les excitations locales (négligée pour calculer l’état fondamental) necessaires pour le calcul des
états excités.
+
+
=α
+α’
+
+β
 E1

 E 12
1
E 12  α 
÷  =ε 
E 2  β 
 S12
+
** +
* +
+α’’ +α’’ +α’’
+β’
S12  α 
÷  
1  β 
+ α’’
+β’’ +α’’ +α’’
+ ….
*
** +
* +
*
+ β’’
+
+ ….


0 E12 E12'
E12'' ÷


 E1 0
1 0 0 S12 S12'
S12'' ÷ α
α 



E1* 0 E12* E12*' E12*'' ÷ α ' ÷
*
*'
*'' ÷ α '


÷

1
0
S
S
S

÷
12
12
12
**
**
**'
**''

÷

÷
 α ''

E1 E12 E12 E12 ÷ α ''
**
**'
**''

÷
1
S
S
S
=
E
÷ i

12
12
12

÷

÷ β
E2
0
0 ÷ β ÷

1
0
0 ÷

 β' ÷
 β'

÷
*
E2
0 ÷
1
0 ÷

÷


β '' 
β ''
** ÷
1 ÷

 la molécule chargée.

E
2
trou
dans
les
orbitales
locales
de


Excitations locales obtenues en transferant le
Energie recalculée en DFTB avec cette nouvelle
Dontot et al., submitted, JCTC

÷
÷
÷
÷
÷
÷
÷

Validation DFTB vs ab initio
CASPT2
CASPT2
DFTB-CI
+
+
DFTB-CI
Dontot et al., submitted, JCTC
Validation DFTB vs ab initio
DFTB-CI
CASPT2
Force d’oscillateur (au)
+
Intensité des transitions vers
le 1er et second états excités
angle
Dontot et al., submitted, JCTC
• DFT + TB = DFTB
• Traitement des agrégats moléculaires
 neutres
 chargés
• Algorithme et implémentation
Vers les très gros systèmes
•
Orbitales locales (definition de centres)
•
Produits de matrices utilisant la sparsité
•
Diagonalisation partielle de l’Hamiltonien
•
Seulement quelques étapes de diagonalisation à chaque cycle SCF
•
Combine calculs simple/double precision
•
Elements de matrices calculés en vol
•
Parallelisation OpenMP (acces memoire rapide)
•
•
On ne coupe aucune interaction
Au final diagonalisation exacte
Scemama, Renon, Rapacioli, soumis, JCTC
arXiv:1402.2880
Benchmarks sur des boites d’eau
1.5 million d’atomes
16 coeurs
CALMIP
ALTIX UV
128 coeurs
1 T RAM
Description d’un million d’atomes au niveau quantique
Efficace dès les tailles intéressantes pour les chimistes
Erreur absolue sur l’energie < 10-8 Hartree (taille <3000 molecules)
Scemama, Renon, Rapacioli, soumis, JCTC
arXiv:1402.2880
Perspective
Modélisation de l’ADN dans l’eau
Projet N. Tarrat, J.M. Escudier
Resumé
 DFTB = TB paramétrée sur la DFT avec autocohérence
 Corrections pour agrégats moléculaires :
 Dispersion : empirique
 Amélioration charges atomiques : dipole de liaison (CM3)
 Polarisation : empirique → équivalent terme ordre 2 à n-corps
 Resonance de charge : constrained DFTB + petite IC (valence bond)
•
Implémentation en matrices creuses : vers les très gros systèmes
Remerciements
- LCPQ : L. Dontot, F. Spiegelman, A. Simon, L. Oliveira, A. Scemama, J. Cuny, K.
Korchagina, C. Iftner, N. Suaud, …
- CEMES : N. Tarrat, M. Benoit, J. Morillo …
- CALMIP : N. Renon
- IRAP : C. Joblin, D. Kokkin
- ISMO : P. Parneix, A. Gamboa, C. Falvo
- LCAR : D. Lemoine
- TU Dresden : G. Seifert, T. Heine
Financements : ANR GASPARIM, ANR PARCS, LABEX NEXT, ERC Nanocosmos
Développements réalisés dans le code deMon-Nano (D. Salahub, T. Heine …)