Dynamique moléculaire ab initio en base locale :

principes et applications.
Christophe Raynaud

To cite this version:

Christophe Raynaud. Dynamique moléculaire ab initio en base locale : principes et applications..
Autre. Université Paul Sabatier - Toulouse III, 2005. Français. �NNT : �. �tel-00079109�

HAL Id: tel-00079109

[Link]
Submitted on 9 Jun 2006

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Université Paul Sabatier
Toulouse III
École Doctorale de Chimie

Thèse
présentée
pour obtenir

le grade de Docteur es Sciences

spécialité Physico-Chimie Théorique

par

Christophe Raynaud

DYNAMIQUE MOLÉCULAIRE AB INITIO EN BASE

LOCALE : PRINCIPES ET APPLICATIONS

Soutenue le 19 juillet 2005 devant la commission d'examen :

M. J. Alberto Beswick Professeur des Universités Président

Université Paul Sabatier Toulouse III

M. Martin J. Field Chercheur CEA Rapporteur

Institut de Biologie Structurale  Grenoble

M. Mark E. Tuckerman Professeur des Universités Rapporteur

New York University

M. Bernard Bigot Haut Commissaire à l'Énergie Atomique Examinateur

Professeur des Universités

CEA  École Normale Supérieure de Lyon

M. Alain Fuchs Professeur des Universités Examinateur

Université Paris Sud XI

M. Alain Vigroux Professeur des Universités Examinateur

Université Paul Sabatier Toulouse III

M. Franck Jolibois Maître de Conférences Codirecteur

Université Paul Sabatier Toulouse III

M. Laurent Maron Maître de Conférences Codirecteur

Université Paul Sabatier Toulouse III

Laboratoire de Physique Quantique UMR 5626  IRSAMC

118, route de Narbonne - 31062 Toulouse Cedex
À Léa et Valentine,

À mes parents.
 Il n'est pas nécessaire que ces hypothèses soient vraies, ou même vraisemblables.
Une chose sut : qu'elles orent des calculs conforment à l'observation.

Osiander
préface d'Osiander, éditeur de Copernic
Remerciements

M erci à toutes les personnes qui ont pu m'aider, m'encourager, me guider et me

soutenir pendant ces trois années. Citer toutes ces personnes serait une vaine
entreprise, j'en oublierai certainement. . . Merci encore de ne pas m'en tenir rigueur !

Je tiens à remercier tout d'abord les personnes qui ont accepté de juger ce travail.
Je suis honoré que M. Martin Field et M. Mark Tuckerman aient accepté de bien
vouloir rapporter cette thèse et je les en remercie. Je remercie également M. Alberto
Beswick d'avoir présidé le jury de soutenance. Mon goût pour la chimie théorique
vient notamment de certains enseignants que j'ai eu, dont M. Bernard Bigot, je lui
suis extrêmement reconnaissant d'avoir volontiers accepté de faire partie de ce jury
de thèse. Je remercie aussi M. Alain Fuchs et M. Alain Vigroux d'avoir accepté de
juger ce travail.

Si je ne devais remercier que deux personnes, ce serait sans conteste mes deux
directeurs de thèse, Franck et Laurent. J'ai beaucoup appris à leur côté, travailler avec
eux a été un réel plaisir ; leurs qualités d'encadrants ne sont plus à démontrer, et je
les remercie d'avoir été, à tout moment, présents et attentifs. Les mots me manquent
pour exprimer toute la gratitude que je leur porte, ils ont largement fait beaucoup plus
qu'être directeurs de thèse, et nalement je les remercie pour leur amitié.

Un grand merci aussi à mon colocataire de bureau, Lionel, avec qui j'ai eu moult
discussions, scientiques ou profanes, et qui m'a supporté pendant ces années ! Je te
remercie pour ton amitié et j'espère de tout c÷ur avoir l'occasion de retravailler avec
toi.
 Ces trois années m'ont aussi permis d'exercer mes premières armes en tant qu'en-
seignant. J'espère ne pas avoir traumatisé trois générations d'étudiants et je tiens ici
à remercier toutes les personnes qui ont pu m'aider et me conseiller, et avec qui j'ai
pris un grand plaisir à enseigner. Merci à Romuald, Sophie, Fabienne (B.) et (encore !)
Franck et Laurent.
Je tiens aussi à remercier les membres du Laboratoire de Physique Quantique grâce
à qui ces trois années se sont passées dans la bonne humeur. Merci Jean-Pierre, Ro-
muald, Sophie, Manu, Alex, Marie-Catherine, Daniel, Nathalie, Florent, etc. Je tiens
à remercier particulièrement Fabienne (B.) pour son amitié et là encore j'espère que
nous aurons l'occasion dans le futur de travailler ensemble. Merci enn à tous ceux que
j'oublie. . .
Enn je tiens à remercier ma famille et mes amis, pour avoir été présents et m'avoir
encouragé le long de ces trois années.
Table des matières

Introduction 1

A Dynamique Moléculaire : Généralités 5

I Dynamique Moléculaire Quantique . . . . . . . . . . . . . . . . . . . . 5
II Dynamique Moléculaire Ab Initio . . . . . . . . . . . . . . . . . . . . . 7
II.1 Dynamique moléculaire  d'Ehrenfest  . . . . . . . . . . . . . . 9
II.2 Dynamique moléculaire  Born-Oppenheimer  . . . . . . . . . 11
II.3 Dynamique moléculaire  Car-Parrinello  . . . . . . . . . . . . 13
III Intégration des Équations du Mouvement . . . . . . . . . . . . . . . . . 16
III.1 Le point de vue de Hamilton et la mécanique statistique . . . . 17
III.2 La propriété de symplecticité . . . . . . . . . . . . . . . . . . . 21
III.3 Algorithme de Verlet . . . . . . . . . . . . . . . . . . . . . . . . 22
III.4 Les autres propagateurs . . . . . . . . . . . . . . . . . . . . . . 24
III.5 L'opérateur de Liouville et l'intégration numérique . . . . . . . 26
IV Choix des Conditions Initiales . . . . . . . . . . . . . . . . . . . . . . . 30
V Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

B La Structure Électronique des Édices Moléculaires 33

I Approximation Monoélectronique  Déterminant de Slater . . . . . . . 34
II Fonctions de Base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
II.1 Base de type onde plane . . . . . . . . . . . . . . . . . . . . . . 35
II.2 Fonctions de Slater et fonctions gaussiennes . . . . . . . . . . . 37
ii TABLE DES MATIÈRES

III Méthodes Ab Initio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

III.1 Hartree-Fock  équations de Roothan et Hall . . . . . . . . . . 38
III.2 Méthodes  post Hartree-Fock  . . . . . . . . . . . . . . . . . . 40
IV La Théorie de la Fonctionnelle de la Densité . . . . . . . . . . . . . . . 41
IV.1 Théorèmes de Hohenberg et Kohn  Approche de Kohn et Sham 41
IV.2 Les diérentes fonctionnelles . . . . . . . . . . . . . . . . . . . . 43
IV.2.a Le gaz uniforme d'électrons : l'approximation de la den-
sité locale . . . . . . . . . . . . . . . . . . . . . . . . . 43
IV.2.b L'approximation du gradient généralisé . . . . . . . . . 43
IV.2.c Les fonctionnelles hybrides . . . . . . . . . . . . . . . . 44
V Les Pseudo-Potentiels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
V.1 Les pseudo-potentiels atomiques . . . . . . . . . . . . . . . . . . 45
V.2 Les potentiels eectifs de groupe . . . . . . . . . . . . . . . . . . 46
VI Approches Hybrides : Les Méthodes  QM/MM  . . . . . . . . . . . . 47
VII Les Dérivées Premières de l'Énergie . . . . . . . . . . . . . . . . . . . . 49
VIII Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

C Dynamique Moléculaire Car-Parrinello 55

I Aspects Historiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
II Quelques Aspects Techniques . . . . . . . . . . . . . . . . . . . . . . . 57
II.1 Les équations du mouvement . . . . . . . . . . . . . . . . . . . 57
II.2 Les contraintes d'orthonormalisation . . . . . . . . . . . . . . . 59
II.2.a Verlet . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
II.2.b Verlet aux vitesses . . . . . . . . . . . . . . . . . . . . 60
III Propagation Classique de la Matrice Densité . . . . . . . . . . . . . . . 62
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello . . . . . . . . 64
IV.1 Article  C. Raynaud et al., Phys. Chem. Chem. Phys., 6 (2004)
4226 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
IV.2 Pertinence de la méthode . . . . . . . . . . . . . . . . . . . . . . 73
TABLE DES MATIÈRES iii

IV.3 Application en Verlet aux vitesses . . . . . . . . . . . . . . . . . 74

V Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

D Au Delà du Microcanonique 77
I Principe de la Mécanique Statistique non Hamiltonienne . . . . . . . . 78

II Chaînes de Thermostats de Nosé-Hoover . . . . . . . . . . . . . . . . . 81

II.1 Dénition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

II.2 Distribution canonique . . . . . . . . . . . . . . . . . . . . . . . 83

II.3 L'intégration des équations du mouvement . . . . . . . . . . . . 84

III Calculs de Propriétés . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

III.1 Simulation de spectres infrarouge . . . . . . . . . . . . . . . . . 89

III.2 Article  C. Raynaud et al., Chem. Phys. Lett., 414 (2005) 161 89

III.3 Simulation de spectres RMN . . . . . . . . . . . . . . . . . . . . 96

III.4 Article  C. Raynaud et al., accepté à Chem. Phys. Chem. . . 96

III.5 Simulation de spectres UV-visible . . . . . . . . . . . . . . . . . 120

III.6 Article  C. Raynaud et al., soumis à J. Mol. Struct. . . . . . . 120

IV Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

E Le Calcul des Énergies Libres 133

I Intégration Thermodynamique . . . . . . . . . . . . . . . . . . . . . . . 136

II Ensemble Blue-Moon . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

II.1 Relation entre les moyennes restreintes et contraintes . . . . . . 137

II.2 Force thermodynamique . . . . . . . . . . . . . . . . . . . . . . 140

III Mise en Place Pratique de Contraintes Géométriques . . . . . . . . . . 141

III.1 Contraintes holonomes . . . . . . . . . . . . . . . . . . . . . . . 142

III.2 Méthode des paramètres indéterminés . . . . . . . . . . . . . . . 142

III.3 La méthode des paramètres indéterminés combinée avec l'algo-

rithme de Verlet . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

IV Applications à la Réactivité Chimique . . . . . . . . . . . . . . . . . . 148

iv TABLE DES MATIÈRES

IV.1 Article  C. Raynaud , accepté à

et al. J. Phys. Chem. A .... 149
IV.2 Article  C. Raynaud , accepté à
et al. J. Phys. Chem. A .... 165
V Estimation de l'Énergie de Point Zéro . . . . . . . . . . . . . . . . . . . 187
VI Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
Conclusion 197

Annexes 201
Systèmes à Échelles de Temps Multiples . . . . . . . . . . . . . . . . . . . . 201
Intégrateurs d'Ordres Supérieurs . . . . . . . . . . . . . . . . . . . . . . . . 205
Force Thermodynamique et Contraintes Multiples . . . . . . . . . . . . . . . 207
Dérivées de Contraintes Holonomes . . . . . . . . . . . . . . . . . . . . . . . 213
Bibliographie 221
Introduction : une histoire de temps...

Fugit irreparabile tempus1

L e présent étant connu, prévoir ce que sera l'avenir. Comprendre, quantier et

rationaliser l'inuence de paramètres observables d'un système. Tels sont les
enjeux des sciences physiques qui, à chaque échelle, associent un ensemble de lois pré-
dictives. Dans le domaine microscopique, tout système physique est supposé évoluer
dans le respect des sept postulats de la mécanique quantique. Cependant, la résolution
de l'équation maîtresse de la mécanique quantique, l'équation d'évolution, est haute-
ment non triviale.
La formulation de Schrödinger [1] de la mécanique quantique fait apparaître l'indé-
pendance au temps de l'hamiltonien du système. L'équation d'évolution peut alors être
réduite à une forme stationnaire, en oubliant la dimension temporelle. La résolution
précise de cette relation est loin d'être une tâche aisée ; elle est encore l'un des enjeux
majeurs de la chimie quantique. Le problème à n corps, cher aux enseignants de méca-
nique, voit dans la détermination de la structure électronique des édices moléculaires
son paroxysme.
Néanmoins, des modèles simples peuvent être utilisés pour comprendre et prévoir
les géométries des molécules, telles que la théorie de Lewis [2] ou bien l'approche
VSEPR [3]. Des notions comme la symétrie, l'électronégativité et le recouvrement per-
mettent de rationaliser simplement les phénomènes chimiques. La théorie des orbitales
frontières et l'approche des fragments, couronnées par le prix Nobel décerné à Roald
1 Le temps fuit irréparable. Fin d'un vers de Virgile (Géorgiques, III, 284).
2 INTRODUCTION

Homann et Kenichi Fukui en 1981 [4, 5] ont permis de comprendre de façon qualita-
tive un grand nombre de résultats expérimentaux sur la base d'analyses relativement
simples dans un langage accessible à tous les chimistes.

Le niveau quantitatif peut aussi être atteint, par des approches plus précises, telles
que les méthodes semi-empiriques, la théorie de la fonctionnelle de la densité et enn
les méthodes ab initio. Le plus souvent, l'utilisation de ces méthodes pour l'étude théo-
rique des phénomènes chimiques privilégie ainsi une approche  statique  sur un petit
nombre de molécules, voire une seule. Entre autres, ces approches passent par une ex-
ploration précise des surfaces d'énergie potentielle mises en jeu et plus particulièrement
par une caractérisation des points remarquables de ces surfaces, tels que les minima
locaux et les états de transition. La plupart des études théoriques s'intéressant aux
propriétés structurales et spectroscopiques des molécules, ou à la réactivité chimique,
ne prennent uniquement en compte que ces extrema des surfaces d'énergie potentielle.

Cependant, la grande majorité des observations expérimentales, à laquelle les pré-

dictions théoriques sont confrontées, sont constatées à partir d'échantillons de dimen-
sion macroscopique, faisant intervenir un grand nombre de molécules, sous certaines
conditions expérimentales, de température, de pression, de concentration, etc. De plus,
les temps propres des méthodes spectroscopiques sont variables ; certaines techniques
sont  lentes  par rapport aux temps caractéristiques vibrationnels moléculaires,
d'autres sont du même ordre de grandeur, voire plus rapides. Enn, certains méca-
nismes réactionnels sont essentiellement gouvernés par un eet entropique, qui est
dicilement reproductible par une étude théorique statique. Il faut garder à l'esprit
qu'un processus chimique est avant tout un processus dynamique faisant intervenir le
déplacement des molécules et en leur sein le mouvement des noyaux.

Comment tenir compte d'une partie de ces eets d'un point de vue théorique ? L'une
des réponses à cette question apparemment simple serait, entre autres, de réintroduire
la dimension temporelle, éludée lors de la résolution de l'équation de Schrödinger. À
ce jour, il existe une grande variété de méthodes de dynamique moléculaire, mais leur
INTRODUCTION 3

domaines d'application ne sont pas toujours adaptés à l'étude de processus chimiques.

Les techniques quantiques de propagation de paquet d'ondes [6] ou toute méthode
supposant la connaissance a priori de la surface d'énergie potentielle atteignent un
degré de précision très élevé, mais elles ne sont applicables qu'à des systèmes molé-
culaires de petite taille. À l'opposé, on peut distinguer des techniques de dynamique
moléculaire basées sur une approche de la surface d'énergie potentielle modélisée de
manière simpliée par un champ de force ; celles-ci ont la capacité de traiter des sys-
tèmes de grande dimension. Toutefois, la mécanique moléculaire, fondement de ces
méthodes, ne traite pas la densité électronique et il est donc impossible de décrire la
plupart des mécanismes réactionnels d'intérêt chimique. À l'interface entre ces deux ap-
proches existent les méthodes de dynamique moléculaire ab initio qui se caractérisent
par un traitement classique de la dynamique des noyaux et par des considérations
d'ordre quantique pour la partie électronique. En général, ces approches ne nécessitent
pas la connaissance a priori de la surface d'énergie potentielle car les énergies et les
gradients nucléaires sont calculés uniquement lorsque ces informations sont utiles et
nécessaires. Ces méthodes connaissent un grand engouement depuis l'apparition de
l'approche  Car-Parrinello  [7] qui permet la propagation classique simultanée des
noyaux et des orbitales utilisées pour décrire la fonction d'onde électronique.

De nos jours, les domaines d'applications privilégiés de la méthode Car-Parrinello

sont principalement la physique des solides, les interactions avec des surfaces, les agré-
gats ou les liquides. Le traitement théorique dynamique de la réactivité chimique, or-
ganique ou inorganique, ne constitue, quant à lui, qu'une activité marginale. Ceci peut
s'expliquer par le fait que la grande majorité des programmes de type Car-Parrinello
n'utilisent que des fonctions de base de type ondes planes pour décrire la structure
électronique. Bien que cette approche ait plusieurs avantages, elle est confrontée à
un grand nombre de dicultés lorsqu'il s'agit d'étudier certains types de réactions
chimiques. Notamment, l'utilisation des ondes planes interdit l'emploi de certaines mé-
thodes de chimie quantique ayant fait leurs preuves pour l'étude de molécules d'intérêt
4 INTRODUCTION

chimique, telles que les descriptions ab initio [8]. An de résoudre ce problème, des
méthodes incluant des fonctions de base gaussiennes pour décrire la structure électro-
nique sont apparues au cours de la dernière décennie [9, 10, 11]. C'est dans ce contexte
que s'inscrit ce travail de thèse.
La première partie de ce manuscrit décrit le contexte théorique général de la dyna-
mique moléculaire ab initio et le lien de ce genre d'approche avec la physique statistique.
La deuxième partie est quant à elle consacrée aux diérentes méthodes qui permettent
de décrire la structure électronique des édices moléculaires ; cette partie n'est certai-
nement pas exhaustive mais elle décrit brièvement les diérentes approches utilisées
lors de ce travail de thèse. La partie suivante est plus spéciquement consacrée à l'ap-
proche de Car et Parrinello pour la dynamique moléculaire. Notamment, la variante de
cette méthode développée durant cette thèse est précisée. Ce travail sur la dynamique
moléculaire ab initio en base locale nous amène naturellement à prendre en compte
de façon explicite les eets de température ; la quatrième partie est ainsi dévolue aux
simulations dans l'ensemble canonique. Enn, la cinquième et dernière partie de ce
manuscrit s'intéresse à l'estimation de grandeurs thermodynamiques, en particulier les
diérences d'énergie libre, pour l'analyse dynamique de la réactivité chimique.
Notons enn que ce manuscrit présente un point de vue relativement général de
la dynamique moléculaire ab initio en base locale ; ce travail est néanmoins loin d'être
exhaustif et nous verrons au cours des diérentes parties qu'il reste moult améliorations
et illustrations à apporter dans ce domaine.
A Dynamique Moléculaire :

Généralités

Dans cette partie, l'attention est portée sur l'évolution temporelle d'un système mo-
léculaire. Plus particulièrement, les méthodes considérant les noyaux comme des par-
ticules classiques sont décrites. Toutefois, l'objectif n'est pas d'énumérer de façon
exhaustive les diérentes méthodes de dynamique moléculaire mais de préciser le
contexte dans lequel s'inscrit la dynamique moléculaire ab initio.

I Dynamique Moléculaire Quantique

L
temps :
e point de départ de toute méthode de dynamique moléculaire est la mécanique

quantique non relativiste, basée sur l'équation de Schrödinger dépendante du

∂ ~ I }; t) = H Φ({~ri }, {R
~ I }; t)
i~ Φ({~ri }, {R (A.1)
∂t
où Φ désigne la fonction d'onde et H représente l'hamiltonien non relativiste :

X ~2 X ~2 X e2 X e 2 ZI X e2 ZI ZJ
H =− ∇2I − ∇2i + − +
2MI 2me |~
ri − r~j | ~ ~j | I<J |R~I − R~J |
I i i<j I,j |RI − r

(A.2)

X ~2 X ~2
=− ∇2I − ~ I })
∇2i + VN −e ({~ri }, {R (A.3)

I
2M I i
2m e
X ~2
=− ∇2I + He ({~ri }, {R~ I }) (A.4)

I
2M I
6 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

Les jeux de vecteurs {~ri } et {R

~ I } désignent respectivement les positions électroniques

et nucléaires. À partir de l'équation de Schrödinger dépendante du temps (A.1), il est

possible de dériver les équations de la dynamique moléculaire classique [12]. Dans ce
but, les contributions électroniques et nucléaires de la fonction d'onde Φ, qui dépend
à la fois des coordonnées électroniques et nucléaires, doivent être séparées. La forme la
plus simple est le produit suivant :
µ Z t ¶
i
~ I }; t) ≈ Ψ({~ri }; t)χ({R
Φ({~ri }, {R ~ I }; t) exp ′ ′
dt Ee (t ) (A.5)
~ t0

où les fonctions d'ondes électronique Ψ et nucléaire χ sont chacune normées à chaque

instant (t). Le facteur de phase Ee est introduit an de simplier par la suite les
équations nales :
Z
Ee = ~ I }; t)He Ψ({~ri }; t)χ({R
Ψ⋆ ({~ri }; t)χ⋆ ({R ~ I }; t)d~rdR
~ (A.6)

Il faut remarquer que cette approximation dière de l'approximation Born-Oppenheimer

quant à la séparation des variables lentes et rapides [13, 14]. Cette séparation conduit,
à l'aide de l'équation (A.1), aux relations suivantes :
X ~2 ½Z ¾
∂Ψ
i~ =− 2
∇ Ψ+ ⋆ ~ ~
χ ({RI }; t)VN −e χ({RI }; t)dR Ψ~ (A.7)
∂t i
2me i
X ~2 ½Z ¾
∂χ
i~ =− 2
∇I χ + ⋆
Ψ ({~ri }; t)He Ψ({~ri }; t)d~r χ (A.8)
∂t I
2M I

Ce jeu d'équations couplées est à la base de la méthode du champ auto-cohérent dé-

pendant du temps (TDSCF) introduit par Dirac dès 1930 [15, 16]. Les électrons et
les noyaux, décrits de manière quantique, évoluent dans un potentiel dépendant du
temps : les noyaux évoluent dans le potentiel des électrons et réciproquement. Ces po-
tentiels sont évalués à l'aide des valeurs moyennes quantiques de l'autre classe de degré
de liberté, i.e. le potentiel ressentit par les noyaux est le champ moyen créé par les
électrons et celui ressentit par les électrons est le champ moyen créé par les noyaux.
L'équation (A.5) conduit ainsi à une description  en champ moyen  des dynamiques
quantiques couplées des électrons et des noyaux.
II Dynamique Moléculaire Ab Initio 7

II Dynamique Moléculaire Ab Initio

La prochaine étape dans le développement des équations de la dynamique molécu-
laire classique consiste à considérer les noyaux comme des particules ponctuelles. La
fonction d'onde nucléaire correspondant à l'approche TDSCF peut être exprimée en
fonction d'un facteur d'amplitude A et d'un facteur de phase S , selon :

~ I }; t) exp( i S({R
~ I }; t) = A({R
χ({R ~ I }; t)) (A.9)
~

Dans cette représentation polaire, A et S sont réels [17]. Cette transformation de la

fonction d'onde nucléaire, conduit, après avoir séparé les parties réelle et imaginaire,
aux équations TDSCF pour les noyaux :

∂A X 1 A
+ (∇I A)(∇I S) + (∇2I S) = 0 (A.10)
∂t I
MI 2
Z
∂S X 1 X 1 ∇2 A
+ (∇I S)2 + Ψ⋆ He Ψd~r = ~2 I
(A.11)
∂t I
2M I
I
2M I A

en fonction de ces nouvelles variables A et S . La relation (A.10) pour A peut être vue
comme une équation de conservation [17] en identiant la densité nucléaire |χ|2 = A2 .
Cette équation, indépendante de ~, assure la conservation de la probabilité particulaire
en présence d'un ux. La relation (A.11) pour S , quant à elle, contient un terme
dépendant de ~. Cette contribution disparaît dans la limite classique (~ → 0) :
Z
∂S X 1
+ (∇I S)2 + Ψ⋆ He Ψd~r = 0 (A.12)
∂t I
2M I

L'équation résultante (A.12) est ainsi isomorphe aux équations classiques du mouve-
ment dans la formulation de Hamilton-Jacobi :

∂S ~ I }, {∇I S}) = 0
+ H ({R (A.13)
∂t

où H représente la fonction de Hamilton usuelle :

~ I }, {P~I }) = T ({P~I }) + V ({R

H ({R ~ I }) (A.14)
8 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

où T et V représentent respectivement les contributions cinétiques et potentielles de

l'hamiltonien et {P~I } désignent les impulsions associées aux positions nucléaires {R
~ I },

en identiant P~I = ∇I S . Les équations classiques du mouvement sont ainsi obtenues à

partir de la relation (A.12), donnant lieu à :
Z
dP~I
= −∇I Ψ⋆ He Ψd~r (A.15)
dt

soit encore :
Z
~¨ I = −∇I
MI R ~ I })
Ψ⋆ He Ψd~r = −∇I Ve ({R (A.16)

Les noyaux évoluent ainsi classiquement dans un potentiel Ve créé par les électrons.
Ce potentiel, qui dépend de la position à l'instant (t) des positions nucléaires, est le
résultat de la moyenne de l'hamiltonien électronique He sur tous les degrés de liberté
électroniques, c'est à dire la valeur moyenne quantique hΨ|He |Ψi prise avec les po-
sitions nucléaires xes à leurs positions instantanées {R
~ I (t)}. Toutefois, la fonction

d'onde nucléaire est toujours présente dans l'équation TDSCF (A.7) associée aux de-
grés de liberté électroniques. Pour être cohérent, la fonction d'onde nucléaire χ doit
être remplacée par les positions nucléaires. Dans ce cas, le passage quantique-classique
peut être facilement réalisé en remplaçant la densité nucléaire |χ|2 de l'équation (A.7)
Q
par un produit de distributions ~ I (t)) centrées sur les positions nucléaires
~I − R
δ(R
instantanées. Cela donne, par exemple, pour l'opérateur position, la valeur attendue :
Z
~ I }; t) R
χ⋆ ({R ~ I χ({R
~ I }; t) dR
~ −→ R
~ I (t) (A.17)

Cette limite classique conduit à une équation d'évolution électronique dépendante du

temps :
∂Ψ
i~ ~ I (t)})Ψ({~ri }, {R
= He ({~ri }, {R ~ I (t)}) (A.18)
∂t
où les électrons évoluent de façon auto-cohérente avec les noyaux. Il faut souligner
que maintenant, He et Ψ, dépendent paramétriquement des positions nucléaires au
temps (t), via le potentiel d'interaction électrons-noyaux VN −e . Cette approche reliant
les équations (A.16) et (A.18) est appelée dynamique moléculaire  d'Ehrenfest ,
II Dynamique Moléculaire Ab Initio 9

en l'honneur d'Ehrenfest qui, en 1927, fut le premier à poser la question suivante1 :

comment les équations de la dynamique classique newtonienne pouvent être dérivées à
partir de l'équation de Schrödinger ? Cette approche est qualiée d'hybride ou mixte,
puisque les noyaux sont décrits classiquement alors que les électrons sont traités comme
des objets quantiques [18].

II.1 Dynamique moléculaire  d'Ehrenfest 

Les jeux d'équations d'évolution couplées (A.16) et (A.18) peuvent alors être résolus
simultanément. Ce type d'approche ne nécessite donc pas la détermination préalable
de l'hypersurface d'énergie potentielle, celle-ci peut être résolue  en vol  à l'aide de
l'équation de Schrödinger dépendante du temps (A.18). Même si l'approche TDSCF
est une théorie de champ moyen pour le mouvement des noyaux dans le champ des
électrons et réciproquement, les transitions électroniques entre états sont incluses dans
la dynamique moléculaire d'Ehrenfest. La fonction d'onde électronique Ψ peut être
exprimée comme une combinaison linéaire de plusieurs états Ψk :

∞
X
~ I }; t) =
Ψ({~ri }, {R ~ I })
ck (t)Ψk ({~ri }, {R (A.19)
k

où les coecients {ck (t)} peuvent être complexes. Dans ce cas, la norme de ces coe-
cients |ck (t)|2 décrit explicitement l'évolution temporelle des populations des diérents
états électroniques Ψk . Un choix possible pour ces fonctions de base Ψk est la base
adiabatique obtenue à partir des solutions stationnaires de l'équation de Schrödinger
indépendante du temps :

~ I })Ψk ({~ri }; {R
He ({~ri }; {R ~ I }) = Ek ({R
~ I })Ψk ({~ri }; {R
~ I }) (A.20)

1 Es ist wünschenswert, dir folgende Frage möglichst elementar beantworten zu können : Welcher

Rückblick ergibt sich vom Standpunkt der Quantenmechanik auf die Newtonshen Grundgleichungen

der klassischen Mechanik ?

10 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

Les équations d'évolution correspondantes (A.16) et (A.18) peuvent alors être écrites
selon :
X X
~¨ I (t) = −
MI R |ck (t)|2 ∇I Ek − c⋆k cl (Ek − El )dk,l (A.21)
I
k k,l

X
i~ċk (t) = ck (t)Ek − i~ ~˙ I dk,l
cl (t)R (A.22)
I
I,l

I représente le terme de couplage :

où dk,l
Z
dk,l ~
I ({RI (t)}) = Ψ⋆k ∇I Ψl d~r (A.23)

Par dénition, les éléments de couplages diagonaux sont nuls dk,k

I = ~0. L'approche
d'Ehrenfest permet ainsi de décrire rigoureusement les transitions électroniques non-
adiabatiques entre diérents états électroniques Ψk et Ψl au cours de la dynamique,
tout en considérant les noyaux comme des particules classiques.
La restriction à un seul état électronique dans le développement (A.19) (dans la
plupart des cas cet état est le fondamental) conduit à :

~¨ I (t) = −∇I hΨ0 |He |Ψ0 i

MI R (A.24)

∂Ψ0
i~ = He Ψ0 (A.25)
∂t

Il faut souligner ici que la propagation de la fonction d'onde est unitaire, c'est à dire
que la fonction d'onde préserve sa norme et que le jeu d'orbitales, utilisé pour le déve-
loppement de cette dernière, reste orthonormé. La dynamique moléculaire d'Ehrenfest
est certainement la première approche de dynamique moléculaire  en vol . Toutefois,
hormis quelques exceptions [19, 20], celle-ci n'est que peu utilisée pour des simulations
de dynamique moléculaire. En eet, l'échelle de temps et par conséquent le pas de
temps utilisé pour intégrer les équations du mouvement sont gouvernés par la dyna-
mique intrinsèque des électrons. Le mouvement des électrons étant plus rapide que
celui des noyaux, le pas de temps le plus grand utilisable est celui qui permet d'in-
tégrer l'équation d'évolution électronique. Une telle approche demande alors un eort
II Dynamique Moléculaire Ab Initio 11

calculatoire important puisque les dérivées premières de l'énergie par rapport aux co-
ordonnées nucléaires  le gradient de l'énergie  doivent être calculées à chaque pas
de la simulation.

II.2 Dynamique moléculaire  Born-Oppenheimer 

Une alternative à la précédente approche permettant d'inclure la description quan-
tique des électrons dans les simulations de dynamique moléculaire consiste à résoudre,
de façon statique, le problème de la structure électronique à chaque pas de la simula-
tion pour un jeu de coordonnées nucléaires xées à leurs positions instantanées. Ainsi,
le problème de la structure électronique revient à résoudre un problème quantique
indépendant du temps. Les positions nucléaires évoluent ainsi selon la dynamique mo-
léculaire classique et la structure électronique est obtenue en résolvant l'équation de
Schrödinger indépendante du temps. La dépendance temporelle de la structure élec-
tronique devient alors une conséquence des mouvements nucléaires et n'est plus intrin-
sèque comme dans le cas de la dynamique d'Ehrenfest. L'approche résultante, dite de
 Born-Oppenheimer , est dénie selon les relations suivantes :

~¨ I = −∇I min{hΨ0 |He |Ψ0 i}

MI R (A.26)
Ψ0

He Ψ0 = E0 Ψ0 (A.27)

pour l'état fondamental Ψ0 . La diérence principale par rapport à la dynamique d'Eh-

renfest est que l'énergie électronique E0 doit être minimisée à chaque pas du mouve-
ment nucléaire. Dans l'approche d'Ehrenfest, si la fonction d'onde initiale minimise
l'énergie, alors au cours de la simulation, la fonction d'onde propagée minimise natu-
rellement l'énergie électronique. Une extension immédiate de l'approche précédente est
l'application du même schéma à un état excité Ψk , sans considérer les couplages  non
adiabatiques . En particulier, ceci conduit à négliger les termes suivants :
Z
~ I (t)}) (A.28)
′
DIkk ({R =− Ψ⋆k ∇2I Ψk′ d~r
12 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

Ces couplages apparaissent lors du développement de l'équation de Schrödinger élec-

tronique indépendante du temps. Cette approximation, qui consiste à négliger les cou-
plages dénis par les relations (A.23) et (A.28) est dite approximation de  Born-
Oppenheimer . Une variante de cette approximation qui tient compte uniquement
des termes diagonaux des couplages non adiabatiques (A.28), i.e. k = k ′ , est dite ap-
proximation de  Born-Huang . Ces couplages deviennent importants lorsque les états
électroniques k et k ′ se rapprochent, c'est-à-dire lorsqu'ils sont presque dégénérés.
Il est intéressant de traduire les équations du mouvement Born-Oppenheimer dans
le cas d'un hamiltonien à une particule. Dans le cas, par exemple, de l'approximation
Hartree-Fock, dénie par le minimum variationnel de la valeur moyenne de l'énergie
hΨ0 |HeHF |Ψ0 i, la fonction d'onde Ψ0 est développée sur un seul déterminant de Slater
Ψ0 = det{ψi } et les orbitales ψi sont soumises aux contraintes d'orthonormalisation
hψi |ψj i = δij . La minimisation contrainte de l'énergie par rapport aux orbitales
¯
¯
min{hΨ0 |He |Ψ0 i}¯¯
HF
(A.29)
{ψi } hψi |ψj i=δij

peut être traduite selon le formalisme de Lagrange :

X
L = −hΨ0 |HeHF |Ψ0 i + Λij (hψi |ψj i − δij ) (A.30)
i,j

où les Λij représentent les multiplicateurs de Lagrange associés aux contraintes d'or-
thonormalisation. La diérenciation de ce lagrangien (A.30) par rapport aux orbitales
δL
=0 (A.31)
δψi⋆
conduit aux équations Hartree-Fock usuelles :
X
HeHF ψi = Λij ψj (A.32)
j

La forme canonique HeHF ψi = εi ψi de ces équations peut être obtenue après une trans-
formation unitaire. Les équations du mouvement correspondantes aux relations (A.26)
et (A.27) deviennent alors :
¯
¨ ¯
~ I = −∇I min{hΨ0 |He |Ψ0 i}¯
MI R HF
¯ (A.33)
{ψi } hψi |ψj i=δij
II Dynamique Moléculaire Ab Initio 13

HeHF ψi = εi ψi (A.34)

Un jeu d'équations similaires peut être obtenu dans le cadre de la théorie de la fonction-
nelle de la densité (DFT). Dans ce cas HeHF est remplacé par l'hamiltonien Kohn-Sham
HeKS .
Les premières applications de la dynamique moléculaire Born-Oppenheimer ont été
eectuées dans le cadre d'une description semi-empirique du problème de la structure
électronique [21, 22]. La première implémentation ab initio vint plus tard [23] dans
le cadre de l'approximation Hartree-Fock. Mais ce type d'approche, très coûteuse en
temps de calcul, ne sortit de la condentialité que récemment (durant les années 1985
1990), grâce au développement de capacités informatiques et de programmes de calcul
de structure électronique plus performants. Indubitablement, l'utilisation des méthodes
de la théorie de la fonctionnelle de la densité en chimie quantique  qui s'est déve-
loppée durant la même période  a beaucoup contribué à l'essor de ces méthodes.
Enn, l'approche proposée par Car et Parrinello [7] en 1985 a certainement favorisé
l'utilisation de la dynamique moléculaire ab initio en chimie quantique.

II.3 Dynamique moléculaire  Car-Parrinello 

Le but avoué de la méthode  Car-Parrinello  est de réduire, de façon signi-
cative, le coût computationnel des approches de dynamique moléculaire. Elle peut
être considérée comme combinant les avantages des dynamiques de type Ehrenfest et
Born-Oppenheimer. Dans les dynamiques d'Ehrenfest, l'échelle de temps et ainsi le
pas de temps utilisé pour intégrer les équations de mouvement sont gouvernés par la
dynamique intrinsèque des électrons. Le mouvement des électrons étant plus rapide
que le mouvement des noyaux, le pas de temps le plus grand utilisable sera celui qui
permet d'intégrer l'équation d'évolution électronique. Au contraire, la dynamique des
électrons n'est pas considérée dans l'approche de type Born-Oppenheimer, permettant
ainsi d'intégrer les équations de mouvement sur l'échelle de temps des noyaux. Tou-
tefois, cela signie que le problème de la structure électronique doit être résolue de
14 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

manière auto-cohérente à chaque pas de la dynamique moléculaire ; ce n'est pas néces-

saire dans l'approche de type Ehrenfest puisqu'il est possible de propager la fonction
d'onde en appliquant l'hamiltonien dépendant du temps à la fonction d'onde initiale,
elle même obtenue initialement de façon auto-cohérente. L'idée de base de l'approche
Car-Parrinello consiste à exploiter la séparation entre les échelles de temps du mouve-
ment rapide des électrons et celui, plus lent, des noyaux. Cette diérence d'échelle de
temps est alors vue comme une diérence d'échelle d'énergie d'un point de vue de la
mécanique classique. Ainsi, le problème mixte classique/quantique devient un problème
purement classique avec deux échelles d'énergie bien séparées. Cependant d'un point
de vue dynamique, le système électronique ne dépend plus explicitement du temps
de la même manière. De plus, la quantité essentielle, l'énergie électronique du système
hΨ|He |Ψi, dépend de manière paramétrique des positions nucléaires {R
~ I } mais aussi de

la fonction d'onde Ψ. En mécanique classique, la force exercée sur les noyaux s'exprime
comme la dérivée du lagrangien par rapport aux positions nucléaires. Ceci suggère une
dérivée analogue par rapport aux orbitales moléculaires, qui peut être interprétée, d'un
point de vue classique, comme la force s'exerçant sur ces orbitales. Ainsi, R. Car et M.
Parrinello postulèrent le lagrangien suivant :
Z
1 X 1X ~˙ I
2
LCP = µ d~r|ψ̇i (~r)|2 + MI R
2 i 2 I
X µ Z ¶ h i (A.35)
+ Λij ∗ ~I}
d~rψi (~r)ψj (~r) − δij − Ee {ψi }, {R
i,j

où µ est une  masse ctive  associée aux orbitales moléculaires. L'avant dernier
terme du lagrangien Car-Parrinello est associé aux contraintes d'orthonormalisation des
orbitales moléculaires. Les équations de Newton correspondantes sont obtenues à partir
des équations d'Euler-Lagrange associées, appliquées aux noyaux et aux orbitales :

d ∂LCP ∂LCP
= (A.36)
dt ∂ R~˙ I ~I
∂R
d δLCP δLCP
⋆ = (A.37)
dt δ ψ̇i δψi
II Dynamique Moléculaire Ab Initio 15

Les équations du mouvement génériques Car-Parrinello ont ainsi la forme suivante :

µZ ¶
¨ ∂Ee ∂ X
~
MI R I = − + Λij ⋆
d~rψi (~r, t)ψj (~r, t) − δij (A.38)
~ I ∂R
∂R ~I
ij
X µZ ¶
δEe δ
µψ̈i (~r, t) = − ⋆ + ⋆ Λij ⋆
d~rψi (~r, t)ψj (~r, t) − δij (A.39)
δψi δψi ij

Il faut noter que pour des raisons d'homogénéité, la dimension du paramètre de masse
ctive associé aux orbitales moléculaires est une énergie multipliée par un temps au
carré.
Ainsi, les noyaux évoluent dans le temps à une certaine température  instanta-
P ~˙ 2 alors qu'une température ctive proportionnelle à
née  proportionnelle à I MI R I
P
µ i hψ̇i |ψ̇i i est associée aux degrés de liberté électroniques. Dans cette terminologie,
des  électrons froids  signie que le sous-système électronique est proche de la sur-
face  exacte  Born-Oppenheimer. La fonction d'onde de l'état fondamental optimisée
pour la conguration nucléaire initiale restera alors proche de la solution exacte du-
rant son évolution temporelle si elle  reste  à une température susamment basse.
Finalement, l'idée de base de l'approche Car-Parrinello consiste à traiter les orbitales
moléculaires comme des variables dynamiques classiques.
L'approche originale de Car et Parrinello fut initialement développée et implémentée
dans le cadre de la théorie de la fonctionnelle de la densité en utilisant des ondes
planes comme fonctions de base pour la structure électronique [8]. L'utilisation de ce
lagrangien étendu avec une description de la structure électronique développée sur un
jeu de fonctions gaussiennes a initialement été proposée par M.J. Field en 1991 [9,
24]. Cette idée fut rapidement reprise par E.A. Carter [25, 26, 27, 28, 29, 30, 31,
32]. Malheureusement, ces travaux conclurent quant à l'inecacité de la dynamique
moléculaire Car-Parrinello en base locale [29]. Récemment, une alternative envisageant
la propagation classique  non pas des orbitales moléculaires  mais de la matrice
densité, a été proposée [33]. Les détails des diverses implémentations de dynamique
moléculaire Car-Parrinello seront donnés dans la partie C.
16 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

III Intégration des Équations du Mouvement

En dynamique moléculaire classique, la trajectoire des diérentes composantes du

système moléculaire est générée par intégration des équations du mouvement, qui pour
chaque particule I , s'écrivent :
~I
d2 R ~
∂Ep ({R})
MI = − (A.40)
dt2 ∂R~ I (t)

~ est la fonction énergie potentielle du système à N particules, qui ne dépend que

Ep ({R})
~ I }. Les équations (A.40) sont intégrées numériquement
des coordonnées cartésiennes {R
en utilisant un pas de temps, δt, innitésimal, garantissant la conservation de l'énergie
du système. Le même type de schéma d'intégration peut être appliqué à la propagation
des orbitales moléculaires dans une approche de type Car-Parrinello.
Il est cependant parfaitement illusoire d'espérer générer une trajectoire exacte aux
temps longs, dès lors que les équations newtoniennes du mouvement sont résolues
numériquement, avec un pas d'intégration ni. L'exactitude de la solution des équa-
tions (A.40) n'est, cependant, pas aussi cruciale qu'il n'y paraît au premier abord. Ce
qui importe, en revanche, est que le comportement statistique de la trajectoire soit,
quant à lui, correct, pour garantir la reproduction des propriétés dynamiques et ther-
modynamiques du système avec susamment de précision. Cette condition peut être
remplie si l'intégrateur utilisé pour propager le système possède la propriété de sym-

plecticité. Un propagateur dit symplectique conserve la métrique invariante de l'espace

des phases (cette notion sera précisée par la suite). Il en résulte que l'erreur associée à
ce propagateur est nécessairement bornée :
n ¯ ¯
1 X ¯¯ E (kδt) − E (0) ¯¯
lim ¯ ¯≤² (A.41)
n→∞ n E (0)
k=1

où n désigne le nombre de pas de la simulation, E (t) l'énergie totale du système au

temps (t) et ² la borne supérieure de conservation de l'énergie. En supposant que le pas
d'intégration soit limité, l'intégration numérique des équations du mouvement ne donne
pas lieu à une croissance erratique de l'erreur de conservation de l'énergie, qui pourrait
III Intégration des Équations du Mouvement 17

aecter signicativement le comportement statistique de la dynamique moléculaire aux

temps longs. An de préciser l'équivalence statistique entre la solution numérique et la
solution exacte des équations du mouvement, il est utile de rappeler quelques points
de l'approche de Hamilton et le lien avec la mécanique statistique.

III.1 Le point de vue de Hamilton et la mécanique statistique

L'hamiltonien d'un système constitué de N particules interagissant entre elles est

déni selon :
X P~ 2
~ I }, {P~I }) =
H ({R I ~ I })
+ V ({R (A.42)
I
2M I

où {P~I } sont les impulsions des particules dénies tel que P~I = MI V~I et V ({R
~ I })

désigne le potentiel interparticulaire. Les forces {F~I } agissant sur ces particules sont
conservatives puisqu'elles dérivent de ce potentiel :

∂V
F~I = − (A.43)
~I
∂R

Les équations du mouvement (A.40) peuvent être obtenues à partir de la relation (A.42)
et des relations de Hamilton :
~
~˙ I = ∂H = PI
R (A.44)
∂ P~I MI
˙ ∂H ∂V
P~I = − =− = F~I ({R
~ I }) (A.45)
~I
∂R ∂R~I

En dérivant par rapport au temps la première relation de Hamilton puis en substi-

tuant ce résultat dans la seconde, les équations du mouvement (A.40) sont facilement
retrouvées. Ainsi, l'état classique d'un système à n'importe quel instant (t) est totale-
~ I } et toutes les impulsions {P~I }
ment déterminé en spéciant toutes les positions {R
associées. Autrement dit, on peut dénir un vecteur Γ constitué du jeu complet des
positions et des impulsions Γ = ({P~I }, {R
~ I }). Ce vecteur est un élément de l'espace

des phases. L'espace des phases à 2N dimensions est donc la réunion de tous les états
classiques accessibles au système.
18 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

Deux importantes propriétés des équations du mouvement doivent être précisées.

La première est la réversibilité temporelle, i.e. ces équations du mouvement gardent la
même forme si on eectue la transformation t → −t. La conséquence de cette symétrie
du temps est que le comportement physique microscopique est indépendant du sens
de l'écoulement du temps. La seconde propriété importante de ces équations est que
l'hamiltonien (A.42) est un invariant du système. Ceci se vérie facilement en dérivant
par rapport au temps cet hamiltonien H et en utilisant les relations de Hamilton (A.44)
et (A.45) pour les positions et les impulsions :
dH X · ∂H ˙ ∂H ˙
¸ X·
∂H ∂H ∂H ∂H
¸
= ~ +
R P~ = − =0 (A.46)
dt ∂ ~I I
R ∂ ~I I
P ∂ ~ I ∂ P~I
R ∂ ~I ∂ R
P ~I
I I

L'invariance de l'hamiltonien est équivalente à la conservation de l'énergie totale du

système et fournit un lien fondamental avec la mécanique statistique. Rappelons que
la mécanique statistique relie les détails microscopiques d'un système aux observables
physiques telles que les propriétés thermodynamiques à l'équilibre, les coecients de
diusion, etc. La mécanique statistique est basée sur le concept d'ensemble de Gibbs :
un grand nombre de congurations individuelles microscopiques d'un très grand sys-
tème conduit à des propriétés macroscopiques ; ce qui implique notamment qu'il n'est
pas nécessaire de connaître précisément le mouvement de chaque particule pour prédire
les propriétés macroscopiques de ce système. Il est susant de moyenner sur un grand
nombre de systèmes identiques, chacun étant dans un état microscopique diérent ;
i.e. les observables macroscopiques d'un système sont formulées en terme de moyenne
d'ensemble . Les diérents ensembles statistiques sont caractérisés par des variables
thermodynamiques qui peuvent être xées, telles que l'énergie totale E , la tempéra-
ture T , la pression P , le potentiel chimique µ, etc. Un de ces ensembles est l'ensemble
microcanonique , noté N V E , caractérisé par un nombre N constant de particules, un
volume V constant et l'énergie totale E constante. On peut aussi citer l'ensemble ca-

nonique (N V T ), l'ensemble isotherme-isobare (N P T ) et l'ensemble grand canonique

(µV T ). Ces variables thermodynamiques doivent être vues comme des paramètres ex-
périmentaux de contrôle qui spécient les conditions sous lesquelles l'expérience est
III Intégration des Équations du Mouvement 19

conduite.
Considérons maintenant un système de N particules occupant un volume V et
évoluant selon les équations du mouvement de Hamilton. Selon la relation (A.46), l'ha-
miltonien H est un invariant du système et l'énergie totale E doit être conservée. De
plus, le nombre de particules et le volume sont supposés xes. Ainsi, une trajectoire
dynamique de ce système génère une série d'états classiques ayant pour constantes N ,
V et E correspondant à l'ensemble microcanonique. Si cette dynamique génère tous les
états accessibles au système avec N , V et E xés, alors une moyenne sur cette trajec-
toire conduira au même résultat qu'une moyenne dans l'ensemble microcanonique. La
condition de conservation de l'énergie, qui impose une  restriction  sur l'ensemble des
états accessibles au système, dénit une hypersurface dans l'espace des phases appelée
surface d'énergie constante ou d'isoénergie. Un système évoluant en accord avec les
relations de Hamilton restera sur cette surface.
Le fait de supposer qu'un système, au bout d'un temps inni, parcourt la totalité
de cette hypersurface d'énergie constante porte le nom d'hypothèse ergodique. Ainsi,
à l'aide de l'hypothèse ergodique, la moyenne sur une trajectoire d'un système sujet
aux relations de Hamilton conduit à la moyenne de l'ensemble microcanonique. D'un
point de vue mathématique, si A({R
~ I }, {P~I }) est une fonction correspondant à une

observable physique, la moyenne dans l'ensemble microcanonique de cette observable

est :
Z Y h i
1
hAiN V E = ~ I dP~I A({R
dR ~ I }, {P~I }) δ H ({R
~ I }, {P~I }) − E (A.47)
h3N Ω I

où h représente la constante de Planck et Ω désigne la fonction de partition microca-

nonique dénie selon :
Z Y h i
1
Ω = 3N ~ I dP~I δ H ({R
dR ~ I }, {P~I }) − E (A.48)
h I

Cette fonction de partition permet de  compter  le nombre d'états accessibles au

système sujet aux relations de Hamilton. En pratique, les trajectoires en dynamique
moléculaire simulent l'évolution sur un temps ni τ du système et la moyenne de
20 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

l'observable A sur cette trajectoire est donnée par la relation :

Z τ
eτ = 1
A ~ I (t)}, {P~I (t)})
dt A({R (A.49)
τ 0

L'hypothèse ergodique est donc équivalente à l'égalité suivante :

eτ = hAi
lim A (A.50)
τ →∞

La signication de l'équivalence statistique entre les trajectoires numériques obtenues

par des simulations de dynamique moléculaire et les trajectoires exactes est maintenant
claire. Même si la trajectoire numérique s'éloigne de plus en plus de la trajectoire
exacte au cours de la simulation, tant que cette trajectoire conserve l'énergie totale
du système avec une précision numérique satisfaisante, cette trajectoire numérique
génère un ensemble de congurations appartenant à l'hypersurface d'énergie constante.
L'existence d'une borne supérieure de l'erreur numérique, i.e. le caractère symplectique
des trajectoires, est précisée par la suite. Ainsi, à l'aide de l'hypothèse ergodique, la
moyenne d'ensemble d'une observable peut être obtenue à partir d'une ou plusieurs
trajectoires numériques.

En dépit de l'utilité de la dynamique moléculaire hamiltonienne, sa principale res-

triction est évidente : l'ensemble généré est l'ensemble microcanonique. En eet, la
plupart des mesures expérimentales sont eectuées dans des conditions de température
et de pression constantes ou des conditions de température et de volume constants.
Ainsi, an de décrire les propriétés thermodynamiques d'un système libre de taille -
nie sous de telles conditions, il est nécessaire de générer l'ensemble thermodynamique
correspondant. Une des méthodes les plus intéressantes pour générer ces ensembles est
basée sur la dynamique non hamiltonienne. Cette approche est précisée à la partie D.
Auparavant, la propriété de symplecticité des algorithmes d'intégration et l'intégration
numérique des équations du mouvement sont abordées.
III Intégration des Équations du Mouvement 21

III.2 La propriété de symplecticité

L'évolution temporelle d'un système, dans l'ensemble microcanonique, est régie par
les équations de Hamilton qui peuvent se mettre sous la forme condensée suivante :

t
Γ̇ = gH (Γ, t) (A.51)

où Γ = ({R
~ I }, {P~I }) désigne un vecteur de l'espace des phases et g t représente le
H

 champ des vitesses  au point Γ, déni par :

 
∂H
~
t
gH = ∂P  (A.52)
− ∂H
~
∂R

On appelle aussi gH
t
le ot hamiltonien par analogie avec la mécanique des uides et,
par une même analogie avec un écoulement incompressible, on peut démontrer qu'un
volume v quelconque de l'espace des phases est conservé par le ot hamiltonien gH
t
,
soit :
dv
=0 (A.53)
dt
Soit un point quelconque de v , lors d'une évolution temporelle, ce point va suivre une
trajectoire précise dans l'espace des phases, imposée par les relations de Hamilton. Il
va en être de même pour tous les autres points appartenant à v . Ce volume v va donc
être déformé au cours de l'évolution, de manière analogue à un élément de uide pris
dans un écoulement. Ce théorème d'incompressibilité du ot, ou encore de conservation
du volume de l'espace des phases, stipule donc simplement que la déformation due au
ot hamiltonien conserve le volume.
Ce théorème permet, de plus, de retrouver le théorème de Liouville. Dénissons la
densité d'états ρ = dn
dv
d'un système au voisinage d'un point de l'espace des phases ;
dn représente le nombre d'états contenus dans l'élément de volume dv de l'espace
des phases. Le théorème d'incompressibilité du ot permet de préciser que dv reste
constant lors de l'évolution du système. Par ailleurs, le nombre d'états dn reste égale-
ment constant, car aucune trajectoire ne peut traverser la surface frontière dénissant
22 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

le volume dv . En eet, si c'était le cas, alors deux conditions initiales diérentes mè-
neraient au même état (celui situé à la frontière), ce qui contredit le principe d'unicité
des solutions des équations de Hamilton. Finalement, le quotient dn
dv
est une constante,
et la relation (A.53) peut être écrite sous la forme d'une relation de conservation pour
la densité d'état, connue sous le nom d'équation de Liouville :
dρ ∂ρ
= + Γ̇∇ρ (Γ, t) = 0 (A.54)
dt ∂t
L'équation de Liouville et le théorème associé sont donc intimement reliés à l'inva-
riance de la métrique de l'espace des phases. L'espace des phases étant un espace à 2N
dimensions, le produit de diérentielles :
Y
dv = ~ I dP~I
dR (A.55)
I

peut être considéré comme un élément innitésimal de volume. Or, d'après le théorème
d'incompressibilité, tout élément de volume de l'espace des phases pour un système
hamiltonien est un invariant. Par conséquent, la métrique dv de l'espace des phases est
invariante.
La propriété de symplecticité des propagateurs, brièvement introduite auparavant,
apparaît alors fondamentale : un propagateur est dit symplectique s'il conserve le vo-
lume de l'espace des phases, ou bien encore sa métrique ; le théorème de Liouville est
alors vérié et la structure formelle des relations de Hamilton est préservée. L'un des
plus simples propagateurs symplectiques que l'on puisse utiliser est l'algorithme de
Verlet, décrit par la suite.

III.3 Algorithme de Verlet

Il existe plusieurs approches pour intégrer numériquement les équations de Newton

du mouvement. Sans doute la plus simple, l'algorithme de Verlet [34] envisage un
développement en série de Taylor de la position en (t − δt) et en (t + δt), conduisant
à:
R ~ I (t) − R
~ I (t + δt) = 2R ~ I (t) + O((δt)3 )
~ I (t − δt) + (δt)2 A (A.56)
III Intégration des Équations du Mouvement 23

où A
~ I (t) désigne l'accélération au temps (t) de la particule I . L'algorithme de Verlet

est peu coûteux numériquement, tant du point de vue de la mémoire nécessaire que
du point de vue de la simplicité de l'implémentation. De plus il conserve très bien
l'énergie et la quantité de mouvement totale, même pour des pas de temps relative-
ment grands [35]. L'expression (A.56) ci-dessus est symétrique, l'algorithme de Verlet
est donc naturellement réversible dans le temps. Ces qualités en font la méthode de
résolution des équations classiques du mouvement la plus employée dans les simulations
de dynamique moléculaire. Il est à noter que les vitesses n'interviennent pas directe-
ment dans l'algorithme. Bien que non nécessaires à la description des trajectoires, leur
calcul est obligatoire pour évaluer l'énergie cinétique, T (P~ ), qui ne dépend que des
impulsions P~I , et, par conséquent, l'énergie totale E . Ces dernières peuvent néanmoins
être estimées numériquement à l'instant (t) par :
~ I (t + δt) − R
R ~ I (t − δt) ¡ ¢
V~I (t) = + O (δt)2 (A.57)
2δt

La délicate manipulation des vitesses a été à l'origine de plusieurs variantes de cet algo-
rithme. Le premier avatar de la méthode de Verlet est l'algorithme  saute-mouton 
(leap-frog en anglais), qui nécessite de stocker les positions et les accélérations aux
temps (t), et les vitesses aux temps (t − δt/2) :

R ~ I (t) + V~I (t + δt )δt

~ I (t + δt) = R (A.58)
2
δt δt
V~I (t + ) = V~I (t − ) + A ~ I (t)δt (A.59)
2 2

Dans la pratique, la première étape est le calcul de V~I (t + δt2 ) à partir duquel on déduit
V~I (t), nécessaire à l'évaluation du terme cinétique T (P~ ), selon :

V~I (t + δt
) + V~I (t − δt
)
V~I (t) = 2 2
(A.60)
2

L'élimination des vitesses entre les deux équations (A.58) et (A.59) montre l'équivalence
entre l'algorithme de Verlet avec l'algorithme leap-frog. Enn, l'algorithme de Verlet aux
vitesses corrige le défaut principal des précédents, à savoir la dénition des vitesses, dont
24 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

l'erreur associée est en O ((δt)2 ). L'incorporation explicite des vitesses dans l'algorithme
de Verlet peut s'écrire :
2
~ I (t + δt) = R
R ~ I (t) + V~I (t)δt + A~ I (t) (δt) (A.61)
2
~ ~
AI (t) + AI (t + δt)
V~I (t + δt) = V~I (t) + δt (A.62)
2

Les trois algorithmes présentés ci-dessus sont rigoureusement équivalents quant aux
trajectoires produites dans l'espace des congurations. Les diérences se situent dans
le traitement réservé aux vitesses, donc à l'énergie cinétique T (P~ ).

III.4 Les autres propagateurs

Il existe par ailleurs une multitude de propagateurs, plus ou moins performants, à

pas de temps variable ou xe [36, 37, 38]. D'une manière générale, les propagateurs  à
pas de temps xe  plus précis que l'algorithme de Verlet, nécessitent naturellement la
connaissance des dérivées secondes de l'énergie par rapport aux positions nucléaires ou
par rapport au temps. Plusieurs algorithmes sophistiqués ont été proposés pour proter
au mieux de la connaissance du hessien [39]. Malheureusement l'évaluation de ce terme,
particulièrement gourmande en temps de calcul, limite la taille du système étudié ainsi
que le temps total de simulation. Ce type de propagateur est le plus souvent utilisé
pour la réalisation de trajectoires dites de référence. Leur utilisation  routinière 
est encore proscrite du fait de leur coût calculatoire. Un propagateur d'ordre quatre
communément utilisé, ne nécessitant pas l'évaluation du hessien, est l'algorithme de
prédiction-correction de Gear [36, 35]. Cet algorithme utilise une étape de prédiction :
2 3
R ~ I (t) + δtV~I (t) + (δt) F~I (t) + (δt) J~I (t)
~ I∗ (t + δt) = R (A.63)
2MI 6MI
2
δt ~ (δt) ~
V~I∗ (t + δt) = V~I (t) + FI (t) + JI (t) (A.64)
MI 2MI
F~I∗ (t + δt) = F~I (t) + δtJ~I (t) (A.65)

suivie par une étape de correction :

III Intégration des Équations du Mouvement 25

2 h i
~ I (t + δt) = R
R ~ I∗ (t + δt) + (δt) F~I (t + δt) − F~I∗ (t + δt) (A.66)
12MI
~ ~ ∗ 5δt h ~ ~ ∗
i
VI (t + δt) = VI (t + δt) + FI (t + δt) − FI (t + δt) (A.67)
12MI
1 h~ i
J~I (t + δt) = J~I (t) + FI (t + δt) − F~I∗ (t + δt) (A.68)
δt
où F~I (t) représente les forces à l'instant (t) exercées sur la particule I, i.e. F~I (t) =
~ I Ee (t).
−∇ La variable J~I (t) est une approximation de la dérivée première des forces

par rapport au temps :

dF~I (t)
J~I (t) ≈ (A.69)
dt
Ce propagateur, bien que d'ordre quatre, n'est ni réversible dans le temps ni symplec-

tique. G. J. Martyna et M. E. Tuckerman ont récemment proposé [40] un intégrateur

d'ordre quatre, réversible et symplectique, ne nécessitant pas l'évaluation du hessien.

Pour ce faire, une approximation des dérivées temporelles des forces doit être formelle-

ment introduite. Ceci est accompli en étendant l'espace des phases à un jeu de fonctions

{F~ ∗ , ~vF ∗ , J~∗ , ~vJ ∗ }. La nouvelle variable J~∗ représente, comme dans la relation (A.69),

une approximation de la dérivée première des forces par rapport au temps et la variable

~vJ ∗ représente, quant à elle, une approximation de la dérivée seconde des forces par

rapport au temps. Cet algorithme peut se résumer à l'aide des relations suivantes :

2 h i
R~ I (t + δt) = R~ I (t) + (δt)V~I (t) + (δt) 5F~I (t) − F~ ∗ (t)
I
8MI
(A.70)
(δt)3 (δt)4
+ ~vFI∗ (t) + ~vJ ∗ (t)
6MI 48MI
δt h ~ i
V~I (t + δt) = V~I (t) + 4FI (t) + F~I∗ (t) + 3F~I (t + δt)
8MI
(A.71)
2
(δt)
+ ~vF ∗ (t)
8MI I

~ ∗ 1h ~ ~ ∗ ~
i δt
FI (t + δt) = 2FI (t) − FI (t) + FI (t + δt) + ~vFI∗ (t) (A.72)
2 2
1 3 h~ ~ ∗
i
~vFI (t + δt) = − ~vFI (t) +
∗ ∗ FI (t + δt) − FI (t) (A.73)
2 2(δt)
3 3 h~ ~ ∗
i
~vJ (t + δt) = −~vJ (t) − ~vFI (t) +
∗ ∗ ∗ FI (t + δt) − FI (t) (A.74)
δt (δt)2
26 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

Il faut noter la présence d'une faute de frappe dans l'article [40] relative à l'équa-
tion (A.74) et qui concerne le coecient multiplicatif du troisième terme. Cet inté-
grateur, dont l'utilisation n'est que trop condentielle, possède des performances bien
supérieures par rapport au schéma de Verlet [41]. En eet la dérive énergétique observée
pour un pas de propagation donné est bien plus importante si le propagateur est de type
Verlet que dans le cas du prédicteur-correcteur symplectique d'ordre quatre précédent.
L'utilisation de ce dernier propagateur est certainement limitée par la diculté de le
généraliser à d'autres outils de la dynamique moléculaire tels que les contraintes géo-
métriques ou bien les thermostats. Néanmoins son implémentation dans le cadre d'une
approche Born-Oppenheimer dans l'ensemble microcanonique demeure facilement ac-
cessible, autorisant ainsi des simulations avec des pas de temps d'intégration supérieurs
à ceux employés avec l'algorithme de Verlet pour une qualité de simulation identique.
Ce propagateur possède de plus les qualités remarquables de réversibilité temporelle
et de symplecticité, puisque par construction, il est basé sur l'opérateur d'évolution de
Liouville. Cette notion ainsi que le lien avec la propriété de symplecticité sont expliqués
par la suite.

III.5 L'opérateur de Liouville et l'intégration numérique

Auparavant, il a été montré que, même si la trajectoire obtenue numériquement

en intégrant les équations du mouvement diverge de la trajectoire exacte, elle reste
statistiquement équivalente à la trajectoire exacte. Cette équivalence n'est vériée que si
l'énergie totale du système est conservée avec une certaine tolérance. L'erreur commise
pour les trajectoires numériques est d'autant plus petite que le pas de temps utilisé
pour intégrer les équations du mouvement est petit. Cependant, plus le pas de temps
est petit, plus le nombre d'itérations pour simuler un temps τ déni est grand. Ainsi,
l'intégration numérique des équations du mouvement est un compromis entre le plus
grand pas de temps possible et une tolérance acceptable quant à la conservation de
l'énergie.
III Intégration des Équations du Mouvement 27

L'approche présentée ici est basée sur la formulation d'un opérateur d'évolution de
la mécanique classique [42, 43, 44]. Considérons un système régi par la mécanique de
Hamilton. Les équations du mouvement peuvent prendre la forme réduite suivante :

Γ̇ = iLΓ (A.75)

où iL représente l'opérateur de Liouville, déni selon :

iL = {. . . , H }
X · ∂H ∂ ∂H ∂
¸
≡ −
I ∂ ~I ∂ R
P ~I ∂R~ I ∂ P~I (A.76)
" #
X P~I ∂ ∂
= + F~I
MI ∂ R~I ∂ P~I
I

où la notation {f, H } désigne le crochet de Poisson entre l'hamiltonien H et une

fonction f . Il faut remarquer dès à présent que tout propagateur construit à partir de
l'opérateur d'évolution de Liouville est symplectique. En eet, la métrique de l'espace
des phases est conservée puisque l'hamiltonien H est un invariant et le crochet de
Poisson de deux invariants est lui-même un invariant du système.
La relation (A.75) admet comme solution formelle :

Γ(t) = eiLt Γ(0) (A.77)

Cette relation est le point de départ pour la dérivation des procédures d'intégration
numérique. L'opérateur unitaire U (t) = eiLt est le propagateur classique. Son action
sur Γ(0) ne peut être déterminée analytiquement que pour quelques cas simples. Toute-
fois cette relation formelle peut être utilisée pour générer des propagateurs numériques
à l'aide d'approximations simpliant la relation (A.77). Supposons, par exemple, que
l'opérateur de Liouville iL puisse être décomposé en deux parties, iL = iL1 + iL2 , de
telle façon que l'action du propagateur classique sur Γ(0) pour chaque partie puisse
être analytiquement déterminée. Ce propagateur classique peut être ré-écrit à l'aide du
28 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

théorème de Trotter [45], qui stipule :

U (t) = eiLt = exp [(iL1 + iL2 )t]

· µ ¶ µ ¶ µ ¶¸M (A.78)
iL2 t iL1 t iL2 t
= lim exp exp exp
M →∞ 2M M 2M
Posons ensuite t/M = δt. Pour M ni, l'approximation suivante peut être faite :
M
Y µ ¶ µ ¶
iL1 δt iL1 δt ¡ ¢
exp(iLM δt) ≈ exp × exp (iL2 δt) × exp + O (M δt)3 (A.79)
k=1
2 2

soit, avec M = 1 :
µ ¶ µ ¶
iL1 δt iL1 δt ¡ ¢
exp(iLδt) ≈ exp × exp (iL2 δt) × exp + O (δt)3 (A.80)
2 2
On peut donc dénir un propagateur discret en temps :
µ ¶ µ ¶
iL1 δt iL1 δt
G (δt) = exp × exp (iL2 δt) × exp
2 2
µ ¶ µ ¶ (A.81)
δt δt
= U1 × U2 (δt) × U1
2 2
Les trois composantes de ce propagateur G (δt) sont chacune unitaires, le propagateur
G (δt) est donc aussi unitaire, i.e. G −1 (δt) = G † (δt) = G (−δt). Ceci signie que tout
propagateur basé sur la factorisation de Trotter est réversible dans le temps.
Pour mieux préciser l'action de chaque composante du propagateur, dénissons :

Γ1 [δt; Γ(0)] = U1 (δt) Γ(0) (A.82)

Γ2 [δt; Γ(0)] = U2 (δt) Γ(0) (A.83)

correspondant respectivement à l'état au temps (δt) alors que le système, initialement

dans l'état Γ(0) est soumis soit à U1 soit à U2 . En appliquant ainsi l'opérateur G (δt),
déni par la relation (A.81), à un état initial Γ(0), on obtient :
¶ µ µ ¶
δt δt
Γ(δt) = U1 × U2 (δt) × U1 Γ(0)
2 2
µ ¶ · ¸
δt δt
= U1 × U2 (δt) Γ1 ; Γ(0)
2 2
µ ¶ ½ · ¸¾ (A.84)
δt δt
= U1 Γ2 δt; Γ1 ; Γ(0)
2 2
µ ½ · ¸¾¶
δt δt
= Γ1 ; Γ2 δt; Γ1 ; Γ(0)
2 2
III Intégration des Équations du Mouvement 29

Considérons, par exemple, le découpage de l'opérateur de Liouville suivant :

X ∂ X P~I ∂
iL1 = F~I et iL2 = (A.85)
I ∂ P~I I
~I
MI ∂ R

En appliquant la relation (A.81), ceci conduit au propagateur :

à ! à ! à !
δt X ∂ X ∂ δt X ∂
G (δt) = exp F~I exp δt P~I exp F~I (A.86)
2 I ∂ P~I I
~I
∂R 2 I ∂ P~I

De plus, en développant en série de Taylor l'exponentielle, on peut montrer que tout

opérateur de la forme ea∂/∂x agissant sur une fonction f (x), avec a indépendant de x,
s'écrit :
µ ¶ "
∞ ¡ ∂ ¢k
# "∞ #
∂ X a ∂x X ak ∂ k
exp a f (x) = f (x) = k
f (x)
∂x k! k! ∂x
k=0
∞
k=0
(A.87)
X ak ∂ k f (x)
= = f (x + a)
k=0
k! ∂xk

Ainsi, en appliquant successivement les trois composantes du propagateur, déni par

³ ´
~ I (0)}; {P~I (0)} , on obtient :
la relation (A.86), à un état Γ(0) = {R
µ ¶
δt δt
P~I = P~I (0) + F~I (0)
2 2
µ ¶
δt ~ δt
~ ~
RI (δt) = RI (0) + PI (A.88)
MI 2
µ ¶
δt δt
P~I (δt) = P~I + F~I (δt)
2 2
Ces relations peuvent être exprimées sous la forme condensée suivante :
2
~ I (t) + δt P~I (t) + (δt) F~I (t)
~ I (t + δt) = R
R (A.89)
MI 2MI
³
δt ~ ´
~ ~
PI (t + δt) = PI (t) + ~
FI (t) + FI (t + δt) (A.90)
2

On retrouve alors l'algorithme Verlet aux vitesses, présenté précédemment. Cette ma-
nière élégante de retrouver le propagateur de Verlet permet de justier ses propriétés
de symplecticité et de réversibilité temporelle. Un développement de la factorisation de
Trotter à des ordres plus élevés conduit à des propagateurs plus précis. Toutefois, ce
développement fait apparaître les dérivées des forces par rapport aux positions, termes
30 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

trop gourmands en temps de calcul. Le développement d'intégrateurs au moyen de la

factorisation de Trotter s'avère très utile pour des systèmes dont les échelles de temps
caractéristiques sont multiples. C'est en particulier le cas lorsque le système est couplé
à un  réservoir , tel qu'un thermostat. Cet aspect est plus détaillé dans la partie C
ainsi qu'en annexe.
Les diérents algorithmes propagateurs permettent d'estimer l'évolution temporelle
d'un système moléculaire avec plus ou moins d'ecacité et de précision, ceci à partir
d'un jeu de conditions initiales données. Le choix de ce point de l'espace des phases de
départ peut grandement inuencer le résultat des simulations de dynamique molécu-
laire.

IV Choix des Conditions Initiales

À énergie totale E xée, la conguration de départ {R(0)}
~ de la simulation peut
en principe être quelconque dans la mesure où elle appartient à l'espace des phases
autorisé, ce qui ici s'écrit simplement V ({R
~ I }) < E . Ceci suppose que le système a

un comportement ergodique. Il est d'usage de débuter une simulation avec une géomé-
trie pertinente, en général un isomère de basse énergie, si possible le plus bas connu.
Cette règle n'est pas absolue, notamment pour les simulations à haute énergie ou dans
l'ensemble canonique, mais la règle énoncée plus haut impose de toute façon que la
conguration de départ soit de plus en plus proche du minimum global à mesure que
l'énergie totale décroît.
Le problème des vitesses initiales se pose en termes diérents selon que l'on simule
un système ni ou ayant des conditions limites périodiques. Dans ce dernier cas, il sut
de vérier que la quantité de mouvement totale est nulle à l'instant initial. En général
une distribution des vitesses aléatoires ou suivant une loi de Maxwell est choisie, la
contribution du centre de masse est retranchée et l'ensemble des vitesses est multiplié
par un facteur constant pour obtenir l'énergie totale recherchée.
IV Choix des Conditions Initiales 31

Dans le cas d'un système ni, la contrainte sur le moment cinétique est plus stricte,
mais le principe reste similaire. En plus de l'énergie totale E , de la quantité de mou-
vement totale P~ , le moment cinétique L
~ doit également être imposé. Pour simplier,

les valeurs de P~ et L
~ choisies sont en général toutes deux nulles. On commence à nou-

veau par choisir des vitesses aléatoires en direction et en norme pour chaque particule,
éventuellement suivant une loi de Maxwell. La quantité de mouvement totale est alors
supprimée par un simple décalage du vecteur requis. Le nouvel ensemble de vitesses
ne satisfait pas à la contrainte sur le moment cinétique, ce dernier est annulé de la
manière suivante.

 Le moment cinétique créé par les vitesses initiales (après décalage) est :

X
~ (0) =
L ~ I × V~ (0) = I · ~ω (0)
MI R (A.91)
I
I

 Ce moment peut être compensé en enlevant la contribution de rotation globale :

(1) (0)
V~I = V~I − ~ω (0) × R
~I (A.92)

où I représente la matrice d'inertie. Finalement, les vitesses sont multipliées par un

facteur constant α pour atteindre l'énergie cinétique initiale souhaitée :
" #1/2
(1) E − V ({R~ I })
V~I (t = 0) = αV~I α= P (A.93)
1 ~ (1) 2
2 I MI (VI )

~ non nul, la vitesse de chaque

Si le système étudié possède un moment cinétique L
particule doit comporter une composante de rotation globale :

V~Irot = ~ω rot × R
~I avec ~ = I · ~ω rot
L (A.94)

Dans le calcul de l'énergie totale, il est alors nécessaire de tenir compte de l'énergie
~ I (t = 0)} doit
de rotation induite par le moment cinétique. La géométrie de départ {R
donc vérier à présent :

~ I (0)}) + 1 L
V ({R ~ † · I−1 ({R
~ I (0)}) · L
~ <E (A.95)
2
32 DYNAMIQUE MOLÉCULAIRE : GÉNÉRALITÉS

Notons cependant qu'à moment cinétique nul, une autre façon d'engendrer des condi-
tions initiales d'énergie totale E sans moment cinétique consiste à déplacer aléatoi-
rement la géométrie jusqu'à ce que l'énergie potentielle atteigne E . Il sut alors de
relâcher le système sans vitesse initiale.

V Conclusion
Il existe plusieurs façons de simuler l'évolution temporelle d'un système molécu-
laire. À partir d'une approche quantique, la dynamique moléculaire ab initio s'inscrit
comme une vision mixte où les noyaux sont des particules classiques alors que la struc-
ture électronique est décrite dans une approche quantique. Cette approche se justie
donc lorsque les eets quantiques nucléaires et les couplages adiabatiques peuvent être
négligés. Les équations du mouvement des noyaux résultantes peuvent être intégrées
dans le temps au moyen d'algorithmes propagateurs symplectiques, tel que le schéma
de Verlet, sans pour autant nécessiter l'évaluation des dérivées secondes du potentiel
électronique. La réalisation de trajectoires à l'aide de ces outils permet de générer un
jeu de points de l'espace des phases, qui peut s'identier à l'ensemble microcanonique
au regard de l'hypothèse ergodique. Ceci constitue la base de la dynamique molécu-
laire ab initio. Néanmoins un point n'a pas été abordé au cours de cette partie : la
description de la structure électronique. Ces méthodes, qui sont au c÷ur de la chimie
quantique, sont présentées dans la partie suivante.
B La Structure Électronique des
Édices Moléculaires
Dans la partie précédente, l'attention a été portée sur les diérentes possibilités d'es-
timer l'évolution temporelle de la structure nucléaire et éventuellement celle des élec-
trons. Le but de cette partie est de souligner les grandes lignes des diérentes ap-
proches de détermination de la structure électronique ainsi que le calcul des dérivées
premières de l'énergie par rapport aux positions des noyaux, les gradients nucléaires.

P our la dynamique moléculaire Born-Oppenheimer, nous avons vu que la struc-

ture électronique est déterminée à chaque pas de la propagation nucléaire comme
solution stationnaire. Cela revient à résoudre l'équation de Schrödinger indépendante
du temps :
He Ψ = Ee Ψ (B.1)

où Ψ désigne la fonction d'onde électronique et He l'hamiltonien électronique :

X ~2 X e2 X e2 ZI
He = − ∇2i + − (B.2)
2me |~ri − ~rj | ~
i i<j I,j |RI − ~
rj |

D'un autre côté, pour les méthodes de dynamique moléculaire où la propagation tem-
porelle de la structure électronique est explicite, la fonction d'onde initiale doit être
déterminée. Là encore, la structure électronique est souvent obtenue en résolvant l'équa-
tion de Schrödinger indépendante du temps (B.1). Il faut noter que les diérentes
approches qui sont exposées par la suite tiennent compte de l'approximation de Born-
Oppenheimer, énoncée précédemment (cf. partie A, Ÿ II.2 p. 11).
34 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

I Approximation Monoélectronique  Déterminant de

Slater
La principale diculté dans la résolution de l'équation (B.1) provient du terme de
répulsion biélectronique, second terme de la relation (B.2). Cependant, si on considère
que les électrons évoluent indépendamment les uns des autres, l'hamiltonien électro-
nique He peut alors s'écrire comme la somme d'hamiltoniens monoélectroniques :
n
X
He = hi (~ri ) (B.3)
i=1

L'hamiltonien monoélectronique hi rend compte du terme d'énergie cinétique de l'élec-

tron i, du potentiel attractif nucléaire qu'il ressent mais aussi de la répulsion moyenne
entre l'électron i et tous les autres. Cette approximation est dite approximation mono-
électronique. La fonction d'onde Ψ peut alors s'écrire comme un produit de fonctions
d'ondes monoélectroniques φi , elles-mêmes fonctions propres des opérateurs hi . Une
spin-orbitale φi est le produit d'une orbitale spatiale ψi avec une fonction propre de
spin α ou β . La forme la plus simple que la fonction d'onde Ψ peut alors prendre est
le produit suivant, dit produit de Hartree :
Y
Ψ= φi (B.4)
i

Toutefois, cette expression ne respecte pas le principe de Pauli : en eet la fonction

d'onde résultante n'est pas antisymétrique par rapport à l'échange de deux électrons.
Une solution, proposée par J. C. Slater, est d'exprimer la fonction d'onde électronique
sous la forme d'un déterminant [46, 47], dit déterminant de Slater :
¯ ¯
¯ ¯
¯ φ1 (1) φ2 (1) . . . φn (1) ¯
¯ ¯
¯ ¯
1 ¯¯ φ1 (2) φ2 (2) . . . φn (2) ¯¯
Ψ= √ ¯ . .. ... .. ¯¯ (B.5)
n! ¯¯ .. . . ¯
¯ ¯
¯ ¯
¯φ1 (n) φ2 (n) . . . φn (n)¯

La résolution de l'équation de Schrödinger (B.1) conduit à la détermination de ces

spin-orbitales φi . Au lieu de rechercher directement ces fonctions, il est plus commode
II Fonctions de Base 35

de les exprimer comme une combinaison linéaire de fonctions de base, connues :

φi (~r) =
X
ci,α χα (~r) (B.6)
α

L'indétermination repose alors sur les coecients du développement c . Il faut néan-

i,α

moins noter que, d'un point de vue mathématique, l'espace vectoriel des fonctions est
de dimension innie. Par conséquent ce développement sur un nombre ni de fonctions
de base  seule possibilité pratique  est une approximation qui sera d'autant plus
justiée que la base sera étendue. On parle alors de complétude de base.

II Fonctions de Base
On peut distinguer trois grands types de fonctions de base pour la description de la
structure électronique. Chacun ayant ses avantages, leurs principales caractéristiques
sont brièvement énoncées par la suite.

II.1 Base de type onde plane

Très utilisées pour l'étude des systèmes en phase condensée, les ondes planes sont
issues de la physique du solide. La périodicité d'un arrangement cristallin produit un
potentiel périodique et impose la même périodicité pour la densité électronique. À
partir du théorème de Bloch et des conditions aux limites périodiques de Born-Von
Karman [48], dénissons une (super-) cellule, selon :

a3

a2

a1
36 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

Cette  boîte , dénie par la matrice h = [~a1 , ~a2 , ~a3 ] a pour volume Ω = det h. Ainsi
le système est invariant selon n'importe quelle translation de vecteur déni selon :

~ = i · ~a1 + j · ~a2 + k · ~a3

L ∀ (i, j, k) ∈ N3 (B.7)

De plus, il est utile de dénir le réseau de l'espace réciproque :

h i
2π(h)−1 = ~b1 , ~b2 , ~b3 (B.8)

tel que ~bi · ~aj = 2πδij . Ainsi le système est invariant dans l'espace réciproque selon
n'importe quelle translation de vecteur dénie par :

~ = i · ~b1 + j · ~b2 + k · ~b3

G (B.9)

Les fonctions de base de type onde plane ont la forme suivante :

1 ~ · ~r)
fG (~r) = √ exp(iG (B.10)
Ω
Les ondes planes possèdent des propriétés remarquables ; en particulier elles sont :
 périodiques par rapport à la cellule h : fG (~r + L)
~ = fG (~r)

 orthonormées :hfG |fG′ i = δG,G′

Il est important de noter que les ondes planes sont des fonctions sans origine, i.e. elles
ne dépendent pas de la position des noyaux. Ceci implique plus particulièrement que
les contributions de Pulay (cf. ŸVII) soient exactement nulles, facilitant ainsi l'évalua-
tion des gradients nucléaires. Ceci implique aussi que ces fonctions soient totalement
délocalisées, ne favorisant pas particulièrement certains atomes ou certaines régions de
l'espace. Ainsi, la manière d'améliorer la qualité de la base est d'augmenter  l'énergie
de coupure , i.e. d'augmenter le plus grand vecteur G
~ à inclure dans le développement

des orbitales. Le grand désavantage des ondes planes est que, du fait de leur nature
non locale, il faut employer un grand nombre de fonctions pour correctement décrire les
zones de l'espace où la densité électronique varie grandement, notamment au voisinage
des noyaux. Cette raison explique l'utilisation systématique de pseudo-potentiels, pour
tous les atomes, lors de calculs en base d'ondes planes, éludant ainsi la description de
la structure nodale complexe des orbitales dans la région de c÷ur.
II Fonctions de Base 37

II.2 Fonctions de Slater et fonctions gaussiennes

D'autres fonctions que les ondes planes peuvent être utilisées pour développer les
orbitales selon la relation (B.6). Un choix possible pour ces fonctions de base sont
les orbitales atomiques. Ce choix se voit justié si l'on considère, qu'au sein d'une
molécule, les atomes gardent en partie leur identité et donc que les molécules sont des
assemblages d'atomes légèrement  distordus . Cette vision des orbitales moléculaires
comme une combinaison nie d'orbitales atomiques est connue comme l'approximation
LCAO (de l'anglais Linear Combination of Atomic Orbital ) [49, 50]. D'un point de vue
mathématique, de nombreuses fonctions peuvent être choisies pour décrire les orbitales
atomiques. En pratique, deux types de fonctions sont communément employés. Les
fonctions de Slater sont caractérisées par une discontinuité à l'origine (prise à la position
du noyau) :
S(~r, ζ) = N~r n−1
exp(−ζ~r)Ylm (θ, φ) (B.11)

où Ylm désigne une harmonique sphérique, N une constante de normalisation, ζ l'ex-

posant de l'orbitale et n, l et m respectivement les nombres quantiques principal, or-
bitalaire et magnétique. Les fonctions de Slater décrivent de manière très satisfaisante
la densité électronique, ainsi que le comportement à l'origine des orbitales. Néanmoins,
l'évaluation des intégrales biélectroniques multicentriques s'avère délicate et coûteuse
en temps de calcul.
Un autre choix de fonctions de base, proposé par S. Boys [51], est l'utilisation de
fonctions gaussiennes, elles-aussi centrées au noyau. La forme cartésienne générale de
ce type de fonctions est :

G(~r) = N xi y j z k exp(−α~r 2 ) (B.12)

où α désigne l'exposant de la fonction et (i, j, k) sont des nombres entiers. Cet exposant
α est à relier directement avec l'extension spatiale de la fonction. Le triplet (i, j, k)
est quant à lui responsable de la nature angulaire de l'orbitale. Plusieurs fonctions
gaussiennes sont nécessaires pour décrire, avec la même qualité qu'une seule fonction
38 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

de Slater, une orbitale atomique. Une orbitale atomique est donc souvent développée sur
plusieurs fonctions gaussiennes. Toutefois les intégrales biélectroniques multicentriques
s'avèrent plus faciles à calculer : en eet le produit de deux fonctions gaussiennes est
une nouvelle fonction gaussienne, l'intégrale peut alors être évaluée analytiquement.

III Méthodes Ab Initio

1
Ex nihilo nihil

III.1 Hartree-Fock  équations de Roothan et Hall

À partir des approximations précédentes, c'est-à-dire l'approximation monoélectro-
nique et l'approximation LCAO, on peut appliquer le principe variationnel et obtenir
les équations dites de Roothan et Hall [52, 53] :
Q
X
∀ i ∈ [1, Q] (Fsr − εi Ssr )cri = 0 (B.13)
r=1

où Q est le nombre de fonctions de base, εi désigne l'énergie orbitalaire, i.e. la va-

leur propre associée à l'orbitale ψi . La reformulation matricielle de ce jeu d'équations
conduit à :
FC = SCε (B.14)

où ε désigne la matrice diagonale des énergies orbitalaires εi , C est la matrice des coe-
cients de la combinaison linéaire des fonctions de base, S est la matrice de recouvrement
et F représente la matrice de Fock. Les éléments de ces matrices s'écrivent :

Srs = hχr |χs i (B.15)

X · ¸
core 1
Frs = Hrs + Pµν (rs|νµ) − (rµ|νs) (B.16)
µν
2
1 Rien (ne vient) de rien. Célèbre aphorisme résumant la philosophie de Lucrèce et d'Epicure, mais

tiré d'un vers de Perse (Satires, III, 84), qui commence par De nihilo nihil (Rien ne vient de rien,

c'est-à-dire Rien n'a été tiré de rien. Rien n'a été créé, mais tout ce qui existe existait déjà en quelque

manière de toute éternité).

III Méthodes Ab Initio 39

Hcore est la matrice de l'hamiltonien de c÷ur, c'est-à-dire l'hamiltonien monoélec-

tronique ne contenant que les termes d'énergie cinétique électronique et d'attraction
nucléaire. P désigne la matrice densité, dénie selon :
occ
X
Pµν = 2 c⋆µi cνi (B.17)
i

et le terme (ij|kl) représente une intégrale de répulsion biélectronique :

Z
1
(ij|kl) = d~r1 d~r2 χ⋆i (~r1 )χj (~r1 ) χ⋆k (~r2 )χl (~r2 ) (B.18)
|~r1 − ~r2 |

La résolution des équations de Roothan et Hall s'eectue de manière itérative. En eet

la matrice de Fock F dépend du résultat  les orbitales  via la matrice densité. À
partir d'un jeu arbitraire de coecients pour la combinaison linéaire des orbitales, la
matrice de Fock est construite et la résolution de l'équation (B.14) conduit à un nouveau
jeu de coecients, point de départ pour un nouveau cycle, et ce jusqu'à convergence.
Cette procédure est la méthode dite du champ auto-cohérent (Self Consistent Field

en anglais). Notons que le coût calculatoire d'un calcul d'énergie dans l'approxima-
tion Hartree-Fock, est proportionnel à Q4 , où Q représente le nombre de fonctions
de base pour le développement des orbitales. Ainsi, dans la méthode Hartree-Fock,
chaque électron interagit avec les autres électrons de façon moyenne, c'est-à-dire qu'il
ressent un champ moyen. Il n'y a pas d'interaction instantanée entre deux électrons.
Ce phénomène est responsable de la corrélation électronique. Les méthodes qui sont
brièvement exposées par la suite sont dites post Hartree-Fock puisqu'elles ont pour but
de décrire  du moins en partie  la corrélation électronique. Ces méthodes  post

Hartree-Fock  ont donc pour objet de calculer l'énergie de corrélation, dénie comme
la diérence entre l'énergie exacte et l'énergie de la limite Hartree-Fock.
40 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

III.2 Méthodes  post Hartree-Fock 

Il existe un grand nombre de méthodes ab initio permettant de prendre en compte
les eets de corrélation électronique. Parmi les approches variationnelles, on peut citer
les méthodes d'interaction de congurations [54] et plus particulièrement l'approche
CASSCF [55] qui envisage une interaction de congurations dans un espace orbita-
laire restreint. On peut aussi citer les approches perturbatives de la corrélation élec-
tronique, telles que les méthodes de Møller-Plesset [56] basées sur un traitement par
la méthode des perturbations de Rayleigh-Schrödinger d'un hamiltonien d'ordre zéro
Hartree-Fock. Il existe aussi des méthodes mixtes qui traitent de manière variation-
nelle les congurations de poids les plus importants et de façon perturbative les autres
congurations, telles que l'approche CIPSI [57, 58, 59] ou bien CASPT2 [60]. Ces ap-
proches basées sur une description multicongurationnelle de la fonction d'onde sont,
entre autres, pertinentes pour la description de systèmes à couche ouverte ou bien
lorsque une quasi-dégénérescence apparaît. On peut enn citer l'approche monocon-
gurationnelle coupled-cluster [61, 62, 63, 64] basée sur une dénition diérente et non
linéaire de la fonction d'onde.
Les méthodes ab initio précédentes ne constituent pas la totalité des approches
existantes pour décrire la structure électronique des édices moléculaires. En eet ces
approches sont basées sur diérents traitements de la fonction d'onde électronique
dépendant des 3n coordonnées électroniques pour pouvoir ainsi accéder à l'énergie
électronique du système. Une alternative consiste à étudier directement la densité élec-
tronique : la théorie de la fonctionnelle de la densité. Notons enn que ces méthodes
ab initio n'ont pas été utilisées de manière combinée avec la dynamique moléculaire
au cours de cette thèse, néanmoins ces approches multicongurationnelles seraient per-
tinentes, par exemple, pour la réalisation de dynamiques moléculaire sur des états
excités.
IV La Théorie de la Fonctionnelle de la Densité 41

IV La Théorie de la Fonctionnelle de la Densité

Les théories de la fonctionnelle de la densité permettent d'obtenir un gain de temps
de calcul substantiel par rapport aux méthodes  post Hartree-Fock  puisque l'opéra-
teur d'échange est modélisé à l'aide d'une densité à une particule. Leur coût calculatoire
est proportionnel à Qm où Q est le nombre de fonctions de base et m = 3 − 4. Cette
théorie a été proposée et utilisée en physique du solide sous l'impulsion de J.C. Slater
qui a proposé de remplacer le terme d'échange Hartree-Fock par le potentiel de Dirac
ρ 3 [15]. En 1964, P. Hohenberg et W. Kohn ont prouvé que les propriétés de l'état
1

fondamental sont dénies de façon biunivoque par la densité électronique ρ(~r) [65].
Ces travaux marquent le début des méthodes issues de la théorie de la fonctionnelle de
la densité.

IV.1 Théorèmes de Hohenberg et Kohn  Approche de Kohn

et Sham
Credo quia absurdum2
Le premier théorème de Hohenberg et Kohn stipule que, pour l'état fondamental,
la densité électronique ρ(~r) détermine le potentiel extérieur [65]. Ce théorème peut
facilement se démontrer ab absurdo. L'énergie de l'état fondamental du système Ee
peut alors s'écrire comme une fonctionnelle de la densité :

Ee [ρ] = T [ρ] + Ve−e [ρ] + Ve−N [ρ] (B.19)

où T [ρ] désigne le terme d'énergie cinétique, Ve−N [ρ] le terme d'interaction électrons-
noyaux et Ve−e [ρ] le terme d'interaction électrons-électrons. Le second théorème de
Hohenberg et Kohn stipule que la densité peut être exactement calculée grâce à un
principe variationnel [65]. Ces équations ne sont pas utilisables directement pour eec-
tuer des calculs atomiques ou moléculaires. Pour faire de la relation formelle (B.19) un
2 Je le crois parce que c'est absurde. Paroles faussement attribuées à Saint Augustin.
42 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

outil pratique, W. Kohn et L.J. Sham ont proposé de décomposer la fonctionnelle de

l'énergie totale en divers termes [66] :
Z
Ee [ρ] = Ts [ρ] + (Vext (~r) + Vcoul (~r)) d~r + Exc [ρ] (B.20)

où Ts [ρ] désigne le terme d'énergie cinétique d'un système ayant la même densité
électronique que le système réel mais au sein duquel il n'y a pas d'interaction entre les
électrons, Vcoul est le terme d'interaction coulombienne entre les électrons :
Z
1 ρ(~r ′ )
Vcoul (~r) = d~r ′
(B.21)
2 |~r − ~r ′ |

Vext est le terme associé au potentiel extérieur des noyaux :

X ZI
Vext (~r) = − (B.22)
~I|
|~r − R
I

et enn Exc qui désigne la fonctionnelle  d'échange-corrélation . Ce terme inclut

toutes les contributions énergétiques non prises en compte par les termes précédents :
l'échange électronique, la corrélation électronique, une partie de l'énergie cinétique ainsi
que des corrections au potentiel classique de Coulomb. Enn Kohn et Sham proposèrent
l'utilisation d'une fonction d'onde déterminantale, constituée de n particules indépen-
dantes dans n orbitales φKS r) [66]. À l'aide de ces orbitales, la densité électronique
i (~

s'écrit :
n
X ¯ KS ¯2
ρ(~r) = ¯φi (~r)¯ (B.23)
i

Ces orbitales sont alors solutions d'un jeu d'équations aux valeurs propres de la forme :
· ¸
1 2
− ∇i + Vext (~r) + Vcoul (~r) + Vxc (~r) φKS r) = εi φKS
i (~ i (~
r) (B.24)
2

où Vxc est le potentiel d'échange corrélation déni comme la dérivée fonctionnelle de

l'énergie d'échange corrélation Exc :
δExc
Vxc (~r) = (B.25)
δρ(~r)

Cette équation aux valeurs propres (B.24) est analogue à celle de la méthode Hartree-
Fock. Cependant, au contraire de l'opérateur de Fock qui est non local, l'opérateur de
IV La Théorie de la Fonctionnelle de la Densité 43

Kohn-Sham est local : il ne dépend que de la position dans l'espace de l'électron ~r et

non de l'électron lui-même. De la même manière que pour l'approche Hartree-Fock, ces
orbitales sont développées sur un jeu de fonctions de base, et ces équations sont résolues
de manière itérative par un processus auto-cohérent. Il faut enn noter que la densité,
et donc l'énergie, seront déterminées exactement dès lors que le potentiel d'échange-
corrélation exact est connu. La principale diculté des méthodes de la théorie de la
fonctionnelle de la densité réside dans la dénition de cette fonctionnelle d'échange-
corrélation. Dans ce but, plusieurs approximations ont été réalisées pour conduire à
diérentes fonctionnelles.

IV.2 Les diérentes fonctionnelles

IV.2.a Le gaz uniforme d'électrons : l'approximation de la densité locale

Les premières mises en ÷uvre de l'approche de Kohn et Sham utilisent pour la

fonctionnelle d'échange-corrélation l'approximation de la densité locale (LDA, de l'an-
glais Local Density Approximation ). Celle-ci est basée sur l'analyse du gaz homogène
d'électrons. Ces fonctionnelles sont connues pour sous-estimer l'énergie d'échange et
surestimer l'énergie de corrélation. D'un point de vue structural, ces fonctionnelles
conduisent, en général, à des longueurs de liaisons trop courtes, une surestimation des
énergies de liaisons et des liaisons hydrogènes trop faibles.

IV.2.b L'approximation du gradient généralisé

L'approximation de la densité locale est appropriée à l'étude des systèmes dont la

densité varie lentement dans l'espace. Une façon d'améliorer la fonctionnelle d'échange-
corrélation est de la rendre dépendante à la fois de la densité locale et des variations
de cette dernière, c'est-à-dire du gradient de la densité :
Z
Exc = ρ(~r)Vxc [ρ(~r), ∇ρ(~r)] d~r (B.26)
44 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

La plupart des fonctionnelles corrigées du gradient (GGA, de l'anglais Generalised

Gradient Approximation ) est construite comme l'addition d'une correction à une fonc-
tionnelle LDA. L'une des plus populaires fonctionnelles d'échange est celle proposée
par Becke en 1988 [67] ; celle-ci est paramétrée sur les énergies d'échanges connues
d'atomes de gaz rares. Cette dernière est très souvent associée à la fonctionnelle de
corrélation proposée par C. Lee, W. Yang et R.G. Parr [68], paramétrée sur l'énergie
de corrélation de l'atome d'hélium. Il existe bien d'autres fonctionnelles corrigées du
gradient, obtenues de diérentes manières, comme l'ajustement des paramètres sur des
données moléculaires par exemple ou bien la régularisation formelle du comportement
asymptotique. Les fonctionnelles corrigées du gradient s'avèrent beaucoup plus ecaces
pour les calculs moléculaires. Toutefois des problèmes persistent sur des données éner-
gétiques, telles que les énergies de liaisons ou bien les énergies d'activation des réactions
chimiques.

IV.2.c Les fonctionnelles hybrides

Les fonctionnelles hybrides incluent dans leur forme une fraction de l'échange exact,
comparable à l'échange Hartree-Fock mais calculé à l'aide des orbitales Kohn-Sham.
L'une des plus utilisée est la fonctionnelle proposée par A.D. Becke [69], ajustée sur des
énergies d'atomisation via trois paramètres. Les performances de ce type de fonction-
nelles sont bonnes, ce qui font d'elles les plus utilisées en chimie. Toutefois, il faut noter
que le calcul de l'échange exact augmente signicativement le temps de calcul, ce qui
rend l'utilisation de ce type de fonctionnelle prohibitive lorsque la base est développée
sur des ondes planes.
Le principal avantage des méthodes basées sur la théorie de la fonctionnelle de
la densité est un gain substantiel de temps de calcul par rapport aux méthodes ab

initio. Outre les diérentes approximations pour l'hamiltonien décrivant la structure

électronique, il existe d'autres outils permettant un gain de temps de calcul conséquent,
notamment en réduisant la taille de la base : les pseudo-potentiels.
V Les Pseudo-Potentiels 45

V Les Pseudo-Potentiels

Ambitiosa recidet ornamenta3

V.1 Les pseudo-potentiels atomiques

Notons que les propriétés énoncées par la suite à propos des pseudo-potentiels s'ap-
pliquent dans le cas où la base atomique est développée sur des fonctions gaussiennes.
Même si les notions suivantes sont analogues à celles employées pour les pseudo-
potentiels dans une base d'ondes planes, certaines propriétés sont diérentes.

Les similitudes chimiques et physiques des éléments classés dans le tableau de Men-
deleïev sont dues à la structure électronique de valence : ce sont les électrons de valence
qui déterminent les propriétés chimiques et physiques des atomes et des molécules. Les
électrons de c÷ur ne sont que très légèrement aectés par l'environnement molécu-
laire. C'est pourquoi des approximations où le c÷ur des atomes demeure invariant
dans l'environnement moléculaire ont été imaginées. Une de ces approximations repose
sur une théorie introduite en physique du solide en 1935 et propose l'utilisation de
pseudo-potentiels [70, 71], en vue de simuler les eets des électrons de c÷ur tout en
préservant les propriétés des électrons de valence. Cette méthode commencera à être
véritablement utilisée en chimie quantique qu'à partir des années 1970 [72]. Dans cette
approche, les électrons de c÷ur ne sont plus traités explicitement lors du calcul mais
leur présence est simulée à l'aide d'un pseudo-potentiel tel que les électrons de valence
aient le même comportement, dans le champ du pseudo-potentiel que dans le champ
réel. Les principaux avantages de cette méthode sont :

 la réduction du nombre d'électrons à traiter lors du calcul,

 la réduction du nombre de fonctions de base et ainsi un allègement du calcul mais
aussi la possibilité d'élargir la base de valence,
 la possibilité d'inclure certains eets relativistes.

Les pseudo-potentiels de c÷ur ne conservent pas la structure nodale des orbitales de

3 Il retranchera les ornements ambitieux. Précepte d'Horace, dans l'Art poétique (v. 447).
46 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

valence dans la région de c÷ur. Le terme monoélectronique de l'hamiltonien pour un

électron de valence devient :

1 z
hi (ri ) = − ∇2i − + Wps,i (B.27)
2 ri

où z = Z − nc désigne la charge nette du c÷ur et Wps,i l'opérateur pseudo-potentiel.

Cet opérateur pseudo-potentiel va directement agir sur les électrons de valence pour
que leur comportement mime celui qu'ils auraient dans le champ réel uniquement dans
la zone de valence. On peut distinguer deux grandes classes de pseudo-potentiel : les
pseudo-potentiels cohérents de forme (shape-consistent en anglais) et ceux cohérents

en énergie (energy-consistent en anglais). La distinction provient de la méthode d'ex-

traction du pseudo-potentiel où, dans le premier cas, on cherche avant tout à reproduire
la forme des orbitales dans la région de valence, alors que, dans le second cas, on cherche
avant tout à reproduire certaines observables telle que l'énergie orbitalaire.

V.2 Les potentiels eectifs de groupe

Dans la grande majorité des réactions chimiques, seuls quelques atomes de l'édice
moléculaire interviennent directement dans le processus réactionnel, le reste de la mo-
lécule étant principalement  spectateur . Diérencier les électrons actifs et inactifs
d'un système moléculaire est une approximation majeure en chimie quantique. Dans le
cadre de calculs ab initio ou DFT, les électrons de c÷ur sont bien souvent remplacés
par un pseudo-potentiel. Dans ce cadre, la séparabilité des électrons de c÷ur et de
valence est justiée selon deux critères : les diérences spatiales et énergétiques entre
ces électrons. Toutefois, ce concept de séparabilité n'est pas exclusif aux électrons de
c÷ur, il peut être généralisé au niveau d'un groupement chimique, diérenciant ainsi
deux parties, l'une  active  et la seconde  spectatrice . Un pseudo-potentiel de
groupe est donc l'analogue, pour un groupement chimique, à un potentiel eectif de
c÷ur, pour l'atome. Néanmoins le critère de séparabilité énergétique n'est plus ap-
plicable : cette utilisation s'appuie uniquement sur la séparabilité spatiale. Les bases
VI Approches Hybrides : Les Méthodes  QM/MM  47

théoriques des potentiels eectifs de groupe ont été établies dès 1987 [73], développées
par la suite [74] et leur pertinence a été démontrée dans diverses études [75, 76, 77].
Les principaux avantages des pseudo-potentiels de groupe sont similaires à ceux des
pseudo-potentiels atomiques, c'est-à-dire la réduction du nombre d'électrons explicites
lors du calcul ainsi que la réduction du nombre de fonctions de base. Par contre, un
nouvel avantage de cet outil est son utilisation pour la dénition et le traitement de la
frontière dans les méthodes mixtes de type QM/MM [78]. Ce type de méthode hybride,
qui a aussi pour but la réduction du temps de calcul, est brièvement exposé par la suite.

VI Approches Hybrides : Les Méthodes  QM/MM 

Alors que les champs de force de la mécanique moléculaire sont intrinsèquement
incapables de décrire la structure électronique des molécules, le traitement par des
méthodes de mécanique quantique des systèmes de grande taille est impraticable, mal-
gré les progrès des algorithmes à croissance linéaire [79, 80]. Cependant, les processus
chimiques concernent le plus souvent un nombre réduit d'atomes, le centre actif. Le
reste du système, constituant l'environnement des atomes actifs, a tout de même une
inuence essentielle. Les méthodes hybrides de type  QM/MM , dont les premiers
calculs furent introduits dès 1976 [81], proposent de traiter le centre actif grâce aux mé-
thodes de mécanique quantique et de décrire l'environnement par une méthode moins
coûteuse comme la mécanique moléculaire. Deux grands schémas d'approches hybrides
peuvent être distingués : les schémas additifs et les schémas soustractifs.
Les méthodes relevant d'un schéma additif envisagent l'écriture de l'hamiltonien du
système S selon :

H (S) = HQM (I) + HQM/M M (I/O) + HM M (O) (B.28)

où I désigne l'ensemble des particules décrites à haut niveau de la région QM, O pour
les particules de la région MM décrites à un plus bas niveau de théorie et (I/O) les
particules communes aux deux domaines. L'hamiltonien d'interaction HQM/M M entre
48 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

les deux régions peut s'exprimer simplement dans le cas d'un découpage des domaines
QM et MM sans  rupture  de liaisons. Par exemple un système soluté/solvant peut
être décrit en traitant à haut niveau le soluté, à un plus bas niveau de théorie le solvant ;
l'hamiltonien d'interaction entre les deux régions décrit alors l'interaction électrosta-
tique entre des charges ponctuelles de la partie solvant et la structure électronique du
soluté, polarisant ainsi cette dernière. L'expression de l'interaction entre les deux do-
maines s'avère plus délicate si la frontière est intramoléculaire. Il est alors indispensable
de prendre en compte les interactions covalentes. Ceci a fait l'objet de nombreuses pro-
positions, telles que la méthode LSCF [82] (de l'anglais Local Self Consistent Field ).

Ces diérents traitements ne seront pas plus développés ici.

La méthode illustrant un schéma soustractif pour les approches de type QM/MM
est la méthode ONIOM [83, 84] (de l'anglais Our N-layered Integrated molecular Or-

bital/molecular Mechanics approach ). En conservant les mêmes notations que précé-

demment, ce schéma soustractif est déni par l'hamiltonien suivant :

H (S) = HQM (I) + HM M (S) − HM M (I) (B.29)

Ici, un hamiltonien d'interaction entre les deux domaines n'est pas dénie ; ONIOM est
une méthode d'extrapolation. On réalise ainsi plusieurs calculs successifs : le système
modèle à un  haut niveau  de théorie (HQM (I)), puis deux calculs à un niveau de
théorie moindre, l'un pour le système modèle (HM M (I)) et l'autre pour le système réel
(HM M (S)). Les énergies et leurs dérivées sont ensuite combinées pour engendrer une
surface d'énergie potentielle associée au système entier. ONIOM est la généralisation
de la méthode  IMOMM/IMOMO  [85, 86], laquelle est limitée à une partition du
système en deux couches. ONIOM étend ce concept à n couches, les calculs avec deux
ou trois couches étant, en pratique, les plus répandus.
Avec l'ensemble des méthodes décrites précédemment, nous avons à disposition
les outils nécessaires pour déterminer l'énergie potentielle d'un système moléculaire.
An de réaliser le traitement dynamique du système, les gradients nucléaires sont
indispensables et peuvent être dérivés de l'énergie potentielle.
VII Les Dérivées Premières de l'Énergie 49

VII Les Dérivées Premières de l'Énergie

Pour la réalisation de dynamique moléculaire classique, les dérivées premières de

l'énergie par rapport aux positions nucléaires sont nécessaires. Celles-ci jouent le rôle
de l'accélération des noyaux. Le but de ce paragraphe est d'expliquer, succinctement,
comment ces dérivées sont calculées.
Une approche simpliste pour calculer les dérivées premières de l'énergie serait de
diérencier numériquement l'énergie. Il est cependant possible d'obtenir ces gradients
analytiquement. Leur expression dans l'approximation Hartree-Fock a initialement été
proposée par P. Pulay en 1969 [87, 88]. Les idées essentielles de ce développement sont
illustrées par la suite. Notons que des expressions analogues sont obtenues dans le cadre
de la théorie de la fonctionnelle de la densité.
L'énergie totale d'un système  à couche fermée  dans l'approximation Hartree-
Fock peut s'exprimer selon :
X 1X
Ee = core
Pµν Hµν + Pνµ Pλσ (µν||σλ) + VN −N (B.30)
µν
2 µνλσ

où VN −N représente le terme de répulsion internucléaire déni par :

XX Z I ZJ
VN −N = (B.31)
~I − R
|R ~J|
I J>I

et où les intégrales sont notées simplement :

1
(µν||σλ) = (µν|σλ) − (µλ|σν) (B.32)
2

La dérivée partielle de l'énergie dénie par la relation (B.30) par rapport à la position
du noyau R
~ I conduit à :
X core
∂Ee ∂Hµν 1X ∂(µν||σλ) ∂VN −N
= Pµν + Pνµ Pλσ +
~I
∂R ∂ ~I
R 2 ∂ ~I
R ~I
∂R
µν µνλσ
X ∂Pνµ X ∂Pνµ (B.33)
core
+ Hµν + P (µν||σλ)
~ ~ λσ
µν ∂ RI µνλσ ∂ RI

où Pνµ désigne un élément de la matrice densité, déni selon la relation (B.17). La

relation (B.33) suggère que les dérivées des coecients du développement LCAO sont
50 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

requises. Toutefois, le développement des deux derniers termes de la relation (B.33)

donne :
X ∂Pνµ X ∂Pνµ
core
Hµν + P (µν||σλ) = (B.34)
~I
∂R ~ I λσ
∂R
µν µνλσ
X X ∂Cµi
4 core
Hµν Cνi (B.35)
~
µν i ∂ RI
X X ∂Cµi
+4 P (µν||σλ)Cνi (B.36)
~ λσ
µνλσ i ∂ RI
" #
X X ∂Cµi X
core
=4 Hµν + Pλσ (µν||σλ) Cνi
~
∂ RI
µν i λσ

(B.37)
X X ∂Cµi
=4 F C (B.38)
∂ ~ I µν νi
R
µν i
X X ∂Cµi
=4 εi S C (B.39)
~ µν νi
i µν ∂ RI

Pour évaluer les dérivées de ces coecients, la condition d'orthonormalité des orbitales
est employée :
X
Cµi Sµν Cνj = δij (B.40)
µν

La diérenciation de cette condition conduit à :

X ∂Cµi X ∂Sµν
2 Sµν Cνi = − Cµi Cνi (B.41)
~I
∂R ~I
∂R
µν µν

À l'aide de cette expression et de la relation (B.33), la dérivée partielle de l'énergie

totale peut s'exprimer selon :
X core
∂Ee ∂Hµν 1X ∂(µν||σλ)
= Pνµ + Pνµ Pλσ
~I
∂R ∂R~I 2 µνλσ ∂R~I
µν
(B.42)
∂VN −N X ∂Sµν
+ − Qµν
∂R ~I ∂R~I
µν

où Qµν est déni par :

X
Qµν = 2 εi Cµi Cνi (B.43)
i
VII Les Dérivées Premières de l'Énergie 51

Il faut remarquer ici que, lors du développement entre les relations (B.38) et (B.39),
les orbitales moléculaires sont supposées être solutions des équations de Roothan. Ceci
n'est vrai qu'à la condition d'avoir minimisé l'énergie par rapport à ces orbitales. Cette
condition n'est, a priori, pas remplie lors d'une dynamique Car-Parrinello. L'expression
de ces dérivées [89] est dans ce cas :
X core
∂Ee ∂Hµν 1X ∂(µν||σλ)
= Pνµ + Pνµ Pλσ
~I
∂R ∂R~I 2 µνλσ ∂R~I
µν
X ∂Mµσ (B.44)
∂VN −N
+ +4 u F c
∂R ~I ∂ ~ I σi µν νi
R
µνσi

Cette expression dière de l'expression (B.42) par le dernier terme où apparaît la

dérivée de la transformation en base orthogonale des orbitales, dénie par :

X¡ ¢
uσi = M−1 σν
cνi (B.45)
ν

Dans le cas d'une transformation symétrique [90] (appelée aussi transformation de

Löwdin), la matrice M est :
¡ ¢ X
Mµν = S−1/2 µν = Tµα s−1/2
α Tνα (B.46)
α

où sα et Tµα désignent respectivement les valeurs et vecteurs propres de la matrice de

recouvrement S. La dérivée de cette matrice de passage s'écrit alors :
¡ ¢ ∂S
∂ S−1/2 µν X −Tµα ∂ R~αβ Tνβ
= ³ I ´ (B.47)
~I
∂R 1/2 1/2 1/2
s α sβ s α + s β
1/2
αβ

Les dérivées premières de l'énergie Hartree-Fock peuvent donc être calculées analyti-
quement, selon la relation (B.42) si l'énergie électronique est minimisée par rapport
aux orbitales moléculaires. Ce sera le cas lors d'une dynamique de type  Born-
Oppenheimer , mais lors d'une dynamique  Car-Parrinello , ces gradients seront
calculés selon la relation (B.44).
Dans ces dernières relations apparaissent les dérivées des intégrales électroniques,
il faut donc préciser la dérivée des fonctions de base. Rappelons l'expression générale
52 LA STRUCTURE ÉLECTRONIQUE DES ÉDIFICES MOLÉCULAIRES

des fonctions gaussiennes :

Gijk (~r) = Nijk xi y j z k exp(−α~r 2 ) (B.48)

où ~r = (x, y, z) désigne la position relative d'un électron par rapport au noyau sur
lequel est centrée cette fonction gaussienne, i.e. ~ I . La position de ce noyau
~r = ~ri − R
est notée R
~ I = (XI , YI , ZI ), l'expression de la dérivée de cette fonction, par exemple

par rapport à XI , est :

∂Gijk (~r)
= 2α Gi+1,jk (~r) − i Gi−1,jk (~r) (B.49)
∂XI

Ainsi, la dérivée  par rapport à X  d'une orbitale de type px (i = 1, j = k = 0)

conduira à deux orbitales, l'une de type s et l'autre de type dx2 .
Bien que dans les relations (B.42) et (B.44) il n'apparaisse pas explicitement de
dérivée de la matrice densité, il ne faut néanmoins pas confondre ces relations avec le
théorème d'Hellmann-Feynman. Notons l'énergie totale électronique selon :

Ee = hΨ|H |Ψi (B.50)

la dérivée première de l'énergie par rapport à la position du noyau I s'écrit alors :

¿ À ¿ ¯ ¯ À ¿ À
∂Ee ∂Ψ ¯ ∂H ¯ ∂Ψ
= ¯
|H | Ψ + Ψ ¯ ¯ Ψ + Ψ |H | (B.51)
~I
∂R ~I
∂R ~I ¯
∂R ~I
∂R

La condition de Hellmann-Feynmann est la suivante :

¿ À ¿ À
∂Ψ ∂Ψ
|H | Ψ + Ψ |H | =0 (B.52)
~I
∂R ~I
∂R

Ces contributions sont souvent appelées termes de Pulay. Cette relation n'est vériée
que si la fonction d'onde Ψ est exacte. Ceci conduit alors à l'expression suivante pour
les gradients nucléaires :
¿ ¯ ¯ À
∂Ee ¯ ∂H ¯
= ¯
Ψ¯ ¯Ψ
¯ (B.53)
~I
∂R ~
∂ RI
En pratique, pour des calculs moléculaires, la fonction d'onde Ψ obtenue n'est jamais
la fonction d'onde exacte. Le théorème d'Hellmann-Feynmann est donc caduque et ne
VIII Conclusion 53

doit pas être utilisé. Il peut néanmoins être utilisé si la base est développée sur des
fonctions d'ondes planes. En eet dans ce cas, même si la fonction d'onde totale n'est
pas la solution exacte, la condition (B.52) est vériée puisque les fonctions de base ne
dépendent pas de la position des noyaux.

VIII Conclusion
Il existe diérentes méthodes, plus ou moins précises, pour décrire, de manière quan-
tique, la structure électronique des édices moléculaires. Néanmoins, le choix d'une
méthode pertinente pour décrire un problème particulier est à concilier avec son coût
calculatoire. Une approche de dynamique moléculaire devient très rapidement  gour-
mande  en temps de calcul. Par conséquent, l'utilisation des méthodes de la théorie de
la fonctionnelle de la densité en dynamique moléculaire est bien souvent le compromis
le plus judicieux, à condition que la structure électronique du problème étudié ne pré-
sente pas de caractère multicongurationnel trop marqué. De plus la taille des systèmes
moléculaires qui peuvent être décrits de manière quantique reste encore relativement
modeste. Les méthodes hybrides, de type QM/MM peuvent s'avérer pertinentes pour
l'étude de système de plus grande taille. Une méthode permettant de réduire, de ma-
nière signicative, le temps de calcul dans les approches de dynamique moléculaire ab
initio est l'approche de Car et Parrinello [7]. En eet cette méthode permet d'éviter la
procédure de minimisation de l'énergie par rapport aux orbitales moléculaires en propa-
geant, de manière classique, la fonction d'onde. Une telle approche permet donc un gain
de temps de calcul signicatif par rapport à la dynamique de type Born-Oppenheimer
dans laquelle l'étape limitante est cette minimisation de l'énergie potentielle. Dans
la méthode Car-Parrinello, l'étape limitante du point de vue du temps de calcul est
alors l'estimation des gradients de l'énergie potentielle par rapport aux positions nu-
cléaires. Dans la partie suivante, l'attention est portée sur la dynamique moléculaire
Car-Parrinello et plus particulièrement sur la variante développée durant cette thèse.
C Dynamique Moléculaire

Car-Parrinello

Le but de cette partie est de détailler la mise en ÷uvre pratique de l'approche  Car-

Parrinello , en particulier la variante qui a été développée durant cette thèse. Ces

diverses approches se distinguant principalement par le type de fonctions de base

utilisées, le principe général de la méthode est tout d'abord rappelé.

D epuis l'article de Car et Parrinello [7] proposant la propagation classique de

la fonction d'onde, leur approche a fait l'objet de plusieurs implémentations,
variant notamment par le type de fonctions de base utilisées pour la description de la
structure électronique. Toutefois, aujourd'hui, la dynamique moléculaire Car-Parrinello
est synonyme d'un développement de la fonction d'onde au moyen d'ondes planes. Pour
comprendre pourquoi cette  suprématie  est si agrante, un bref rappel historique
est nécessaire.

I Aspects Historiques
Le développement initial de l'approche Car-Parrinello a été proposé, en 1985, dans
le cadre de la théorie de la fonctionnelle de la densité, dans un schéma de description
de la fonction d'onde par des ondes planes [7]. Ce type de fonctions de base s'avère
notamment pertinent pour l'étude de systèmes en phase condensée. En eet, ces fonc-
tions ne sont pas rattachées à un centre particulier, comme les noyaux dans le cas de
56 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

fonctions gaussiennes ; ces fonctions sont ainsi délocalisées et elles ne favorisent pas
particulièrement une zone de l'espace par rapport à une autre. Leur utilisation s'avère
donc  naturelle  pour décrire des systèmes où, a priori, la densité électronique est
bien délocalisée dans l'espace, voire périodique, comme dans le cas d'un solide cristallin.
Toutefois, un grand nombre d'ondes planes est nécessaire pour décrire correctement les
fortes variations spatiales de densité électronique. Ceci est manifeste notamment pour
décrire la structure nodale des orbitales de valence dans la région de c÷ur, cepen-
dant l'utilisation de pseudo-potentiels permet de réduire signicativement le nombre
de fonctions d'ondes planes. D'un autre côté, de fortes inhomogénéités électroniques
peuvent apparaître dans la région de valence des atomes. C'est notamment le cas lors
d'une réaction chimique, où des liaisons se rompent et d'autres se forment, ou bien
lorsque ces liaisons possèdent un caractère ionique marqué. Un eort calculatoire sera
donc nécessaire pour correctement décrire, à l'aide d'ondes planes, ce type de système.
L'utilisation d'une base locale semble alors plus pertinente pour étudier des systèmes
moléculaires où de fortes inhomogénéités électroniques apparaissent. En eet, même si
les intégrales sont plus coûteuses en terme de temps de calcul, le nombre de fonctions
de base nécessaire pour décrire correctement la structure électronique est bien moindre.

L'utilisation de ce type de base pour la réalisation de simulations de dynamique

moléculaire dans le cadre d'un lagrangien étendu, i.e. de type Car-Parrinello, a été
proposé dès 1990 par M. J. Field [9, 24]. Ces travaux, visant des simulations de  recuit
simulé  ou de minimisation directe, ont été réalisés dans un cadre semi-empirique ainsi
que dans l'approximation Hartree-Fock pour la description de la structure électronique.
Très rapidement, l'idée de combiner la dynamique moléculaire Car-Parrinello avec une
description en base locale de la structure électronique a été reprise par E. A. Carter [25].
Ce schéma fut appliqué dans des approches mono- et multi-congurationnelles de la
fonction d'onde (Hartree-Fock, Generalized Valence Bond et interaction de congura-
tions) [25, 26, 27, 28, 29, 30, 31, 32]. Après avoir utilisé à maintes reprises ce schéma, ces
travaux conclurent quant à l'inecacité de la dynamique moléculaire Car-Parrinello en
II Quelques Aspects Techniques 57

base locale [31]. En eet, l'approche proposée sourait d'un épineux problème : l'éner-
gie totale du système n'était pas conservée au cours de la simulation, entraînant des
 pertes  de l'ordre de 10−2 u.a. en 1 ps de simulation [30]. Ce fait assez surprenant
peut notamment s'expliquer puisque les termes  de Pulay  (cf. Ÿ VII p. 49), non
négligeables pour le calcul des gradients nucléaires, avaient été omis. Ceci est d'au-
tant plus surprenant puisque le schéma proposé par M. J. Field tenait compte de ces
termes : aucune composante de ces gradients n'était oubliée, en particulier les termes
associés au fait que les orbitales propagées ne minimisent pas l'énergie électronique.
De plus, d'autres termes intervenant dans les équations de mouvement associés aux
noyaux avaient été omis : les termes associés aux contraintes d'orthonormalisation (cf.
ŸII.3 éq. A.38 p. 15). Depuis ces développements, l'approche de Car et Parrinello en
base locale était considérée comme vaine et sans avenir. Toutefois H. B. Schlegel et al.

ont récemment proposé un schéma de dynamique moléculaire en base locale à l'aide

d'un lagrangien étendu [33]. Cette méthode est brièvement discutée par la suite ; néan-
moins, auparavant, le principe général et les détails de la mise en ÷uvre de l'approche
Car-Parrinello sont rappelés.

II Quelques Aspects Techniques

Les techniques de mise en ÷uvre de la méthode Car-Parrinello ont été principale-
ment développées par M. E. Tuckerman et M. Parrinello [91]. Notamment, ces derniers
détaillèrent l'intégration des équations du mouvement à l'aide de l'algorithme de Ver-
let aux vitesses. Ce paragraphe reprend donc les principaux points de l'intégration des
équations du mouvement Car-Parrinello.

II.1 Les équations du mouvement

La méthode Car-Parrinello, initialement développée dans le cadre de la théorie de

la fonctionnelle de la densité, est basée sur un lagrangien étendu, où les électrons,
58 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

représentés par un jeu d'orbitales {ψi (~r)}, exécutent une dynamique classique ctive,
leur permettant de  suivre  le mouvement nucléaire. Les équations du mouvement du
système dynamique complet sont dérivées du lagrangien postulé par Car et Parrinello,
dont l'expression est rappelée :
Z
1 X 1X ~˙ 2
LCP = µ d~r|ψ̇i (~r)|2 + MI R I
2 i 2 I
X µZ ¶ h i (C.1)
+ Λij ∗ ~I}
d~rψi (~r)ψj (~r) − δij − Ee {ψi }, {R
i,j

On retrouve dans ce lagrangien les termes habituels de la dynamique moléculaire

classique, i.e. les termes d'énergie cinétique nucléaire et d'énergie potentielle, auxquels
viennent s'ajouter un terme d'énergie cinétique des orbitales moléculaires ainsi qu'un
terme de contrainte, permettant d'assurer l'orthonormalité de ces orbitales :
Z
d~rψi⋆ (~r)ψj (~r) = δij (C.2)

Les équations du mouvement associées à ce lagrangien sont obtenues à l'aide des rela-
tions d'Euler-Lagrange et conduisent aux relations suivantes :

MI R ~¨ I = − ∂Ee (C.3)
∂R~I
δEe X
µψ̈i (~r, t) = − ⋆ + Λij ψj (~r, t) (C.4)
δψi j

La quantité δEe
δψi⋆
peut être ré-écrite de manière équivalente :

δEe
⋆
= −fi HeKS ψi (C.5)
δψi

où HeKS désigne l'hamiltonien Kohn-Sham.

Cette approche découle de la mécanique classique, plusieurs invariants peuvent donc
être dénis, notamment l'énergie totale du système, qui doit être conservée au cours
d'une simulation :
XZ 1X h i
Et = µ d~r|ψ̇i (~r)|2 + ~˙ 2 + Ee {ψi }, {R
MI R ~I} (C.6)
I
i
2 I
II Quelques Aspects Techniques 59

On retrouve dans cette expression l'énergie totale réelle, c'est-à-dire l'énergie cinétique
nucléaire et l'énergie potentielle, à laquelle s'ajoute le terme d'énergie cinétique ctive
associée aux orbitales moléculaires. La conservation de l'énergie totale réelle représente
une mesure de la qualité de la trajectoire, puisque l'énergie totale réelle est l'invariant
associé à une dynamique Born-Oppenheimer de référence. Ces deux quantités seront
ainsi fondamentales pour décrire l'ecacité et la pertinence de la méthode. L'intégra-
tion des équations du mouvement (C.3) et (C.4) au moyen d'un propagateur ne pose
pas de diculté particulière, seule la détermination des forces associées aux contraintes
d'orthonormalisation est particulière.

II.2 Les contraintes d'orthonormalisation

La méthode pour satisfaire les contraintes d'orthonormalisation des orbitales mo-

léculaires (C.2) dépend de l'algorithme choisi pour intégrer les équations du mouve-
ment (C.3) et (C.4). L'algorithme le plus utilisé est le propagateur de Verlet. Le détail
de l'intégration des équations contraintes du mouvement électronique est précisé par
la suite pour les algorithmes Verlet et Verlet aux vitesses.

II.2.a Verlet

L'approche exposée par la suite pour résoudre l'intégration contrainte des équations
du mouvement des orbitales est inspirée de la proposition de J. P. Ryckaert [92]. Dans
cette procédure, l'algorithme de Verlet prédit l'évolution non contrainte des orbitales
selon :
(δt)2
|ψei (t + δt)i = 2|ψi (t)i − |ψi (t − δt)i + |Φi (t)i (C.7)
µ
où |Φi (t)i = −fi HeKS |ψi (t)i est la  force  non contrainte au temps (t) agissant sur
l'orbitale ψi . Ces orbitales partiellement prédites sont ensuite corrigées en ajoutant les
forces associées aux contraintes d'orthonormalité :
X
|ψi (t + δt)i = |ψei (t + δt)i + Xij |ψj (t)i (C.8)
j
60 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

(δt)2
où Xij = µ
Λij . Les multiplicateurs de Lagrange inconnus sont déterminés en appli-
quant la condition d'orthonormalité (C.2) aux orbitales au temps (t + δt) : la substi-
tution de la relation (C.8) dans l'équation (C.2) conduit à une relation que la matrice
X doit satisfaire :
XX† + XB + B† X† = I − A (C.9)

où Aij = hψei (t+δt)|ψej (t+δt)i et Bij = hψi (t)|ψej (t+δt)i. En notant que A = I+O(δt2 )
et B = I + O(δt), l'équation (C.9) peut être résolue de manière itérative selon :
1
X(n+1) = [I − A + X(n) (I − B) + (I − B† )X(n) − X2(n) ] (C.10)
2

en partant d'un jeu initial X(0) = 12 (I − A). L'équation (C.9) peut ainsi être résolue
avec une tolérance de 10−6 en 4 à 6 itérations [91]. La matrice X obtenue est alors
utilisée pour obtenir les orbitales au temps (t + δt) selon la relation (C.8).

II.2.b Verlet aux vitesses

Une autre méthode, équivalente pour résoudre des éventuelles contraintes est l'al-
gorithme  RATTLE  [93], basé sur le propagateur de Verlet aux vitesses. Dans cette
approche, la contrainte (C.2) et sa dérivée première par rapport au temps sont satis-
faites :
hψ̇i (t)|ψj (t)i + hψi (t)|ψ̇j (t)i = 0 (C.11)

Il faut noter que cette condition supplémentaire par rapport au cas précédent est en
principe nécessaire pour démontrer une stricte conservation de l'énergie totale à partir
des équations du mouvement. En pratique toutefois, en utilisant l'algorithme de Verlet
original, la relation (C.2) est satisfaite à la précision numérique près. Dans l'algorithme
Verlet aux vitesses, les orbitales |ψi (t)i et leurs vitesses associées |ψ̇i (t)i sont d'abord
corrigées simultanément par les forces |Φi (t)i :
(δt)2
|ψei (t + δt)i = |ψi (t)i + δt|ψ̇i (t)i + |Φi (t)i (C.12)
2µ
δt
|ψėi (t + δt)i = |ψ̇i (t)i + |Φi (t)i (C.13)
2µ
II Quelques Aspects Techniques 61

Les orbitales |ψi (t + δt)i sont ensuite obtenues en suivant la même procédure que dans
le cas de l'algorithme de Verlet original, dénie par les relations (C.8), (C.9) et (C.10),
excepté le fait que la matrice X est ici dénie par Xij = (δt)2
2µ
Λij . Une fois que les
multiplicateurs de Lagrange Λij sont connus, les vitesses des orbitales sont prédites
selon :

′ ė (t + δt)i + δt |Φ (t + δt)i + δt X Λ |ψ (t)i

|ψ̇i (t + δt)i = |ψ (C.14)
i i ij j
2µ 2µ j

où les forces |Φi (t + δt)i représente la force au temps (t + δt) calculée à partir des or-
bitales |ψi (t + δt)i. La correction nale apportée aux vitesses provient des contraintes
d'orthonormalisation à (t + δt) en utilisant un jeu diérent de multiplicateurs de La-
grange Λ′ij , assurant que la relation (C.11) est satisfaite au temps (t + δt). Les vitesses
des orbitales sont ainsi ré-écrites selon :

′ X
|ψ̇i (t + δt)i = |ψ̇i (t + δt)i + Yij |ψj (t + δt)i (C.15)
j

où Yij = δt ′
Λ
2µ ij
. En substituant l'expression (C.15) dans la relation (C.11) exprimée au
temps (t + δt), les nouveaux multiplicateurs de Lagrange sont dénis par :

1¡ ¢
Y=− C + C† (C.16)
2
′
où Cij = hψi (t + δt)|ψ̇j (t + δt)i.
Finalement, en utilisant l'expression des multiplicateurs de Lagrange (C.16) dans l'ex-
pression (C.15), on obtient les vitesses des orbitales corrigées |ψ̇i (t + δt)i. Le fait que
la matrice Y peut être obtenue sans l'aide d'un schéma itératif montre que les vi-
tesses des orbitales vérient exactement la contrainte d'orthonormalisation éq. (C.11).
Bien que l'algorithme de Verlet aux vitesses semble plus contraignant et plus coûteux,
son utilisation est préférée par rapport aux autres algorithmes puisque l'incorporation
d'autres outils, tels que les thermostats ou les contraintes, est plus aisée. Le détail de
l'implémentation de thermostats ou de contraintes géométriques sera précisé plus tard
dans ce manuscrit.
62 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

III Propagation Classique de la Matrice Densité

Un schéma de dynamique moléculaire, de type Car-Parrinello, en base locale, a
récemment été proposé par H. B. Schlegel et al. [33]. À la diérence de la méthode
originale de Car et Parrinello où la propagation classique des orbitales moléculaires est
eectuée, cette approche considère l'évolution classique de la matrice densité. Cette
méthode porte le nom de ADMP (de l'anglais Atom-centered Density Matrix Propa-

gation ) et est basée sur le lagrangien étendu suivant :

1 1X ~˙ I2
LADM P = µTr (WW) + MI R
2 2 I
h i (C.17)
− Tr [Λ (PP − P)] − Ee P, {R ~I}

où W représente la matrice des  vitesses  associée à la matrice densité P. Les

équations de mouvement découlant de ce lagrangien étendu (C.17) sont intégrées dans
le temps au moyen du propagateur Verlet aux vitesses (cf. Ÿ III.3 p. 22) et la structure
électronique peut être décrite dans une approche monocongurationnelle, i.e. Hartree-
Fock ou DFT. Les détails de cette approche ne seront pas discutés par la suite, toutefois
quelques points particuliers sont rappelés.
Une des particularités de ce schéma est que, à chaque pas de la propagation, la
matrice densité est  puriée  grâce à la transformation de Mac Weeny [94] : P
e =

3P2 − 2P3 . Cette transformation permet, en particulier, d'assurer l'idempotence de la

matrice densité. Un autre point particulier est l'utilisation de masses ctives diérentes
pour les diérents éléments de la matrice densité [11]. La masse ctive µ devient alors
une matrice :
p
1/2
µii = si Fii ≥ −2 a.u.
Me (C.18)
p h p i
1/2
µii = Me 2 |Fii + 2| + 1 si Fii < −2 a.u. (C.19)

µij = 0 si i 6= j (C.20)
1/2

où Me désigne ainsi la masse ctive associée aux orbitales de valence et Fii représente
un élément de la matrice de Fock. Ce schéma permet d'aecter des masses plus impor-
III Propagation Classique de la Matrice Densité 63

tantes aux orbitales de c÷ur, ce qui se justie assez naturellement puisque ces orbitales
possèdent un fort caractère atomique. Ces orbitales sont donc, a priori, plus inertes,
au cours d'une dynamique, que les orbitales de valence.

Cette méthode a été éprouvée et comparée à des trajectoires de type  Born-

Oppenheimer  [95]. Ces travaux ont montré l'ecacité de ce schéma, notamment
en terme de conservation de l'énergie totale du système ; de plus les redistributions
énergétiques rotationnelles et vibrationnelles aux cours des trajectoires ont été analy-
sées, montrant ainsi une reproduction satisfaisante des trajectoires Born-Oppenheimer
de référence. Néanmoins, un point n'est pas abordé dans les diérents articles discutant
de cette méthode. La pertinence de la méthode de Car et Parrinello est habituellement
justiée grâce à la séparabilité énergétique entre les noyaux et les orbitales moléculaires :
l'énergie cinétique ctive associée à la fonction d'onde est, en général, inférieure  de
plusieurs ordres de grandeur  à l'énergie cinétique nucléaire. L'évolution temporelle
du système est alors vue comme l'évolution adiabatique de deux sous-systèmes où les
échanges d'énergie sont minimes. Les exemples illustrant la méthode ADMP montrent
un comportement tout autre [11] : l'énergie cinétique ctive associée à la matrice den-
sité est relativement  importante  (de l'ordre de 10−3 à 10−1 u.a.) mais surtout elle
croît rapidement au cours du temps, mettant ainsi en évidence des échanges d'énergie
non négligeables entre les noyaux et la fonction d'onde. À la vue de ce comportement,
on peut donc s'interroger quant à la validité de cette approche sur des temps de simu-
lation longs. Il s'agit toutefois d'un exemple particulier et il serait dangereux de tirer
des conclusions trop hâtives ; des tests sur des temps longs et des systèmes diérents
permettraient certainement d'éclaircir ce point.
64 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

IV Une Variante : entre Born-Oppenheimer et Car-

Parrinello

1
Inter utrumque tene, medio tutissimus ibis

La méthode de dynamique moléculaire développée durant cette thèse est fortement

inspirée de l'approche proposée par Car et Parrinello. Au contraire de la méthode
développée par H.B. Schlegel qui envisage la propagation classique de la matrice densité,
l'approche proposée ici considère la propagation classique de la fonction d'onde via les
orbitales moléculaires. Cette approche n'est toutefois pas une dynamique moléculaire
de type Car-Parrinello au sens strict du terme puisque quelques cycles SCF (entre
deux et trois) sont eectués après la propagation classique des orbitales moléculaires.
Cette approche se situe donc  entre  les méthodes de type Car-Parrinello et Born-
Oppenheimer. L'intégration des équations du mouvement a été d'abord réalisée grâce
au propagateur de type Verlet original, puisque cet algorithme est :

 simple et ecace, car seules les forces sont nécessaires, pas les dérivées d'ordre
supérieur de l'énergie
 relativement précis, d'ordre trois
 explicitement réversible par rapport au temps
 symplectique, c'est-à-dire qu'il conserve le volume de l'espace des phases.

Toutes ces qualités contribuent à la stabilité sur des temps de simulation longs, notam-
ment en terme de conservation de l'énergie, de l'algorithme de Verlet. Cette implémen-
tation ainsi que les premiers tests concluants quant à l'ecacité de ce schéma ont fait
l'objet d'une publication, insérée par la suite.

1 Reste entre les deux, au milieu tu chemineras en sûreté. Vers d'Ovide (Métamorphoses, II, 137).
C'est le conseil par lequel le Soleil, conant à regret son char à Phaéton, son ls, termine la recom-
mandation qu'il vient de lui faire de n'approcher trop ni du ciel ni de la terre.
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello 65

IV.1 Article  C. Raynaud et al., Phys. Chem. Chem. Phys.,

6 (2004) 4226

Cet article traite de l'approche de dynamique moléculaire développée au cours de

cette thèse. Basée sur une approche de type Car-Parrinello, les orbitales moléculaires
développées sur une base de fonctions gaussiennes sont classiquement propagées, à l'aide
de l'algorithme original de Verlet, puis quelques cycles SCF (entre deux et trois) sont
eectués. Cet article traite plus précisément du comportement du schéma développé
notamment en s'intéressant à la conservation de l'énergie totale réelle au cours des tra-
jectoires simulées. Ces premiers tests ont été réalisés sur trois systèmes diérents : une
molécule d'eau isolée, l'acide monothiooxalique et enn un agrégat de vingt molécules
d'eau. Le comportement énergétique de ces diérentes trajectoires est comparé à des
simulations de référence de type Born-Oppenheimer, où la minimisation de l'énergie
potentielle est eectuée à chaque pas de la dynamique moléculaire. Ces premiers tests
nous permettent de conclure quant à un comportement satisfaisant de cette approche
mixte, que l'on pourrait situer  entre  Born-Oppenheimer et Car-Parrinello.
66 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello 67
68 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello 69
70 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello 71
72 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello 73

IV.2 Pertinence de la méthode

L'une des premières questions que l'on pourrait formuler vis à vis de la méthode
développée durant cette thèse est à propos de l'utilité de ces quelques cycles SCF eec-
tués après la propagation classique des orbitales moléculaires. En eet, dans l'approche
de type Born-Oppenheimer, à chaque pas de la propagation, le vecteur d'essai naturel-
lement choisi pour la fonction d'onde est la fonction d'onde optimisée au pas de temps
précédent. On peut alors se demander comment se comporterait une approche de type
Born-Oppenheimer  hybride , où le nombre de cycles SCF serait restreint à deux ou
trois.
0,04

0,03

0,02
∆Et (u.a.)

0,01

0 100 200
Temps (fs)

Fig. C.1  Déviation de l'énergie totale réelle par rapport à l'énergie totale initiale pour
trois trajectoires de conditions initiales identiques. L'approche Born-Oppenheimer est
en trait plein, celle développée durant cette thèse est en tirets et celle de type Born-
Oppenheimer  restreint  à quelques cycles SCF est en pointillés.

La gure C.1 présente la déviation de l'énergie totale réelle par rapport à l'éner-
gie totale initiale pour trois trajectoires, de conditions initiales identiques, associées
à trois approches diérentes : Born-Oppenheimer, Born-Oppenheimer  restreint  à
trois cycles SCF et enn la méthode développée durant cette thèse. Le système choisi
est l'un des exemples de l'article précédent : l'acide monothiooxalique. La structure
74 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

électronique de ce composé est décrite dans le cadre de la DFT au moyen de la fonc-

tionnelle B3LYP ; la base utilisée pour les atomes d'hydrogène est de qualité triple-ζ ,
6-311++G(d,p) [96], augmentée d'une fonction de polarisation et d'une fonction dif-
fuse ; des pseudo-potentiels de Stuttgart ont été utilisés pour les hétéroatomes avec
leurs bases associées [97]. Les équations du mouvement sont intégrées avec un pas de
temps de 0.1 fs. Dans le cas d'une dynamique de type Born-Oppenheimer restreint,
très rapidement, l'énergie totale n'est pas conservée ; ce type d'approche n'est donc
pas valide. Au contraire, l'approche développée présente une conservation de l'énergie
totale honorable au regard de la trajectoire Born-Oppenheimer de référence. Bien que
cette étape ait initialement été introduite par simplicité, elle permet de se rapprocher
en qualité de la surface exacte Born-Oppenheimer et d'assurer ainsi une bonne conser-
vation de l'énergie totale réelle sur des temps de simulation longs, tout en économisant
du temps de calcul.

IV.3 Application en Verlet aux vitesses

La méthode de dynamique moléculaire développée durant cette thèse a ensuite été

intégrée à l'aide du propagateur Verlet aux vitesses (cf. partie A Ÿ III.3 p. 22). En eet,
bien que cette variante de l'algorithme de Verlet génère des trajectoires identiques dans
l'espace des congurations par rapport au propagateur original, il évite de calculer par
diérence nie les vitesses des particules. Cet algorithme permet donc une meilleure
évaluation de la contribution cinétique à l'énergie totale. De plus, cet algorithme per-
met une implémentation plus aisée d'outils de la dynamique moléculaire, tels que les
thermostats ou bien les contraintes géométriques.
Cette manipulation explicite des vitesses lors de la propagation conduit à des ré-
sultats honorables de la méthode. Nous illustrons ici cette approche par des tests sur
la molécule d'hexatriène, représentée à la gure C.2. La structure électronique de cette
molécule a été décrite dans l'approximation Hartree-Fock, dans une base 6-31G(d,p) de
qualité double-ζ , augmentée d'une fonction de polarisation, pour les hydrogènes [96] ;
IV Une Variante : entre Born-Oppenheimer et Car-Parrinello 75

Fig. C.2  Molécule d'hexatriène

des pseudo-potentiels de Stuttgart [97] associés à leur base optimisée, augmentée d'une
fonction de polarisation d, ont été utilisés pour les atomes de carbone. L'utilisation de
pseudo-potentiels permet ici d'éviter le traitement explicite des orbitales 1s de c÷ur
des carbones et de ne propager classiquement que les orbitales de valence. Chaque
trajectoire a été conduite à l'aide d'un pas de temps de 0.2 f s sur un temps total
de simulation de 5 ps, avec une énergie cinétique initiale de 24 × 10−3 u.a., soit une
 température équivalente  de 200 K . Le comportement du schéma de propagation
est satisfaisant puisque l'énergie totale est ici conservée avec une précision honorable :
l'énergie totale moyenne dévie de 3.1 × 10−5 u.a. au bout de 5 ps de simulation et la
déviation maximale rapportée à l'énergie totale initiale est de 5.4 × 10−5 u.a.. Ces ré-
sultats sont donc honorables au regard du temps relativement long des trajectoires. Le
gain de temps observé pour ce système moléculaire est d'un facteur trois par rapport à
une dynamique moléculaire de type Born-Oppenheimer. Bien que ce facteur de temps
ne soit pas très important, il est néanmoins relativement satisfaisant au regard du coût
énorme des approches de dynamique moléculaire ab initio, compte tenu des capacités
de calcul des ordinateurs actuels.

Il est important de noter que nous avons fait le choix de minimiser l'énergie électro-
nique par rapport à la fonction d'onde, si la propagation classique des orbitales n'est
pas ecace à un instant donné. Cette procédure permet en quelque sorte de réinitia-
liser la propagation classique de la fonction d'onde au cours d'une trajectoire. Il est
néanmoins dicile de prédire de façon générale la fréquence de ces re-convergences. En
eet, pour les systèmes les plus  récalcitrants  que nous avons traités au cours de
cette thèse, il se produit deux à cinq procédures de re-convergence sur quelques fem-
tosecondes de simulation, ce qui est largement négligeable du point de vue du temps
76 DYNAMIQUE MOLÉCULAIRE CAR-PARRINELLO

de calcul. Tous les résultats présentés auparavant sont issus de trajectoires où aucune
procédure de re-convergence n'a été eectuée.
Notons enn que ce schéma de dynamique moléculaire ne présente pas de sensibilité
particulière vis à vis de la masse ctive associée aux orbitales moléculaires, au contraire
de l'approche Car-Parrinello stricto sensu en ondes planes. On peut néanmoins dégager
une gamme pour ce choix, entre 100 et 800 amu.bohr2 , dans laquelle la propagation
classique des orbitales moléculaires est ecace. En dehors de cette gamme, l'algorithme
montre une trop grande fréquence de minimisation de l'énergie potentielle.

V Conclusion
La méthode de dynamique moléculaire développée au cours de cette thèse semble
donc ecace et adaptée à l'étude de systèmes moléculaires. En particulier le comporte-
ment de ce schéma s'avère globalement satisfaisant au regard de la conservation, sur des
temps longs, de l'énergie totale réelle. Même si le gain en temps de calcul par rapport
à une approche de type Born-Oppenheimer n'est pas de plusieurs ordres de grandeur,
il reste néanmoins tout à fait acceptable, compte tenu des coûts rapidement prohibitifs
de ces approches. L'ensemble des outils présentés précédemment permet donc de si-
muler des trajectoires de dynamique moléculaire dans l'ensemble microcanonique. Or,
les expériences auxquelles les résultats théoriques sont confrontés, sont le plus sou-
vent réalisées dans des conditions diérentes. En particulier les eets de température
ne peuvent pas être directement pris en compte dans l'ensemble microcanonique. La
partie suivante explique comment ces eets peuvent être inclus lors de simulations de
dynamique moléculaire. Plus particulièrement, les thermostats qui ont été utilisés et
implémentés avec les précédents outils sont présentés.
D Au Delà du Microcanonique
Le but de cette partie est de préciser comment, en dynamique moléculaire, les eets de
température peuvent être pris en compte de manière explicite. Après un bref rappel
à propos des principes de la mécanique statistique non hamiltonienne, un des outils
permettant de simuler l'ensemble canonique est précisé. Enn, des exemples illus-
trant la mise en application de ces outils avec l'approche de dynamique moléculaire
développée durant cette thèse sont présentés.

L a réalisation de simulations de dynamique moléculaire au moyen des outils pré-

sentés dans les précédentes parties permet de réaliser des trajectoires dans l'en-
semble microcanonique, i.e. les eets de température ne sont pas pris en compte de
manière explicite. La dynamique moléculaire non hamiltonienne est typiquement uti-
lisée pour générer les ensembles canonique, isotherme-isobare ou bien pour décrire des
systèmes sujets à des contraintes holonomes. L'idée de générer par dynamique molé-
culaire de tels ensembles débuta grâce aux travaux d'Andersen [98] qui proposa un
schéma de dynamique à pression constante. Ces travaux furent ensuite étendus pour
générer des distributions canoniques et grand-canoniques [99, 100]. Ces formulations
originales étaient basées sur une extension de l'espace des phases de systèmes hamil-
toniens mais souraient de certains traits indésirables liés à la dénition du temps.
Toutefois, W. G. Hoover en 1984 apporta une correction à l'aide d'une reformulation
non hamiltonienne [101]. Enn, ce n'est que récemment que des travaux proposèrent
de justier la validité statistique de la dynamique moléculaire non hamiltonienne [102],
dont quelques principes sont brièvement introduits par la suite.
78 AU DELÀ DU MICROCANONIQUE

I Principe de la Mécanique Statistique non Hamilto-

nienne

Considérons un système dynamique dont les équations du mouvement pour un

vecteur x de l'espace des phases s'écrivent selon :

ẋ = ζ(x, t) (D.1)

où ζ(x, t) représente la force généralisée, i.e. cette force n'est plus a priori conservative
et peut avoir une dépendance explicite en temps. La solution de cette équation dépend
des conditions initiales x0 = (x10 , ..., xn0 ). Cette solution conduit à un jeu de n vecteurs
qui sont chacun fonctions du temps et des conditions initiales :
¡ ¢
xit = xit t; x10 , ..., xn0 (D.2)

Le système est donc non hamiltonien, et on peut dénir la compressibilité de l'espace

des phases selon :
κ(x, t) = ∇ẋ = ∇ζ(x, t) (D.3)

où ∇ représente le gradient sur l'espace des phases. Cette compressibilité, qui s'annule
pour les systèmes hamiltoniens (cf. partie A, Ÿ III.2 p. 21), est dans ce cas non nulle
et la métrique de l'espace des phases dx n'est plus invariante. Pour préciser cela, on
peut considérer le jacobien de la transformation de coordonnées de la relation (D.2),
i.e. d'un vecteur initial de l'espace des phases x0 au vecteur xt représentant l'état au
temps (t) :
∂ (x1t , ..., xnt )
J(xt ; x0 ) = (D.4)
∂ (xt0 , ..., xn0 )
On peut montrer [102] que ce jacobien satisfait une relation de la forme suivante :

d
J(xt ; x0 ) = J(xt ; x0 )∇ẋ = J(xt ; x0 )κ(x, t) (D.5)
dt

avec pour condition initiale J(0) ≡ J(x0 , x0 ) = 1. Cette relation implique notamment
que le jacobien est égal à 1 à n'importe quel instant (t) si la compressibilité est nulle,
I Principe de la Mécanique Statistique non Hamiltonienne 79

comme dans le cas de systèmes hamiltoniens. Ce jacobien peut être utilisé pour ex-
primer un élément de volume élémentaire (i.e. la métrique) de l'espace des phases au
temps (t) en fonction de l'élément de volume élémentaire au temps initial :

dxt = J(xt ; x0 )dx0 (D.6)

On voit ainsi que cette métrique n'est plus invariante pour les systèmes non hamilto-
niens ; en eet, si J 6= 1 alors dxt 6= dx0 . La théorie statistique complète développée
pour les systèmes non hamiltoniens [102] montre qu'il existe une métrique invariante,
dénie selon :
p
dµ = g(x, t)dx (D.7)
p
où le facteur  de mesure  g(x, t) est donné par la relation (D.8) :

p
g(x, t) = e−w(x,t) (D.8)

et la fonction w(x, t) est reliée à la compressibilité κ(x, t) par :

dw
=κ (D.9)
dt

Cette métrique invariante permet de dénir une équation de Liouville généralisée de la

forme : ¡ √ ¢
∂ ρ g √
+ ∇ (ρ g ẋ) = 0 (D.10)
∂t
où ρ désigne la fonction de distribution, i.e. la densité d'états. Supposons ensuite que le
jeu d'équations du mouvement (D.1) possède nc invariants Kλ , i.e. le système possède
nc lois de conservation de la forme :

dKλ
=0 (D.11)
dt

Ainsi, une trajectoire générée par les équations du mouvement (D.1) n'explore pas
la totalité de l'espace des phases mais le sous-espace déterminé par l'intersection des
hypersurfaces dénies par {Kλ (x) = K
e λ } où {K
e λ } est un jeu de constantes. Pour ce
80 AU DELÀ DU MICROCANONIQUE

système dynamique, la densité d'états microcanonique peut alors être dénie à partir
d'un produit de distributions :
nc
Y
ρ(x) = eλ)
δ(Kλ (x) − K (D.12)
λ=1

Il faut noter qu'il est important de déterminer toutes les lois de conservation satisfaites
par les équations du mouvement d'un système non hamiltonien. En eet, une densité
d'états peut être construite à partir d'un produit de distributions d'un sous-ensemble
d'invariants : ′
nc
Y
ρ(x)red = eλ)
δ(Kλ (x) − K (D.13)
λ=1

où n′c < nc . Une telle densité d'états réduite ρ(x)red satisfait l'équation de Liouville
généralisée (D.10), elle ne décrit cependant pas correctement la distribution de l'espace
des phases d'un système à nc lois de conservation.
À partir de la densité d'états dénie par la relation (D.12), la fonction de partition
microcanonique de ce système non hamiltonien peut être dénie selon :
Z p nc
Y
Ω(N, V, K e nc ) =
e 1 , ..., K dx g(x) eλ)
δ(Kλ (x) − K (D.14)
λ=1

Cette fonction de partition pour un système non hamiltonien doit permettre de retrou-
ver la distribution canonique si l'on souhaite se placer dans cet ensemble. Le paragraphe
suivant précise les outils utilisés lors de cette thèse pour générer l'ensemble canonique :
les chaînes de thermostats de Nosé-Hoover.
II Chaînes de Thermostats de Nosé-Hoover 81

II Chaînes de Thermostats de Nosé-Hoover

Par la suite, nous précisons comment les principes de base de la mécanique statis-
tique non hamiltonienne présentés précédemment peuvent être utilisés pour construire
des équations du mouvement pour la dynamique moléculaire générant l'ensemble ca-
nonique. La méthode des thermostats de Nosé-Hoover est un schéma de dynamique
moléculaire non hamiltonienne permettant de générer un tel ensemble. Cette méthode
est basée sur l'approche initiale de Nosé [99, 100], reformulée par Hoover [101]. Néan-
moins, dans le cas de systèmes possédant un faible nombre de degrés de liberté, l'ajout
d'une variable de thermostat unique ne sut pas à créer assez de  chaos  pour que
la dynamique du système étendu échantillonne correctement l'espace des phases. Ce
problème pathologique frappe en particulier l'oscillateur harmonique [103]. Parmi les
solutions proposées pour résoudre ce problème et pour rendre la dynamique moléculaire
 davantage  ergodique, la méthode des chaînes de thermostats de Nosé-Hoover [103]
se distingue par son élégance et sa simplicité.

II.1 Dénition
La méthode des chaînes de Nosé-Hoover étend l'espace des phases habituel à un jeu
de M variables de thermostat ξ1 , ..., ξM et à leurs impulsions associées. Ces variables
de thermostat jouent le rôle d'un réservoir de chaleur couplé au système. Les équations
du mouvement prennent la forme suivante :

~
~˙ I = PI
R (D.15)
MI

˙
P~I = F~I − ξ˙1 P~I (D.16)

à !
1 X P~ 2
ξ¨1 = I
− gkB T − ξ˙1 ξ˙2 (D.17)
Q1 I
M I
82 AU DELÀ DU MICROCANONIQUE

1 ³ ´
ξ¨ν = Qν−1 ξ˙ν−1
2
− kB T − ξ˙ν ξ˙ν+1 ν = 2, ..., M − 1 (D.18)
Qν

1 ³ ´
¨
ξM = ˙2
QM −1 ξM −1 − kB T (D.19)
QM
La première variable de la chaîne de thermostats ξ˙1 est directement couplée au sys-
tème. Elle joue le rôle d'un coecient de friction qui module l'impulsion des particules
selon que la température  instantanée  est inférieure ou supérieure à la température
imposée T . En eet si la température instantanée est supérieure à la température im-
posée, alors, d'après la relation (D.17),  l'accélération  de ce premier thermostat est
positive et sa  vitesse  augmentera, ce qui aura pour eet, d'après la relation (D.16),
de diminuer la vitesse des particules et donc de rapprocher la température instantanée
de la température imposée. De plus, d'après la relation (D.17), ce premier thermostat
est couplé au second, le second est lui même couplé au troisième et ainsi de suite, selon
la relation (D.18). Enn, d'après l'équation (D.19) le dernier thermostat n'est couplé
qu'avec la variable précédente. Cette chaîne de thermostats couplés permet de prévenir
les uctuations incontrôlées du premier thermostat et de rendre  plus ergodique  la
dynamique moléculaire. Ces variables de thermostat sont aublées de masses Qν , dont
le choix est crucial pour l'échantillonnage de l'ensemble canonique. Si de très grandes
masses sont choisies, les thermostats seront trop inertes et la distribution générée sera
proche d'une distribution microcanonique. Au contraire, si de trop petites masses sont
choisies, les uctuations des impulsions du système sont prohibées. Le choix judicieux
de ces masses permettant de générer une distribution canonique est le suivant :

Q1 = gkB T τ 2 (D.20)

Qν = kB T τ 2 ν = 2, . . . , M (D.21)

où g désigne le nombre de degrés de liberté du système moléculaire, i.e. g = 3N pour

un système de N particules à trois dimensions. La variable τ représente un temps
II Chaînes de Thermostats de Nosé-Hoover 83

caractéristique du système moléculaire. Ce choix permet aux thermostats d'être en

quelque sorte en  résonance  avec le système moléculaire.

II.2 Distribution canonique

L'utilisation de chaînes de thermostats de Nosé-Hoover introduit de nouvelles va-

riables, non physiques, au sein des équations du mouvement. Le principal invariant du
système est l'énergie totale, qui tient compte de termes non physiques associés aux
variables des thermostats et qui a pour expression :
M M
1X ~ I }] + 1
~˙ 2 + Ep [{ψi }, {R
X X
E = MI R I Q ξ˙
ν ν
2
+ gk B T ξ 1 + kB T ξν (D.22)
2 I 2 ν=1 ν=2

Cette énergie totale n'a donc pas de sens physique direct, elle doit néanmoins être
conservée lors des simulations pour justier la validité de la trajectoire générée. La
compressibilité de l'espace des phases étendu, dénie par la relation (D.3), peut ici
s'exprimer :
N h
" #
X i X M
∂ ˙
ξ ∂ ¨
ξ
κ= ~˙ I + ∇ ~ P~˙ I +
∇R~ I R
ν
+
ν
PI ˙
I=1 µ=1
∂ξ ν ∂ ξν
M
(D.23)
X
= −3N ξ˙1 − ξ˙ν
ν=2

À partir de la relation (D.9) reliant le facteur de mesure et la compressibilité de l'espace

des phases, ce facteur s'exprime dans le cas des chaînes de thermostats de Nosé-Hoover
selon :
M
X
w = −3N ξ1 − ξν (D.24)
ν=2

On peut ainsi dénir la métrique invariante de l'espace des phases étendu à ce jeu de
variables non physiques :
à M
!
X Y Y
dµ = exp 3N ξ1 + ξν dP~I dR
~I dξν dpξν (D.25)
ν=2 I ν

De plus, si la dynamique du système n'est pas gouvernée par des forces externes, i.e.
P
lorsque I F~I = ~0, il existe trois (pour un système moléculaire à trois dimensions)
84 AU DELÀ DU MICROCANONIQUE

lois de conservation en plus de la conservation de l'énergie totale, dénie par la rela-

tion (D.22). Ces lois de conservation adoptent la forme suivante :
N
X
~ = eξ1
K P~I = eξ1 P~ (D.26)
I=1

où P~ désigne l'impulsion du centre de masse du système moléculaire. Ainsi, à partir

de la métrique dénie par la relation (D.25), des lois de conservation dénies par les
relations (D.22) et (D.26), on peut dénir la fonction de partition microcanonique du
système global grâce à l'équation (D.14). Dans le cas où la longueur de la chaîne de
thermostats est égale à deux, cette fonction de partition s'exprime selon :
Z Y
~ =
Ω(N, V, E , K) dP~I dR
~ I dξ1 dξ2 dpξ1 dpξ2 exp (3N ξ1 + ξ2 )
I
µ ³ ´ ¶
p2ξ1 p2ξ2 (D.27)
~
× δ H R, P + ~ + + 3N kB T ξ1 + kB T ξ2 − E
2Q1 2Q2
¡ ¢ ¡ ¢ ¡ ¢
× δ eξ1 Px − Kx δ eξ1 Py − Ky δ eξ1 Pz − Kz

An de démontrer que cette distribution est équivalente à une distribution canonique,
l'intégrale sur les variables non physiques associées aux thermostats doit être réalisée.
Après un développement quelque peu fastidieux mais pas dicile, on peut simplier
cette expression :
Z Y ³ ³ ´´
~ ∝
Ω(N, V, E , K) ~ I dP~I exp −βH R,
dR ~ P~ (D.28)
I

où β = 1/kB T représente la température inverse. On démontre ainsi que la distribution

générée par la dynamique d'un système moléculaire couplé à une chaîne de thermostats
de Nosé-Hoover conduit à une fonction de partition qui est proportionnelle à la fonction
de partition de l'ensemble canonique.

II.3 L'intégration des équations du mouvement

Les variables non physiques associées aux thermostats ont une échelle de temps
caractéristique en général petite par rapport à celle du système moléculaire. Par consé-
quent, un petit pas de temps d'intégration doit être utilisé pour correctement décrire la
II Chaînes de Thermostats de Nosé-Hoover 85

dynamique de ces variables non physiques et donc de bien décrire les échanges d'énergie
entre la chaîne de thermostats et le système moléculaire. Une procédure inecace serait
de simplement réduire le pas de temps d'intégration. En eet, les forces agissant sur le
système moléculaire seraient calculées  trop souvent  alors qu'elles ne seraient pas
modiées de manière signicative. Une approche plus ecace est de considérer le sys-
tème global comme la réunion de deux sous-systèmes aux échelles de temps diérentes
et d'intégrer les équations du mouvement de chaque sous-système avec le pas de temps
adapté à chacun. Une telle procédure permet alors de ne pas calculer les forces agissant
sur le système moléculaire à chaque pas d'intégration de la chaîne de thermostats. Les
approches pour l'intégration des équations du mouvement des systèmes aux échelles
de temps multiples peuvent être développées à partir de la factorisation de Trotter de
l'opérateur de Liouville. Cette méthode est brièvement décrite en annexe (cf. annexe
I p. 201). L'intégration des équations du mouvement d'un système moléculaire couplé
à une chaîne de thermostats de Nosé-Hoover peut être eectuée en utilisant le même
principe. L'opérateur de propagation de Liouville peut tout d'abord être décomposé en
somme de trois opérateurs selon :

iL = iL1 + iL2 + iLN HC (D.29)

Les opérateurs iL1 et iL2 sont dénis de la même manière que précédemment (cf.
partie A Ÿ III.5 p. 26) dont l'expression est rappelée :
XN
iL1 = ~˙ I ∂
R (D.30)
~I
∂R
I=1

N
X ∂
iL2 = F~I (D.31)
I=1 ∂ P~I
Les opérateurs iL1 et iL2 sont donc les opérateurs de  translation  associés respecti-
vement aux positions et aux impulsions nucléaires. L'opérateur iLN HC est quant à lui
associé à la chaîne de thermostats de Nosé-Hoover. Son expression est la suivante :
X XM M
X −1 ³ ´ ∂
iLN HC = − ˙ξ1 P~I ∂ + ˙ξν ∂ + ξ¨ν − ξ˙ν ξ˙ν+1 + ξ¨M
∂
(D.32)
I ∂ P~I ν=1 ∂ξν ν=1 ∂ ξ˙ν ∂ ξ˙M
86 AU DELÀ DU MICROCANONIQUE

Un développement analogue aux relations (A.79) et (A.80) conduit à la dénition du

propagateur discret en temps suivant :

µ ¶ µ ¶
δt δt
G(δt) = exp iLN HC × exp iL1
2 2
× exp (iL2 δt) (D.33)
µ ¶ µ ¶
δt δt ¡ ¢
× exp iL1 × exp iLN HC + O (δt)3
2 2

L'opérateur associé aux variables des thermostats est alors développé de telle façon
qu'un pas d'intégration plus petit soit utilisé pour intégrer les équations du mouvement
de la chaîne de thermostats. Cet opérateur prend la forme suivante :

µ ¶ n
Y µ ¶
δt δt
exp iLN HC = exp iLN HC (D.34)
2 j=1
2n

Cette approche revient donc à propager une fois avec un pas de temps (δt) le système
moléculaire et 2n fois avec un pas de temps (δt)/2n les variables de la chaîne de
thermostats. De plus, d'après les relations (D.17), (D.18) et (D.19), les forces agissant
sur les variables des thermostats sont simples à évaluer, cette approche ne pénalise
donc pas l'eort calculatoire.

Grâce à la relation (D.32), on voit nettement que les variables de thermostat sont
couplées entre elles. L'opérateur d'évolution iLN HC ne peut donc pas être utilisé en
l'état, néanmoins on peut utiliser encore une fois le théorème de Trotter pour simplier
son expression. Ainsi, l'opérateur iLN HC déni par la relation (D.32) peut être exprimé
II Chaînes de Thermostats de Nosé-Hoover 87

sous la forme :
µ ¶ µ ¶ µ ¶
δt δt ¨ ∂ δt ˙ ˙ ∂
exp iLN HC = exp ξM × exp − ξM ξM −1
2n 4n ∂ ξ˙M 8n ∂ ξ˙M −1
µ ¶ µ ¶
δt ¨ ∂ δt ˙ ˙ ∂
× exp ξM −1 × exp − ξM ξM −1
4n ∂ ξ˙M −1 8n ∂ ξ˙M −1
× ...
à ! à M
!
δt X ˙ ∂ δt X ˙ ∂
× exp − ξ 1 PI × exp ξν (D.35)
2n I ∂PI 2n ν=1 ∂ξν

× ...
µ ¶ µ ¶
δt ˙ ˙ ∂ δt ¨ ∂
× exp − ξM ξM −1 × exp ξM −1
8n ∂ ξ˙M −1 4n ∂ ξ˙M −1
µ ¶ µ ¶
δt ˙ ˙ ∂ δt ¨ ∂
× exp − ξM ξM −1 × exp ξM
8n ∂ ξ˙M −1 4n ∂ ξ˙M

L'approche décrite ici semble bien indigeste et compliquée. Elle est cependant simple à
implémenter et ne demande pas d'eort calculatoire supplémentaire : le facteur limitant
est toujours l'évaluation des forces agissant sur le système moléculaire. Elle permet de
plus une bonne précision quant à la dynamique intrinsèque des thermostats et donc de
bien décrire les échanges d'énergie entre le système moléculaire et la chaîne.
La précision quant à l'intégration des équations du mouvement des variables de
la chaîne de thermostats peut encore être améliorée grâce au schéma d'intégration de
Yoshida-Suzuki [104, 105]. Cette approche, brièvement décrite en annexe (cf. annexe II
p. 205), permet de réduire l'erreur associée au développement de l'opérateur d'évolution
sans eectuer de développement à des ordres supérieurs. L'expression (D.34) associée
à l'opérateur d'évolution de la chaîne de thermostats adopte la forme :
µ ¶ n
"m µ ¶#
δt Y Y ω δt
exp iLN HC = exp iLN HC k (D.36)
2 j=1 k=1
2n

où ωk représente les facteurs de Yoshida-Suzuki  simpliés  (cf. annexe II p. 205).

Les chaînes de thermostats de Nosé-Hoover peuvent nalement être incluses dans
un schéma de dynamique moléculaire au moyen de l'algorithme d'intégration présenté
ci-dessus. Ces thermostats permettent alors la simulation de dynamique moléculaire
88 AU DELÀ DU MICROCANONIQUE

dans l'ensemble canonique. Ces simulations sont utiles pour l'estimation de certaines
propriétés. En eet, ces mêmes propriétés sont le plus souvent mesurées expérimen-
talement à une température donnée. Ces propriétés peuvent être en particulier des
observables spectroscopiques. Par la suite nous présentons des exemples d'utilisation
du schéma de dynamique moléculaire développé au cours de cette thèse, combiné avec
les chaînes de thermostats de Nosé-Hoover, pour l'estimation de caractéristiques spec-
troscopiques.

III Calculs de Propriétés

Les paramètres uniquement accessibles à la dynamique moléculaire sont ceux qui
s'expriment explicitement en fonction du temps. L'utilisation de chaînes de thermostats
de Nosé-Hoover permet d'identier, à l'aide de l'hypothèse ergodique, les moyennes
temporelles et les moyennes canoniques. En pratique toutefois, pour  vérier  l'hy-
pothèse ergodique, les temps simulés se doivent d'être les plus longs possibles. De plus,
 l'équilibre  doit être atteint, ce qui nécessite au préalable des temps de thermalisa-
tion longs. Ainsi, sous réserve d'une thermalisation ecace et d'un temps simulé assez
long, la moyenne canonique d'une observable s'identie à sa moyenne temporelle sur
les trajectoires générées.

En chimie expérimentale, il existe un grand nombre de méthodes permettant la

caractérisation des molécules. Dans un souci de compréhension des phénomènes expé-
rimentaux mais aussi dans un but prédictif, plusieurs méthodes de la chimie quantique
ont été développées. Parmi toutes ces techniques de caractérisation, certaines sont uti-
lisées de façon quasi-routinière, c'est entre autres le cas de la spectroscopie infrarouge,
de la résonance magnétique nucléaire (RMN) ou de la spectroscopie ultraviolet-visible
(UV).
III Calculs de Propriétés 89

III.1 Simulation de spectres infrarouge

Une des informations, que l'on peut extraire de la dynamique moléculaire est le
spectre de puissance P(ω). Celui-ci permet de caractériser le nombre d'oscillateurs de
fréquence donnée, à l'énergie totale de la simulation si l'ensemble simulé est microcano-
nique, ou bien à une température donnée si cet ensemble est canonique. Ce spectre de
puissance est obtenu via la fonction d'autocorrélation des vitesses, grâce au théorème
de Wiener-Khintchine [42] :
DP E
Z ~ ~
I VI (t) · VI (0)
∞
P(ω) = 2 C(t) cos ωtdt avec C(t) = DP E (D.37)
~ ~
I VI (0) · VI (0)
0

Ce spectre de puissance permet de déterminer si un mode propre de vibration ω est

peuplé. L'activité d'un mode de vibration en spectroscopie infrarouge peut quant à
elle être estimée classiquement par la fonction d'autocorrélation du moment dipolaire
µ [106, 107, 108] selon :
Z +∞
a(ω) ∝ ω 2
hµ(t) · µ(0)i exp(−iωt)dt (D.38)
−∞

En pratique il est toutefois avantageux numériquement [109] d'utiliser l'une des pro-
priétés des transformées de Fourier, permettant de reformuler la relation (D.38) :
Z +∞
a(ω) ∝ hµ̇(t) · µ̇(0)i exp(−iωt)dt (D.39)
−∞

Le schéma de dynamique moléculaire développé durant cette thèse, combiné à des

chaînes de thermostats de Nosé-Hoover, a permis d'étudier plusieurs systèmes molécu-
laires. Nous présentons par la suite un article, s'intéressant plus particulièrement au
spectre infrarouge de petites molécules.

III.2 Article  C. Raynaud et al. , Chem. Phys. Lett. , 414

(2005) 161
Cet article traite de l'estimation de caractéristiques infrarouge pour des petits agré-
gats de uorure d'hydrogène, (HF )n n = 2−7. Ces résultats, obtenus à l'aide du schéma
90 AU DELÀ DU MICROCANONIQUE

de dynamique moléculaire développé, sont comparés avec ceux obtenus par des simu-
lations de dynamique moléculaire de type Car-Parrinello en ondes planes. Cet article
met plus particulièrement en évidence la diculté de traiter correctement la structure
électronique de ce type de composés par des bases d'ondes planes. En eet, pour une
précision comparable, le schéma de dynamique moléculaire en base locale s'avère dans
ce cas plus performant que les approches en base délocalisée.
III Calculs de Propriétés 91
92 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 93
94 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 95
96 AU DELÀ DU MICROCANONIQUE

III.3 Simulation de spectres RMN

La spectroscopie RMN est une spectroscopie  lente  du point de vue des temps
caractéristiques vibrationnels moléculaires. En eet, le temps de relaxation d'un spin
nucléaire est bien plus grand que le temps caractéristique de vibration au sein d'une mo-
lécule. Par conséquent, le spin nucléaire, lors de sa relaxation, ne voit qu'une moyenne
des eets vibrationnels. Il existe plusieurs méthodes de chimie quantique pour déter-
miner le tenseur d'écran électronique [110, 111] et les constantes de couplage spin-
spin [112]. Néanmoins ces méthodes ne prennent pas en compte, a priori, le fait que
le temps de relaxation d'un spin nucléaire est plus grand que les temps vibrationnels
caractéristiques. En eet ce type de calcul est eectué sur une structure géométrique
gée. Nous présentons par la suite un article, à propos des propriétés RMN d'une
molécule où il peut se produire une pseudo-rotation de Berry [113, 114] : P F5 .

III.4 Article  C. Raynaud et al. , accepté à Chem. Phys. Chem.

Une manière pour prendre en compte les eets de température pour l'estimation
de tenseurs d'écran électronique est de moyenner plusieurs calculs réalisés sur un jeu
de structures géométriques diérentes. Cet article traite de l'estimation des déplace-
ments chimiques de la molécule de pentauorophosphore, dans laquelle un échange
de phosphore, par un mécanisme de pseudo-rotation de Berry [113, 114], peut se pro-
duire. L'ensemble des structures obtenues pour l'estimation des constantes d'écran
électronique ainsi que pour les constantes de couplages est issu de trajectoires de dy-
namique moléculaire réalisées à température nie. En particulier l'ensemble des tra-
jectoires conduites à une température de 1000 K permet de retrouver, d'un point de
vue théorique, l'équivalence des cinq atomes de uor expérimentalement observée en
spectroscopie RMN.
III Calculs de Propriétés 97
98 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 99
100 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 101
102 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 103
104 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 105
106 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 107
108 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 109
110 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 111
112 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 113
114 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 115
116 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 117
118 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 119
120 AU DELÀ DU MICROCANONIQUE

III.5 Simulation de spectres UV-visible

La prédiction par les méthodes de la chimie théorique des caractéristiques UV-
visible des molécules relève des méthodes permettant une description relativement
précise des états excités. Néanmoins, comme dans le cas précédent, ces méthodes
ne tiennent pas compte, a priori , des eets de température sur les spectres obtenus.
D'un autre côté, la spectroscopie UV-visible est une spectroscopie  rapide  vis-à-vis
des temps vibrationnels caractéristiques moléculaires. Or les mesures expérimentales
ne sont pas réalisées sur une seule molécule mais sur un ensemble macroscopique, le
spectre résultant est alors la somme des diérentes transitions électroniques observées
pour chaque molécule de l'échantillon. Nous présentons par la suite un article, à propos
de la simulation de spectres UV-visibles d'une molécule : la β−ionone.

III.6 Article  C. Raynaud et al. , soumis à J. Mol. Struct.

De la même façon que dans l'exemple précédent, la dynamique moléculaire ab initio

à température nie peut être utilisée pour générer un ensemble de structures géomé-
triques. Cet article traite de la simulation du spectre UV-visible d'une molécule or-
ganique : la β−ionone. Le spectre simulé est tout à fait satisfaisant par rapport au
spectre d'absorption électronique expérimental, en terme de position, forme et hauteur
relative des diérentes bandes. Il faut noter que la structure électronique de ce com-
posé a été décrite par la méthode hybride ONIOM [83] (cf. partie B Ÿ VI p. 47) : le
chromophore, i.e. le système π , a été traité à  haut  niveau par une méthode de DFT
et le  reste  de la molécule a été représenté à l'aide d'une méthode semi-empirique,
AM1 [115]. Seule la fonction d'onde associée à la partie modèle décrite par la DFT a
été propagée  classiquement  selon la méthode développée durant cette thèse ; lors
des calculs semi-empiriques, l'énergie a été minimisée en chaque point de la trajectoire.
Notons enn que pour réaliser ces trajectoires, les gradients de chaque calcul doivent
être recombinés an d'obtenir les gradients du système global [84].
III Calculs de Propriétés 121
122 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 123
124 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 125
126 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 127
128 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 129
130 AU DELÀ DU MICROCANONIQUE
III Calculs de Propriétés 131
132 AU DELÀ DU MICROCANONIQUE

IV Conclusion
En dynamique moléculaire ab initio, les eets de température peuvent être pris en
compte au moyen de chaînes de thermostats de Nosé-Hoover. Ces outils permettent de
générer une distribution cohérente avec l'ensemble canonique. Sous réserve d'une ther-
malisation ecace et de temps simulés assez longs, on peut identier la moyenne tem-
porelle d'une observable à sa moyenne dans l'ensemble canonique. Ces simulations dans
l'ensemble canonique permettent aussi, entre autres, de bien reproduire certaines ca-
ractéristiques spectroscopiques, telles que les données infrarouge, RMN ou UV-visible.
Néanmoins, certaines grandeurs ne peuvent pas simplement être évaluées par des si-
mulations de ce genre. C'est notamment le cas de l'entropie et des diérences d'énergie
libre pour une réaction chimique. En eet, pour estimer ces grandeurs, il faudrait
a priori simuler un grand nombre d'événements réactionnels. Ceci implique que les
temps simulés soient extrêmement grands, largement supérieurs à quelques dizaines,
voire quelques centaines de picosecondes. Un événement réactionnel, du fait de la hau-
teur relativement importante des barrières d'énergie à franchir, est alors un événement
rare pour la dynamique moléculaire ab initio et il est aujourd'hui impossible de géné-
rer un ensemble statistique permettant d'accéder à ces grandeurs par des dynamiques
 naturelles . Pour réaliser, dans des temps de simulation accessibles, un grand nombre
de ces événements réactionnels, la dynamique naturelle du système doit être modiée.
Dans le but d'estimer ces grandeurs thermodynamiques, plusieurs méthodes ont été
développées. La partie suivante traite de l'estimation des diérences d'énergie libre et
précise en particulier les approches que nous avons implémentées et utilisées au cours
de cette thèse pour l'évaluation des grandeurs thermodynamiques de réaction.
E Le Calcul des Énergies Libres
Le but de cette partie est de préciser comment, en dynamique moléculaire, le calcul des

diérences d'énergies libres peut être envisagé. L'approche de l'ensemble blue-moon

permettant l'évaluation de ces grandeurs est plus particulièrement présentée. Par la

suite, des exemples illustrant la mise en application de ces outils avec l'approche de

dynamique moléculaire développée durant cette thèse sont précisés. Enn, l'estimation

de l'énergie de  point zéro  au-delà de l'approximation harmonique est abordée.

L es simulations de dynamique moléculaire permettent d'accéder aux informations

statistiques d'un système moléculaire. En particulier elles permettent d'accéder
aux propriétés qui sont des fonctions  explicites  des coordonnées de l'espace des
phases. Nous avons présenté quelques exemples de telles propriétés au cours de la pré-
cédente partie : la connaissance de la moyenne d'ensemble ou de la moyenne temporelle
de ces propriétés mécaniques permet d'accéder à des caractéristiques spectroscopiques
ou bien encore à des données thermodynamiques importantes telles que les capacités
caloriques. Néanmoins, il n'est pas possible d'évaluer directement certaines grandeurs,
notamment celles dont la valeur dépend du volume total de l'espace des phases acces-
sible au système. Des exemples typiques de ce genre de quantité sont l'entropie d'un
système ainsi que l'énergie libre F . Par la suite, le terme d'énergie libre est réservé à
l'énergie libre d'Helmholtz. L'énergie libre d'un système moléculaire dépend ainsi de
la totalité du volume de l'espace des phases accessible au système, ceci se retrouve
aisément puisque cette dernière est reliée à la fonction de partition canonique :

F = −kB T ln Ω(N, V, T ) (E.1)

134 LE CALCUL DES ÉNERGIES LIBRES

La connaissance de l'énergie libre absolue d'un système moléculaire ne représente pas

d'intérêt particulier, en revanche la connaissance de l'énergie libre relative, i.e. la dié-
rence d'énergie libre entre deux états est beaucoup plus pertinente. Pour la réactivité
chimique, deux caractéristiques sont essentielles : les grandeurs cinétique et thermody-
namique de réaction. La première caractérise la possibilité pour les réactifs de réagir

∆rF(T) TS

∆rF#

réactifs ∆rF°

produits
C.R.

Fig. E.1  Prol d'énergie libre en fonction d'une coordonnée de réaction.

entre eux et de former une certaine quantité de produit en un temps donné. Cette
donnée cinétique peut être reliée à la diérence d'énergie libre ∆r F ‡ entre l'état de
transition et les réactifs. La seconde, ∆r F o , dénie comme la diérence d'énergie libre
entre les produits et les réactifs, caractérise quant à elle la stabilité thermodynamique
relative des produits par rapport aux réactifs ; elle permet de prédire les proportions
relatives des réactifs et des produits au bout d'un temps de réaction idéalement inni.

A priori , la diérence d'énergie libre entre deux états ξ0 et ξ1 peut être acces-
sible par des simulations de dynamique moléculaire, à condition de simuler un grand
nombre d'événements réactionnels conduisant au passage de ξ0 à ξ1 . Cependant, les
temps de simulation accessibles ne nous permettent pas de simuler un grand nombre
de ces événements rares au regard de la dynamique moléculaire. En eet, la proba-
bilité de trouver le système dans une région au voisinage de l'état de transition est
 135

trop faible. C'est pourquoi plusieurs techniques ont été développées pour estimer ces
diérences d'énergie libre. Parmi celles-ci, on peut citer les méthodes d'intégration
et de perturbation thermodynamiques [35, 12]. On peut aussi citer la méthode um-

brella sampling [116, 117] introduisant un potentiel de biais ou bien encore l'ensemble
blue-moon [118, 119] nécessitant l'introduction de contraintes géométriques. Dans
l'approche umbrella sampling , le potentiel de biais introduit est connu et modie
le potentiel régissant la dynamique du système moléculaire. Le choix d'un potentiel
adapté permet  d'abaisser  articiellement les barrières au cours de la dynamique.
Les diérences d'énergie libre sont alors estimées à partir de la densité d'états is-
sue des trajectoires naturelles sur ce potentiel biaisé, après repondération. On peut
aussi citer les méthodes, apparues plus récemment, de dynamique moléculaire adia-

batique [120, 121] et de méta-dynamique [122, 123], cette dernière ressemblant en de

nombreux points à l'approche umbrella sampling . À partir de deux états ξ0 et ξ1
de l'espace des congurations, ces techniques permettent d'obtenir une  image ther-
modynamique  du chemin de réaction entre ξ0 et ξ1 mais aussi d'estimer la barrière
d'activation associée. Ces techniques nécessitent donc une connaissance préalable du
chemin de réaction ; lorsque la coordonnée de réaction n'est pas bien dénie ou connue,
on peut alors utiliser les approches dites transition path sampling [124].
Au cours de cette thèse, nous nous sommes plus particulièrement intéressés à l'ap-
proche umbrella sampling et celle de l'ensemble blue-moon et nous avons implémen-
tées ces dernières pour pouvoir les utiliser avec l'approche de dynamique moléculaire
développée. Par la suite, seule l'approche de l'ensemble blue-moon est détaillée et
illustrée. Au préalable, les concepts de l'intégration thermodynamique, qui sont à la
base de cette approche, sont brièvement présentés.
136 LE CALCUL DES ÉNERGIES LIBRES

I Intégration Thermodynamique

La méthode de l'intégration thermodynamique utilisait à l'origine comme variable

d'intégration des variables d'état, telles que le volume ou la température [12]. Cette
méthode a par la suite été généralisée à l'utilisation d'autres variables ξ , ces variables
étant des paramètres de l'hamiltonien du système. Ainsi, on peut exprimer l'énergie
libre d'un système restreinte à un paramètre ξ1 :

F (ξ1 ) = −kB T ln Ω(ξ1 ) (E.2)

La fonction de partition canonique Ω(ξ1 ) inclut la restriction à ce paramètre et est

dénie par :
Z Y ³ ³ ´´
Ω(ξ1 ) ∝ ~ I dP~I δ(ξ({R
dR ~ I }) − ξ1 ) exp −βH R,
~ P~ ; ξ1 (E.3)
I

L'énergie libre est une fonction d'état ; par conséquent la diérence d'énergie libre entre
deux états ξ1 et ξ2 ne dépend pas du chemin suivi pour passer de l'un à l'autre, et on
peut écrire cette diérence selon :
Z ξ2
∂F (ξ)
F (ξ2 ) − F (ξ1 ) = dξ (E.4)
ξ1 ∂ξ

On peut alors considérer l'énergie libre comme le potentiel d'une force. Cette force peut
s'exprimer :
∂F (ξ) kB T ∂Ω(ξ)
=− (E.5)
∂ξ Ω(ξ) ∂ξ
À l'aide de la relation (E.3), on peut montrer que l'énergie libre est le potentiel d'une
force moyenne telle que :
Z ξ2 ¿ À
∂H (ξ)
F (ξ2 ) − F (ξ1 ) = dξ ′ (E.6)
ξ1 ∂ξ ξ′

Dans la relation précédente, h. . .iξ′ représente une moyenne statistique, sur un ensemble
à l'équilibre, restreinte à la valeur du paramètre ξ ′ , i.e. à l'hypersurface de l'espace des
n o
phases dénie par ξ({R
~ I }) = ξ ′ . L'expression d'une telle moyenne d'une observable
II Ensemble Blue-Moon 137

³ ´
~
O {R I } , s'exprime par dénition selon :

D ³ ´ ³ ´E
D ³ ´E O {R~ I } δ ξ({R
~ I }) − ξ ′
~I}
O {R = D ³ ´E (E.7)
ξ′ ~
δ ξ({RI }) − ξ ′

La diculté réside maintenant dans l'estimation correcte de cette moyenne. L'ensemble

blue-moon permet d'accéder à cette moyenne statistique via des moyennes tempo-

relles à partir de trajectoires contraintes, i.e. où ~ I })

ξ({R est xé.

II Ensemble Blue-Moon

1
La terre est bleue comme une orange

II.1 Relation entre les moyennes restreintes et contraintes

L'ensemble blue-moon permet d'estimer la moyenne statistique d'une observable

³ ´
restreinte à un paramètre ξ R ~ N donné, à partir de trajectoires où ce paramètre
~ 1 , ..., R

est contraint à la valeur d'intérêt ξ′ :

³ ´
ξ R~ 1 , ..., R
~ N = ξ′ (E.8)

Contraindre ce paramètre à une valeur donnée au cours d'une dynamique impose aussi

que sa dérivée première par rapport au temps soit nulle :

³ ´ X ∂ξ X ∂ξ P~I
ξ˙ R ~N =
~ 1 , ..., R ~˙ I =
·R · =0 (E.9)
∂ ~I
R ∂ ~ I MI
R
I I

La fonction de partition associée à ces trajectoires contraintes adopte donc l'expression

suivante :

Z Y ³ ³ ´´ ³ ´ ³ ´
e ′) ∝
Ω(ξ ~ I dP~I exp −βH R,
dR ~ P~ , ξ δ ξ({R ˙ R
~ I }) − ξ ′ δ ξ({ ~ I }) (E.10)
I

On peut remarquer que cette expression dière de la relation (E.3) correspondant à la

distribution souhaitée. C'est pourquoi la moyenne temporelle d'une observable à partir

1 Paul Éluard, L'Amour, la poésie, 1929

138 LE CALCUL DES ÉNERGIES LIBRES

de trajectoires, où ξ serait contraint, ne correspond pas à la moyenne statistique de

cette même observable à ξ′ donné.

Pour dériver les relations de l'ensemble blue-moon, il faut maintenant eectuer un

changement de variables et passer des coordonnées cartésiennes à un jeu de coordonnées

généralisées. Ces coordonnées généralisées sont dénies par :

~ I }) k = 1, ..., 3N
qk = qk ({R
(E.11)
q1 = ξ

où, par dénition, la première coordonnée q1 correspond au paramètre contraint. Le

lagrangien du système moléculaire peut alors s'exprimer en fonction de ce nouveau jeu

de coordonnées :

1X 2 ³ ´
L = MI R ~˙ I − Ep {R ~I}
2 I
à ! à !
1X X ∂R ~I ~I
X ∂R
= MI q̇k · q̇l − Ep (E.12)
2 I k
∂q k
l
∂q l

1X
= q̇k Gkl q̇l − Ep
2 k,l

où G est la  matrice de passage  associée à ce changement de variables. Les éléments

de cette matrice ont pour expression :

X ∂R ~I
~ I ∂R
Gkl = MI · (E.13)
I
∂qk ∂ql

On peut maintenant dénir les impulsions {pk } associées à ces coordonnées généralisées.
Celles-ci sont dénies telles que :

∂L X
pk = = Gkl q̇l (E.14)
∂ q̇k l

soit encore :
X
q̇k = Zkl pl (E.15)
l

où la matrice Z est la matrice inverse de G :

ZG = I (E.16)
II Ensemble Blue-Moon 139

Les éléments de cette matrice Z ont pour expression :

X 1 ∂qk ∂ql
Zkl = · (E.17)
~ I ∂R
MI ∂ R ~I
I

L'hamiltonien du système peut maintenant s'exprimer uniquement en fonction de ces

coordonnées généralisées :
1X
H = pk Zkl pl + Ep (E.18)
2 k
La fonction de distribution de probabilité souhaitée s'exprime alors :

Z Y
£ ¡ ¢¤
Ω(q) ∝ dpk dqk exp −β pt Zp + Ep δ (q1 − q)
k
Z Y
1
∝ dqk p exp (−βEp ) δ (q1 − q) (E.19)
k
|Z|
Z Y · µ ¶¸
kB T
∝ dqk exp −β Ep + ln |Z| δ (q1 − q)
k
2

où la dernière expression est obtenue après intégration sur les moments conjugués
des coordonnées généralisées. L'expression de la fonction de distribution obtenue à
l'issue de trajectoires contraintes est quant à elle :

Z Y
£ ¡ ¢¤
e
Ω(q) ∝ dpk dqk exp −β pt Zp + Ep δ (q1 − q) δ (q̇1 )
k
Z Y Ã !
£ ¡ t ¢¤ X
∝ dpk dqk exp −β p Zp + Ep δ (q1 − q) δ Z1k pk
k k
Z Y s (E.20)
Z11 £ ¡ ¢¤
∝ dqk exp −β pt Zp + Ep δ (q1 − q)
k
|Z|
Z Y p · µ ¶¸
kB T
∝ dqk Z11 exp −β Ep + ln |Z| δ (q1 − q)
k
2

On peut alors relier la fonction de distribution restreinte à q , dénie par la rela-

tion (E.19), de manière simple à la fonction de distribution générée par des trajectoires
contraintes (E.20).
Rappelons que la moyenne d'une observable O ({qk }), restreinte à q1 = q , s'exprime
par :
140 LE CALCUL DES ÉNERGIES LIBRES

RQ £ ¡ ¢¤
dqk O ({qk }) exp −β Ep + kB2T ln |Z| δ (q1 − q)
k
hOiq = RQ £ ¡ kB T
¢¤ (E.21)
k dq k exp −β E p + 2
ln |Z| δ (q1 − q)
À l'aide de la relation (E.20), la moyenne de cette même observable, obtenue à

partir de trajectoires contraintes, a pour expression :

RQ √ £ ¡ ¢¤
cons
dqk Z11 O ({qk }) exp −β Ep + kB2T ln |Z| δ (q1 − q)
k
hOiq = RQ √ £ ¡ kB T
¢¤ (E.22)
k dq k Z11 exp −β Ep + 2
ln |Z| δ (q1 − q)

On peut nalement relier la moyenne restreinte (E.21) à la moyenne issue de tra-

jectoires contraintes :
D Econs
−1/2
Z11 O
q
hOiq = D Econs (E.23)
−1/2
Z11
q

Ce résultat, initialement démontré en 1989 [118], est à la base de l'ensemble blue-

moon. Pour pouvoir accéder à des diérences l'énergie libre d'après la relation (E.6),

il faut encore estimer la  force thermodynamique  fξ′ , dénie par :

¿ À
∂F ∂H
−fξ′ = = (E.24)
∂ξ ′ ∂ξ ξ′

II.2 Force thermodynamique

L'estimation de la force thermodynamique revient à calculer la dérivée de l'hamil-

tonien par rapport à la contrainte ξ . À l'aide de l'expression (E.18), on peut démontrer

que :

∇ξ F = h∇ξ Ep + kB T ∇ξ ln |J |iξ (E.25)

où |J | représente le déterminant du jacobien associé au passage des coordonnées carté-

siennes aux coordonnées généralisées. On peut démontrer [125, 126] que cette relation

conduit à l'expression suivante :

­ −1/2 ®cons
Z [−λ + kB T G] ξ
∇ξ F = cons (E.26)
hZ −1/2 iξ
III Mise en Place Pratique de Contraintes Géométriques 141

où λ désigne le multiplicateur de Lagrange associé à la contrainte imposée ; sa

détermination est précisée au paragraphe suivant. Z et G ont pour expression :

X 1 µ ∂ξ ¶2
Z= (E.27)
~I
MI ∂ R
I

1 X 1 ∂ξ ∂2ξ ∂ξ
G= 2 · · (E.28)
~ ~ ~ ~J
Z I,J MI MJ ∂ RI ∂ RI RJ ∂ R

Il faut noter que la relation (E.26) a été généralisée dans le cas de contraintes
multiples [127, 128]. Son expression est précisée en annexe III.
Cette relation fait intervenir le multiplicateur de Lagrange associée à la contrainte.
Par la suite la mise en place pratique des contraintes est développée, précisant ainsi la
détermination de ce multiplicateur de Lagrange.

III Mise en Place Pratique de Contraintes Géomé-

triques

On peut distinguer deux grands types de méthode de dynamique moléculaire contrainte :

la méthode analytique et la méthode des paramètres indéterminés [92, 129]. Cette
dernière est essentiellement la méthode analytique modiée de telle façon que les
contraintes soient exactement satisfaites à chaque pas de temps. En eet, la méthode
analytique, en l'absence de schéma correctif, génère des trajectoires où les contraintes
divergent de leurs valeurs imposées [129]. Les équations générales de ces méthodes re-
posent principalement sur deux paramètres : sm et σ . Le premier est relié à l'algorithme
utilisé pour intégrer les équations du mouvement, puisqu'il représente l'ordre maximum
de dérivation des forces par rapport au temps requis par l'algorithme propagateur. Le
second paramètre σ représente la forme générale de la contrainte. En outre, il existe
diérentes techniques pour résoudre les équations issues de ces méthodes, telles que
la méthode matricielle [92] ou bien l'algorithme  SHAKE  [92] et son extension à
142 LE CALCUL DES ÉNERGIES LIBRES

l'algorithme de Verlet aux vitesses, dite  RATTLE  [93]. Nous présentons briève-
ment par la suite la méthode générale des paramètres indéterminés, puis sa mise en
÷uvre pratique avec l'algorithme de Verlet aux vitesses combinée avec la méthode de
résolution RATTLE ; en eet cette approche a été implémentée pour être utilisée avec
les précédents outils décrits dans ce manuscrit.

III.1 Contraintes holonomes

En dynamique moléculaire, les contraintes géométriques choisies pour déterminer

la diérence d'énergie libre le long d'une coordonnée de réaction sont des contraintes
holonomes . Ces l contraintes holonomes σk peuvent se mettre sous la forme générale :
³ ´
~ I (t)} = 0 k = 1, . . . , l
σk {R (E.29)

Une contrainte est dite holonome lorsque son expression ne dépend pas explicitement
des impulsions, i.e. lorsque la forme de Pfa associée à l'expression de la contrainte est
intégrable. Les équations du mouvement des particules, généralisées à la présence de
ces l contraintes holonomes σk , s'expriment alors :

˙
P~I = F~I (t) + G
~ I (t)
l
X (E.30)
= −∇I Ep (t) − λk (t)∇I σk
k=1

où G
~ I (t) représente la force associée aux contraintes et λk désigne les multiplicateurs de

Lagrange associés. La méthode que nous avons choisie et implémentée pour déterminer
ces multiplicateurs de Lagrange est la méthode des paramètres indéterminés [92, 129],
présentée par la suite.

III.2 Méthode des paramètres indéterminés

Un développement en série de Taylor de la solution de la relation (E.30) conduit à

une expression dépendante des forces {F~I (t)}, de leurs dérivées par rapport au temps,
III Mise en Place Pratique de Contraintes Géométriques 143

ainsi que des multiplicateurs de Lagrange {λk } inconnus et de leurs dérivées. Cette
expression adopte la forme suivante :
sX
m +2
~ I (t + δt, {λ(p) (t)}) =R ~˙ I (t) +
~ I (t) + (δt)R 1 ~ (n−2) (δt)n
R FI
n=2
M I n!
sX
m +2
" l n−2
# (E.31)
1 XX p (p) (δt)n
− Cn−2 λk (t) [∇I σk ](n−p−2) (t)
n=2
M I
k=1 p=0
n!

Les termes de la relation précédente peuvent être réordonnés de telle façon que les
termes du même ordre {λ(p) } soient regroupés :
sX
m +2
~ (p) ~ ˙
~ 1 ~ (n−2) (δt)n
RI (t + δt, {λ (t)}) =RI (t) + (δt)RI (t) + F
n=2
MI I n!
~ I (t + δt, {λ(0) (t)}) + δ R
+ δR ~ I (t + δt, {λ(1) (t)})

~ I (t + δt, {λ(2) (t)}) + . . . + δ R

+ δR ~ I (t + δt, {λ(sm −1) (t)})

~ I (t + δt, {λ(sm ) (t)})

+ δR

(E.32)

~ I (t + δt, {λ(0) (t)}) contient sm termes linéaires en {λ(0) (t)}, δ R

où δ R ~ I (t + δt, {λ(1) (t)})

contient sm − 1 termes linéaires en {λ(1) (t)}, etc. Finalement, δ R

~ I (t + δt, {λ(sm ) (t)})

contient un seul terme linéaire en {λ(sm ) (t)}.

La relation (E.32) peut être ré-écrite selon :

~ I (t + δt, {λ(0) (t), . . . , λ(sm −1) (t), λ(sm ) (t)}) =R

R ~ ′ (t + δt, {λ(0) (t), . . . , λ(sm −1) (t)})
I

~ I (t + δt, {λ(sm ) (t)})

+ δR

(E.33)

où le vecteur partiellement contraint R

~ ′ (t+δt, {λ(0) (t), . . . , λ(sm −1) (t)}) est déni comme
I

la somme de tous les termes sauf le dernier de la relation (E.32). Dans le cas où sm = 0,
avec l'algorithme de Verlet par exemple, R
~ ′ (t + δt, {λ(0) (t), . . . , λ(sm −1) (t)}) est réduit
I

à un vecteur purement non contraint. Dans la première étape de la méthode des para-
mètres indéterminés, le jeu de positions {R
~ ′ (t + δt, {λ(0) (t), . . . , λ(sm −1) (t)})} est évalué
I

en calculant uniquement {λ(p) (t)} où p = 0, . . . , sm − 1. Dans une seconde étape, le jeu

144 LE CALCUL DES ÉNERGIES LIBRES

de vecteurs {δ R
~ I } est choisi de manière à satisfaire exactement les contraintes à chaque

pas de temps. En d'autres termes, le jeu {λ(sm ) (t)} est remplacé par un nouveau jeu
de paramètres {γ} :

~ I (t + δt, {λ(0) (t), . . . , λ(sm −1) (t), γ}) =R

R ~ ′ (t + δt, {λ(0) (t), . . . , λ(sm −1) (t)})
I
(E.34)
~ I (t + δt, {γ})
+ δR

Les valeurs de ces nouveaux paramètres {γ} sont déterminées de telle façon que la
relation (E.29) soit exactement satisfaite, soit :

~ I (t + δt, {λ(0) (t), . . . , λ(sm −1) (t), γ})}) = 0 k = 1, . . . , l

σk ({R (E.35)

La comparaison des relations (E.31) et (E.32) conduit à :

l
1 (δt)(sm +2) X
~
δ RI (t + δt, {γ}) = − γk ∇I σk (t) (E.36)
MI (sm + 2)! k=1

Cette relation est l'expression générale pour les déplacements requis à la vérication
~ I évalué au
de la contrainte. De plus, la relation (E.36) montre que ce déplacement δ R
temps (t + δt) est linéaire en {γ} et qu'il dépend des positions {R
~ I } au temps (t).

Le choix d'un algorithme particulier pour l'intégration des équations du mouvement

et la forme des contraintes holonomes déterminent la valeur de sm et la dépendance
fonctionnelle entre les déplacements {δ R
~ I (t+δt, {γ})} et les positions {R
~ I (t)}. On peut

nalement simplier l'expression (E.34) à l'aide de la relation (E.36) pour obtenir :

~ I (t + δt, {λ(0) (t), . . . , λ(sm −1) (t), γ}) =R

R ~ I′ (t + δt, {λ(0) (t), . . . , λ(sm −1) (t)})
l
1 (δt)(sm +2) X (E.37)
− γk ∇I σk (t)
MI (sm + 2)! k=1

La méthode des paramètres indéterminés peut ainsi se résumer en deux étapes. La pre-
mière consiste à évaluer les coordonnées partiellement contraintes et, dans une seconde
étape, les paramètres indéterminés et les coordonnées contraintes. Il faut noter que
cette première étape nécessite le calcul des forces associées aux contraintes et de leurs
dérivées temporelles jusqu'à l'ordre sm , et de même pour les forces {F~I }. Ceci conduit
aux coordonnées partiellement contraintes et les paramètres indéterminés peuvent ainsi
III Mise en Place Pratique de Contraintes Géométriques 145

être évalués. Les contraintes holonomes n'adoptent pas, en général, une forme linéaire
par rapport aux positions des particules. Ainsi, la détermination des paramètres indé-
terminés revient à résoudre un système de l équations non linéaires couplées. Ce type de
système d'équations peut être résolu à l'aide de procédures itératives ; ces procédures
sont elles-mêmes basées sur un développement en série de Taylor de l'expression de la
contrainte holonome suivi d'une linéarisation. Nous présentons par la suite la mise en
÷uvre pratique de la méthode des paramètres indéterminés combinée avec l'algorithme
de Verlet, telle qu'elle a été implémentée durant cette thèse.

III.3 La méthode des paramètres indéterminés combinée avec

l'algorithme de Verlet

L'algorithme de Verlet (cf. partie A Ÿ III.3 p. 22) est du second ordre, il ne nécessite
donc pas les dérivées temporelles des forces, soit sm = 0. Ainsi, les équations pour
l'évolution des coordonnées selon la méthode des paramètres indéterminés combinée
avec l'algorithme de Verlet sont :
l
2 X
~ I′ (t + δt) − (δt)
~ I (t + δt, {γ}) = R
R γk ∇I σk (t) (E.38)
2MI k=1

Cette relation correspond à la relation (E.37) avec sm = 0. Les vecteurs positions

purement non contraints R
~ ′ (t + δt) sont dénis par :
I

2
~ I (t) + (δt) P~I (t) + (δt) F~I (t)
~ I′ (t + δt) = R
R (E.39)
MI 2MI

Le jeu de paramètres {γ} est déterminé de telle façon que les coordonnées au temps
(t + δt) vérient exactement les équations des contraintes. Les nouvelles impulsions
sont dénies selon :
l
δt X
P~I (t + δt, {η}) = P~I′ (t + δt, {η}) − ηk ∇I σk (t + δt) (E.40)
2MI k=1

Les paramètres {η} sont, quant à eux, choisis de telle manière que les impulsions au
temps (t + δt) vérient les équations des contraintes. Ceci est obtenu en dérivant par
146 LE CALCUL DES ÉNERGIES LIBRES

rapport au temps les équations des contraintes (E.29) au temps (t + δt) :

d X 1
~ I (t + δt)}) =
σk ({R ~ I (t + δt)}) = 0
P~I (t + δt) · ∇I σk ({R (E.41)
dt I
M I

On peut utiliser soit la technique d'inversion de matrice ou bien l'approche  SHAKE 

pour résoudre le système de ces l équations couplées (E.41) et obtenir les valeurs de
{η}. La résolution pour {γ} et {η} par la méthode d'inversion de matrice devient
rapidement coûteuse en temps de calcul pour des systèmes avec un grand nombre de
contraintes couplées. Dans de tels cas, il est préférable d'utiliser la procédure SHAKE.
Cette procédure consiste à trouver de manière itérative les déplacements de chaque
particule nécessaires pour satisfaire toutes les contraintes imposées.
La méthode de résolution des jeux d'équations couplées pour {γ} et {η}, combi-
née au propagateur Verlet aux vitesses en utilisant l'approche SHAKE, est dénommée
 RATTLE . La première étape de l'approche RATTLE est analogue à SHAKE :
~ new (t + δt)} obtenues à une itération de cette
les nouvelles positions des particules {R I

procédure sont dénies par :

2
R ~ old (t + δt) − (δt) γ new ∇I σk (t)
~ new (t + δt) = R (E.42)
I I
2MI k
~ old (t + δt)} est donné par la relation (E.39). Ces
où le point de départ des positions {R I

nouvelles positions doivent satisfaire les équations des contraintes, soit :

³n o´ µ½ 2
¾¶
new
~ (t + δt) old
~ (t + δt) − (δt) new
σk R I = σk RI γ ∇I σk (t)
2MI k (E.43)
=0

Le jeu de relations (E.43) n'est pas, en général, linéaire en γknew , même pour des
contraintes simples. Après un développement en série de Taylor de chaque expression
~ old (t + δt)}, cette relation devient :
associée aux contraintes holonomes en {R I
µ½ ¾¶ ³n o´
~ old (δt)2 new ~ Iold (t + δt)
σk RI (t + δt) − γk ∇I σk (t) = σk R
2MI
X (δt)2 ³n o´
− γknew ∇I σk (t) · ∇I σk R ~ old (t + δt)
I
(E.44)
I
2MI

+ ... =0
III Mise en Place Pratique de Contraintes Géométriques 147

où les termes non linéaires ne sont pas explicitement écrits.

Dans un souci d'ecacité calculatoire, tous les termes supérieurs au premier ordre
de la relation (E.44) sont négligés. Ceci se justie naturellement puisque le processus
itératif associé aux contraintes assure que la solution obtenue via cette procédure vérie
la relation non linéaire (E.44). La résolution de cette relation conduit à l'expression
suivante pour γknew :
³n o´
~ old
σk RI (t + δt)
1
γknew = ³n o´ (E.45)
(δt)2 P 1 ∇ σ (t) · ∇ σ R~ old (t + δt)
I 2MI I k I k I

Notons ici que le jeu de paramètres indéterminés {γk } utilisé pour corriger les positions
correspond au multiplicateur de Lagrange nécessaire à la relation (E.26). La validité
de négliger tous les termes non linéaires de la relation (E.44) doit être examinée avec
attention pour chaque forme de contrainte holonome. En eet, la taille de la correction
permise associée à la contrainte sera d'autant plus petite que la non linéarité inhérente
à la contrainte sera  grande , ceci an de justier l'omission des termes non linéaires.
Ce processus itératif est eectué sur l'ensemble des contraintes, jusqu'à ce qu'elles
soient toutes vériées, à une certaine tolérance numérique près. Lorsque toutes ces
dernières sont satisfaites et que les coordonnées au temps (t + δt) sont disponibles, les
forces {F~I (t+δt)} peuvent alors être calculées et utilisées dans une seconde étape. Cette
étape envisage la correction itérative, sur l'ensemble des contraintes, des impulsions.
Leur expression au cours de ce processus est :

δt
P~Inew (t + δt) = P~Iold (t + δt) − ηknew ∇I σk (t + δt) (E.46)
2

où le point de départ des impulsions {P~Iold (t + δt)} est déni par :

" l
#
δt X
P~I′ (t + δt) = P~I (t) + F~I (t) − γk ∇I σk (t) + F~I (t + δt) (E.47)
2 k=1

Les nouvelles impulsions doivent vérier la dérivée par rapport au temps des équations
associées aux contraintes. La combinaison des relations (E.41) et (E.46) conduit à
148 LE CALCUL DES ÉNERGIES LIBRES

l'expression suivante :

X ³n o´ X δt h ³n o´i2
P~Iold (t + δt) · ∇I σk ~ I (t + δt) − ηknew
R ∇I σk R ~ I (t + δt) =0
I I
2
(E.48)

Cette relation est linéaire en ηknew et a pour solution :

P ~ old ³n o´
P (t + δt) · ∇ σ ~
R (t + δt)
2 I I I k I
ηknew = P h ³n o´i2 (E.49)
δt ~
I ∇I σk RI (t + δt)

De la même manière que pour la correction des positions, ce processus itératif par-

court l'ensemble des contraintes imposées au système et ce jusqu'à ce que toutes ces

contraintes soient vériées par les impulsions à une certaine tolérance numérique près.

La détermination des diérents paramètres indéterminés {γk } et {ηk } requiert donc

la connaissance des dérivées des contraintes par rapport aux positions nucléaires. Au

cours de cette thèse, nous avons implémenté diverses contraintes ; leurs détails ainsi

que l'expression de leurs dérivées associées sont précisés en annexe IV.

IV Applications à la Réactivité Chimique

L'approche de l'ensemble blue-moon présentée précédemment permet d'estimer

les diérences d'énergie libre le long d'une coordonnée réactionnelle. Il faut noter

ici que l'étude de réactions chimiques par cette méthode repose sur un choix judi-

cieux et pertinent de cette contrainte. En eet, pour que les diérences d'énergie

libre obtenues par des simulations de dynamique moléculaire soient pertinentes, les

valeurs obtenues selon la relation (E.26), pour chaque valeur discrète de contrainte,

doivent être  convergées . Ceci signie qu'il faut simuler des trajectoires susamment

longues pour que la moyenne temporelle soit équivalente à la moyenne sur l'hypersur-
n ³ ´ ³ ´ o
face ~ I } − ξ = 0, ξ˙ {R
ξ {R ~I} = 0 de l'espace des phases. Un choix inadapté de

cette contrainte peut conduire à des temps de simulation extrêmement grands pour

avoir des valeurs  convergées , rendant ainsi rédhibitoire voire impossible l'étude par

dynamique moléculaire de la réaction choisie.

IV Applications à la Réactivité Chimique 149

Par la suite, nous présentons deux applications de cette méthode.

IV.1 Article  C. Raynaud et al. , accepté à J. Phys. Chem. A

Cet article considère les réactions d'activation de liaisons σ (HH et CH) par un
modèle d'hydrure de lanthanocène : Cl2 LaH . Cette étude nous a en particulier permis
de comparer l'approche de dynamique moléculaire, qui prend en compte de façon expli-
cite les eets de température, avec l'approche couramment utilisée en chimie quantique
pour l'évaluation des grandeurs thermodynamiques. Dans une approche  statique ,
l'estimation des diérences d'énergie libre à température nie repose essentiellement
sur l'approximation harmonique. Cette étude était motivée par la présomption de la
mise en défaut de l'approximation harmonique : en eet les surfaces d'énergies poten-
tielles de ces composés sont le plus souvent  plates . Nous insérons par la suite le
manuscrit de l'article, discutant de cette étude.
150 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 151
152 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 153
154 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 155
156 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 157
158 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 159
160 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 161
162 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 163
164 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 165

Dans l'étude précédente, l'approximation harmonique, couramment utilisée pour

estimer les eets de température, semble être mise en défaut dans le cas où la surface
d'énergie potentielle est assez plate et donc où les eets d'anharmonicité sont marqués.
Nous présentons par la suite une autre étude, considérant des réactions de transfert de
protons intramoléculaires, au sein d'une molécule organique. Cette étude, statique et
dynamique, a fait l'objet d'un article, inséré par la suite.

IV.2 Article  C. Raynaud et al. , accepté à J. Phys. Chem. A

Cet article traite du double transfert de protons au sein de l'acide monothiooxalique.

Dans un premier temps, une étude statique met en évidence un chemin réactionnel plus
favorable que celui qui avait été initialement proposé dans la littérature. Cette étude
statique conduit à un double transfert quasi-concerté où la diérence d'énergie entre les
états de transitions et les intermédiaires réactionnels est très petite. Le prol d'énergie
le long de la coordonnée de réaction présente un caractère de plateau. Dans ce contexte,
ce mécanisme a aussi été étudié d'un point de vue dynamique.
166 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 167
168 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 169
170 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 171
172 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 173
174 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 175
176 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 177
178 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 179
180 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 181
182 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 183
184 LE CALCUL DES ÉNERGIES LIBRES
IV Applications à la Réactivité Chimique 185
186 LE CALCUL DES ÉNERGIES LIBRES
V Estimation de l'Énergie de Point Zéro 187

V Estimation de l'Énergie de Point Zéro

Dans la précédente application concernant le double transfert de protons intramo-

léculaire de l'acide thiooxalique, la contribution d'énergie de point zéro semble être
cruciale. En eet, même s'il est clair que le double transfert de protons, de type 1,4
entre deux hétéroatomes vicinaux, est quasiment concerté, i.e. en une seule étape,
une estimation correcte de cet eet quantique nucléaire est fondamentale. De la même
manière que les eets de température, l'énergie de point zéro est habituellement éva-
luée à l'aide de l'approximation harmonique pour les modes de vibration. Sans pour
autant décrire explicitement les noyaux par une approche quantique, la plupart des
méthodes, se proposant d'estimer l'énergie de point zéro au-delà de l'approximation
harmonique, est basée sur un traitement vibrationnel perturbatif. Cependant, un tel
traitement perturbatif des eets d'anharmonicité, au second ordre par exemple, né-
cessite la connaissance des dérivées, au moins à l'ordre trois, de l'énergie potentielle
par rapport aux positions nucléaires. Ces quantités sont le plus souvent évaluées par
diérence nie des dérivées secondes [130]. Ainsi l'estimation, par ces approches, des
termes cubiques et quartiques de l'énergie de point zéro nécessite, pour un système
moléculaire de N noyaux, 6N − 11 évaluations du hessien. Ces méthodes sont donc
très  gourmandes  en temps de calcul et deviennent rapidement prohibées pour des
systèmes moléculaires où le nombre de degrés de liberté est important.
Au cours de cette thèse, nous nous sommes plus particulièrement intéressés à une
autre approche pour l'évaluation des énergies de point zéro : la méthode dite adiabatic
switching [131]. Cette méthode est, en principe, applicable à n'importe quel système
moléculaire de taille arbitraire possédant un caractère d'anharmonicité plus ou moins
marqué. Cette méthode repose sur une quantication semi-classique, sélectionnant des
trajectoires classiques satisfaisant la condition semi-classique :
µ ¶
1
Ji = ni + h (E.50)
2

où Ji désigne une variable d'action classique et ni représente un nombre quantique

188 LE CALCUL DES ÉNERGIES LIBRES

entier. Cette méthode repose sur l'hypothèse adiabatique d'Ehrenfest, qui stipule que
les nombres quantiques sont conservés au cours d'une transformation lente adiabatique.
L'idée de base de cette approche est que la quantication semi-classique d'un système
anharmonique non séparable peut être eectuée en deux étapes, si l'hamiltonien H du
système peut être exprimé comme la somme de parties séparables et non séparables :
H = H0 +∆H . Ainsi, le mouvement associé à l'hamiltonien séparable H0 est d'abord
quantié, l'hamiltonien régissant l'évolution temporelle du système devient par la suite
dépendant du temps, passant d'une manière douce de H0 à H = H0 + ∆H au
cours d'une longue trajectoire. L'action classique quantiée et associée à l'hamiltonien
H0 est alors transférée adiabatiquement aux quantités correspondantes associées à
l'hamiltonien H .
L'hamiltonien  réel  H permet de décrire le mouvement du système moléculaire
sur la surface d'énergie potentielle complète et s'écrit :
X P~ 2 ³ ´
H = I ~I}
+ V {R (E.51)
I
2MI

La partie séparable H0 représente une approximation de l'hamiltonien réel H et doit

pouvoir décrire le même type de mouvement. Pour cela, on peut choisir l'approximation
harmonique des modes normaux de vibration pour construire cet hamiltonien [132]. Son
expression est alors la somme des termes d'énergie cinétique et potentielle des modes
normaux selon :
X p2 ωi2 qi2
H0 = i
+ (E.52)
i
2 2
où ωi représente la fréquence associée d'un mode normal de coordonnée qi et d'impulsion
pondérée en masse pi . Les mouvements vibrationnels autorisés d'un point de vue semi-
classique sont caractérisés par les variables d'actions Ji et les phases φi , telles que :
I µ ¶
1
Ji = pi dqi = ni + h (E.53)
2
φi ∈ [0, 2π] (E.54)

Les coordonnées et les impulsions associées aux modes normaux sont alors reliées aux
V Estimation de l'Énergie de Point Zéro 189

variables d'actions et aux phases par :

r
ωi Ji
pi = cos φi (E.55)
s π
Ji
qi = sin φi (E.56)
ωi π

L'hamiltonien harmonique peut donc être exprimé selon :

Xµ 1
¶
H0 = ni + ~ωi (E.57)
i
2

Pour un mode propre donné, la connaissance du nombre quantique associé ni et de

la phase φi détermine de manière unique le point de l'espace des phases (qi , pi ). La
quantication du mouvement du système moléculaire apparaît donc ici ; les conditions
initiales de la dynamique sont xées à partir du choix de ces nombres quantiques et
des phases. L'hamiltonien régissant le mouvement du système moléculaire, initialement
harmonique, est dépendant du temps et déni selon :

Ht = H0 + S (t)(H − H0 ) (E.58)

où S (t) désigne la fonction permettant la transformation adiabatique. On peut choisir

1

0,8

0,6
S(t)

0,4

0,2

0
0 0,2 0,4 0,6 0,8 1
t/τ

Fig. E.2  Fonction S (t) utilisée pour rendre l'hamiltonien dépendant du temps.

diverses formes pour cette fonction de passage ; l'une d'entre elles se montre plus par-
ticulièrement ecace [133] quant à la convergence des valeurs propres semi-classiques.
190 LE CALCUL DES ÉNERGIES LIBRES

Cette fonction est représentée sur la gure E.2 et a pour expression :

µ ¶ µ ¶4 µ½· ¸ ¾ ¶
t t t t t
S =− 20 − 70 + 84 − 35 (E.59)
τ τ τ τ τ
où τ représente le temps sur lequel la transformation adiabatique est réalisée. L'hypo-
thèse adiabatique d'Ehrenfest requiert en théorie une transformation inniment lente ;
en pratique un temps de passage correspondant à quelques dizaines de périodes vibra-
tionnelles est susant [134].
La méthode adiabatic switching peut donc être employée pour quantier le mou-
vement vibrationnel sur une surface d'énergie potentielle complète. Les conditions ini-
tiales sont d'abord déterminées et quantiées dans l'approximation harmonique des
modes normaux de vibration, puis la dynamique du système moléculaire est régie par
un hamiltonien dépendant du temps ; cet hamiltonien passe de manière douce de l'ha-
miltonien initial harmonique à l'hamiltonien réel. L'hypothèse adiabatique d'Ehrenfest
suggère ainsi la conservation des nombres quantiques initiaux et donc la quantica-
tion du mouvement vibrationnel nal. Cette approche a démontré son ecacité pour
l'estimation d'énergies vibrationnelles et de point zéro [132, 135].
Au cours de cette thèse, nous avons choisi d'implémenter et d'utiliser cette méthode
pour estimer les énergies de point zéro au-delà de l'approximation harmonique. Plus
particulièrement, cette approche est simple d'utilisation pour les minima locaux des
surfaces d'énergies potentielles. Cependant, son utilisation est plus problématique pour
l'estimation de l'énergie de point zéro au niveau d'un état de transition. En eet, un
état de transition est caractérisé par une courbure négative le long de la coordonnée
réactionnelle et donc par une fréquence associée imaginaire. Ces points particuliers
possèdent donc 3N − 7 degrés de liberté réels intrinsèques, plus un dernier degré cor-
respondant à la coordonnée de réaction. L'utilisation de la méthode de transformation
adiabatique pour ces points ne peut être pertinente que si la coordonnée de réaction
n'est pas couplée aux autres modes de vibration. Or, en pratique, le moindre eet d'an-
harmonicité, combiné avec les erreurs numériques de simulation au cours de la transfor-
mation adiabatique, entraîne le système vers la vallée des réactifs ou celle des produits,
V Estimation de l'Énergie de Point Zéro 191

rendant ainsi caduque la trajectoire. Pour dépasser ce problème, nous proposons de

contraindre au cours de la simulation la coordonnée de réaction, i.e. d'empêcher de
peupler le mode imaginaire. La connaissance des modes normaux de vibration permet
de connaître exactement l'expression du mode imaginaire en fonction des coordonnées
cartésiennes. Il ne reste ensuite plus qu'à interdire tout déplacement le long de cette
coordonnée grâce à la méthode des paramètres indéterminés, présentée précédemment.
Cette approche suppose donc que le couplage entre les diérents modes de vibration
et la coordonnée de réaction ne soit pas trop marqué. Ceci nous laisse penser qu'une
telle approche peut conduire à une estimation relativement satisfaisante de l'énergie de
point zéro, au-delà de l'approximation harmonique, pour les points-selles.

An d'éprouver cette approche, nous avons dans un premier temps considéré le
passage d'un hamiltonien harmonique vers un autre hamiltonien harmonique. Ceci
permet de quantier l'erreur commise et de jauger le schéma de propagation combiné
avec la contrainte du mode imaginaire. La molécule sur laquelle nous avons eectuée
les premiers tests de cette méthode est l'état de transition correspondant à la réaction
d'activation d'hydrogène moléculaire par Cl2 LaH . Le temps de passage considéré ici
pour eectuer la transformation adiabatique est de 5 ps, ce qui est relativement court
vis-à-vis de la plus petite fréquence vibrationnelle de cette molécule. Les équations du
mouvement ont été intégrées avec un pas de temps de 0.1 fs.

Dans le but d'éprouver plus cette méthode, la transformation adiabatique est suivie
par un second passage de l'hamiltonien nal vers l'hamiltonien initial. Ceci permet de
jauger à la fois de la dépendance des conditions initiales et des accumulations d'erreurs
dues à l'intégration numérique des équations du mouvement. Dans ce premier exemple,
l'hamiltonien harmonique nal est relativement proche de l'hamiltonien initial : les
modes propres sont identiques, seules les fréquences dièrent, chacune a été multipliée
par un nombre aléatoire compris entre 0 et 2. La structure initiale correspond à la
structure d'équilibre de l'état de transition, les vitesses initiales sont alors déterminées
d'après les relations (E.53), (E.55) et (E.56). La gure E.3 représente l'évolution au
192 LE CALCUL DES ÉNERGIES LIBRES

cours des deux transformations successives d'un nombre quantique vibrationnel associé
à l'hamiltonien initial.
0,3 0

-0,001
0,2
nombre quantique vibrationnel

Energie totale relative (u.a.)

-0,002
0,1

-0,003

0
-0,004

-0,1
-0,005

-0,2 -0,006
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Temps (fs) Temps (fs)

Fig. E.3  Évolution d'un nombre quantique vibrationnel et de l'énergie totale relative
à l'énergie théorique de point zéro pour l'hamiltonien initial.

On peut remarquer sur la partie gauche de la gure E.3 que ce nombre quantique,
initialement nul, prend des valeurs éloignées de zéro lorsque la dynamique du système
est régie par le second hamiltonien, comme attendu. Néanmoins, ce même nombre quan-
tique reconverge vers zéro lors de la seconde transformation adiabatique. La gure E.3
représente sur sa partie droite l'énergie totale relative à l'énergie totale associée à l'ha-
miltonien initial. De la même manière que le nombre quantique vibrationnel associé
au premier hamiltonien, l'énergie totale reconverge vers la valeur souhaitée après la se-
conde transformation adiabatique. Ces données permettent a priori d'avoir conance
quant à l'ecacité de ces deux transformations successives puisque la valeur initiale de
l'énergie de point zéro est nalement retrouvée.

La gure E.4 représente l'évolution d'un nombre quantique vibrationnel associé au

second hamiltonien harmonique. On peut remarquer à partir de la gure E.4 que ce
nombre quantique vibrationnel, à partir de valeurs non physiques, converge bien vers
zéro comme attendu à la suite de la première transformation, vers 5 ps, pour ensuite
reprendre des valeurs non nulles lors de la seconde transformation.

La gure E.4 représente, sur sa partie droite, l'évolution de l'énergie totale rapportée
V Estimation de l'Énergie de Point Zéro 193

2e-05

0,2
1,5e-05
Nombre quantique vibrationnel

Energie totale relative (u.a.)

0,1

1e-05

5e-06
-0,1

-0,2 0
0 2000 4000 6000 8000 10000 4600 4800 5000 5200 5400
Temps (fs) Temps (fs)

Fig. E.4  Évolution d'un nombre quantique vibrationnel associé à l'hamiltonien har-
monique nal et agrandissement de l'évolution de l'énergie totale relative à l'énergie
théorique totale nale entre 4.5 et 5.5 ps.

à l'énergie théorique de point zéro associée au second hamiltonien. On remarque que

l'énergie totale s'approche, avec une précision très satisfaisante de l'énergie de point
zéro théorique : environ 10−6 u.a.. Ce test semble donc être concluant et bien que le
temps des deux transformations successives soit relativement court, la valeur souhaitée
pour l'énergie de point zéro est obtenue avec une bonne précision.

Pour éprouver plus cette approche, nous l'avons ensuite testée pour deux transfor-
mations successives où, dans ce cas, les modes normaux associés au second hamiltonien
ne sont plus identiques aux premiers. De la même manière que pour le précédent test,
les fréquences ont été multipliées par un nombre aléatoire compris entre 0 et 2. Ceci
permet de tenir compte d'éventuels couplages entre les modes. Dans ce cas, la géomé-
trie initiale ne correspond pas à la structure d'équilibre : les phases φi des diérents
modes normaux ont été choisies de manière aléatoire.

La gure E.5 reporte l'évolution lors des deux passages successifs d'un nombre
quantique vibrationnel associé au premier hamiltonien. De la même manière que pour
le précédent test, on remarque que ce nombre s'éloigne largement de zéro au cours
des deux transformations pour ensuite reconverger vers sa valeur initiale. La gure E.5
représente, sur sa partie droite, l'énergie totale relative à l'énergie totale associée à
194 LE CALCUL DES ÉNERGIES LIBRES

2 0,003

0,0025
1,5
Nombre quantique vibrationnel

Energie totale relative (u.a.)

0,002
1

0,0015

0,5
0,001

0
0,0005

-0,5 0
0 2000 4000 6000 8000 10000 0 2000 4000 6000 8000 10000
Temps (fs) Temps (fs)

Fig. E.5  Évolution d'un nombre quantique vibrationnel et de l'énergie totale relative
à l'énergie totale initiale théorique.

l'hamiltonien initial. Là encore on remarque qu'après les deux transformations adia-

batiques, l'énergie totale reconverge vers la valeur initiale correspondant à l'énergie
de point zéro associée à l'hamiltonien de départ. Dans ce cas, la succession des deux
transformations adiabatiques semble donc être ecace.

3 4,2e-05

2,5
4,1e-05
Nombre quantique vibrationnel

2
Energie totale relative (u.a.)

4e-05

1,5
3,9e-05
1

3,8e-05
0,5

3,7e-05
0

-0,5 3,6e-05
0 2000 4000 6000 8000 10000 4600 4800 5000 5200 5400
Temps (fs) Temps (fs)

Fig. E.6  Évolution d'un nombre quantique vibrationnel associé à l'hamiltonien har-
monique nal et agrandissement de l'évolution de l'énergie totale relative à l'énergie
théorique totale nale entre 4.5 et 5.5 ps.

La gure E.6 représente, comme précédemment, l'évolution d'un nombre quantique

vibrationnel associé au second hamiltonien harmonique. On remarque que ce nombre
vibrationnel converge bien vers zéro après la première transformation adiabatique. An
V Estimation de l'Énergie de Point Zéro 195

de quantier plus précisément l'erreur commise, la gure E.6 représente sur sa partie
droite l'évolution de l'énergie totale rapportée à l'énergie de point zéro associée au
second hamiltonien. De la même manière que pour le premier test, l'énergie totale
converge bien vers la valeur théorique de l'énergie de point zéro associée au second
hamiltonien, avec une très bonne précision : environ 4 × 10−5 u.a..
Ce second test est là encore très satisfaisant ; en eet même si les modes propres ini-
tiaux et naux sont très diérents, la transformation adiabatique permet de transférer
de manière adiabatique la quantication du mouvement vers l'hamiltonien cible.
Les tests qui peuvent naturellement suivre correspondent au passage d'un hamilto-
nien harmonique vers un hamiltonien réel, a priori anharmonique. Toutefois dans ce
cas, les seules possibilités d'estimer la abilité de cette méthode sont d'éprouver dif-
férentes conditions initiales, diérents pas de temps de propagation et enn diérents
temps de passage, puisque l'énergie de point zéro attendue est inconnue. Cet ensemble
de tests n'est malheureusement pas achevé à l'heure de la rédaction de ce manuscrit,
mais les premiers résultats présentés nous laissent penser que cette méthode est une
bonne alternative pour l'estimation de l'énergie de point zéro au-delà de l'approxima-
tion harmonique.
196 LE CALCUL DES ÉNERGIES LIBRES

VI Conclusion
Au regard des capacités informatiques actuelles, l'estimation de grandeurs thermo-
dynamiques par des simulations de dynamique naturelle est aujourd'hui impossible,
puisqu'il faudrait simuler des temps extrêmement longs. Toutefois, à l'aide de trajec-
toires biaisées, l'approche de l'ensemble blue moon , entre autres, permet d'évaluer ces
quantités. L'introduction d'un choix judicieux pour la contrainte représentant la coor-
donnée de réaction permet d'échantillonner correctement l'espace des congurations et
ainsi d'estimer les diérences d'énergie libre associées à un chemin réactionnel. Cette
méthode que nous avons implémentée et combinée avec l'approche de dynamique mo-
léculaire développée au cours de cette thèse permet de rendre compte des grandeurs
thermodynamiques de réaction au-delà de l'approximation harmonique. Plus particu-
lièrement, dans le cas où la surface d'énergie potentielle concernée est  plate , l'esti-
mation usuelle des barrières d'activation à l'aide de l'approximation harmonique peut
conduire à des erreurs relativement conséquentes. L'utilisation de dynamique molécu-
laire ab initio semble alors plus pertinente pour évaluer ces grandeurs. Il faut néanmoins
pondérer ceci en notant que le coût calculatoire de ce type d'approche vis à vis de l'ap-
proximation harmonique habituellement utilisée est énorme. Enn, la détermination
précise de l'énergie de point zéro peut dans certain cas être importante. La méthode de
transformation adiabatique que nous proposons d'utiliser, à la fois pour les minima lo-
caux d'une surface d'énergie potentielle et les points-selles, semble prometteuse. Cette
méthode permettrait, a priori, d'estimer les énergies de point zéro  exactes , à un
niveau de théorie donné, avec une précision satisfaisante. En eet les premiers tests de
cette approche semblent conclure quant à la robustesse et la pertinence de ce schéma.
Conclusion et Perspectives

A u cours de ce travail, nous nous sommes intéressés à la dynamique moléculaire ab

initio, et plus précisément à la mise en ÷uvre de cette technique en utilisant des
fonctions de base locales, gaussiennes, pour la description de la structure électronique
des édices moléculaires. L'approche de la dynamique moléculaire ab initio se situe
entre les méthodes de dynamique quantique et de dynamique moléculaire classique. Ces
premières, qui envisagent une description totalement quantique à la fois des noyaux et
des électrons, sont fortement limitées de par le nombre de degrés de liberté qui peut
être explicitement pris en compte lors de la dynamique. À l'opposé, les techniques de
dynamique moléculaire basées sur une description classique de la structure électronique
par des champs de force, qui permettent de traiter des systèmes de grande dimension,
ne peuvent pas rendre compte de la densité électronique et ainsi décrire la plupart des
mécanismes réactionnels d'intérêt chimique.
La dynamique moléculaire ab initio constitue donc un compromis entre ces deux
approches, en considérant les noyaux comme des particules classiques tout en décrivant
la structure électronique d'un point de vue quantique. Cette technique, qui a connu
un grand essor depuis une vingtaine d'année, est bien souvent synonyme de l'approche
Car-Parrinello. En eet, les approches de type Born-Oppenheimer, qui nécessitent une
convergence complète du problème électronique en chaque point des trajectoires simu-
lées, ne constituent qu'une activité relativement marginale, essentiellement à cause de
leur coût calculatoire rapidement prohibitif. Cependant, l'approche Car-Parrinello est
le plus souvent basée sur une description de la structure électronique sur une base de
198 CONCLUSION

fonctions d'ondes planes dans le cadre de la théorie de la fonctionnelle de la densité.

Ce point de vue pour le traitement de la structure électronique n'est pas toujours le
plus adapté, et dans certains cas on pourrait préférer un traitement par des méthodes
ab initio . Or, ces méthodes ne sont accessibles que dans une base locale d'orbitales
atomiques développées, par exemple, sur un jeu de fonctions gaussiennes.

Ce travail de thèse ne représente que les prémices d'une approche de dynamique

moléculaire ab initio en base locale. Néanmoins, ces premiers résultats semblent en-
courageants. En eet, la technique, fortement inspirée de l'approche Car-Parrinello,
développée durant cette thèse, est a priori viable. Elle permet de diminuer de manière
relativement conséquente l'eort calculatoire nécessaire à ce genre d'approche. Même si
le gain en temps de calcul est modeste par rapport à une approche Born-Oppenheimer,
il reste signicatif au regard du coût énorme des méthodes de dynamique moléculaire
ab initio . Cette approche a été combinée avec diérentes techniques de la dynamique
moléculaire, utiles notamment pour prendre en compte de manière explicite les eets
de température ou encore pour estimer des grandeurs thermodynamiques, telles que
les diérences d'énergie libre. Les quelques illustrations de cette approche conduisent à
des résultats intéressants et montrent que dans certains cas, une approche dynamique
est pertinente pour la description des phénomènes d'intérêt chimique.

Il reste cependant beaucoup à faire pour que la dynamique moléculaire ab initio

en base locale soit véritablement pertinente. Très certainement, le prochain objectif

est d'envisager une description de la structure électronique par une approche multi-
congurationnelle de la fonction d'onde. On pourrait considérer un schéma analogue
à celui proposé dans ce manuscrit, i.e. une accélération de la convergence des orbi-
tales moléculaires, de chaque conguration, à l'aide d'une propagation classique de ces
dernières, combinée avec une propagation classique des coecients pour le développe-
ment multidéterminantal. En eet, même si ce genre de technique est, encore à l'heure
actuelle, très  gourmand  en temps de calcul, des méthodes de localisation [136],
par exemple, permettent de réduire signicativement l'eort calculatoire. De plus, une
CONCLUSION 199

approche multicongurationnelle par fragment [137] serait dans ce cadre très utile et
permettrait la réalisation de simulations de dynamique moléculaire. De plus ce type
de méthodologie, en particulier appliqué à des états excités [138] et combiné avec des
approches de dynamique moléculaire, serait bien adapté à l'étude de mécanismes ré-
actionnels photochimiques [139]. D'une manière générale, la dynamique moléculaire ab
initio combinée avec une méthode de description de la structure électronique précise,
éventuellement diérente du cadre de la théorie de la fonctionnelle de la densité, est
une perspective qu'il ne faut pas négliger.
D'un autre côté, la description complète par des méthodes quantiques en base locale
de systèmes de grande dimension est encore inaccessible. Toutefois, les méthodes de
type  QM/MM  représentent une alternative le plus souvent pertinente pour ce type
d'études. En eet, pour la plupart des cas d'intérêts chimiques, le système étudié est
bien souvent constitué d'une partie  active  de taille relativement réduite, l'environ-
nement est en quelque sorte  spectateur  du phénomène considéré, même s'il joue
un rôle prépondérant. C'est, par exemple, le cas de la photoisomérisation du rétinal
au sein de la bactériorhodopsine [140]. Dans ce cadre, les premiers tests du schéma
de dynamique moléculaire développé durant cette thèse, combiné avec une description
hybride, de type QM/MM, laissent la voie entre-ouverte au traitement de systèmes de
plus grande dimension et éventuellement la simulation de phénomènes chimiques en
phase condensée.
Finalement, le travail accompli semble bien maigre vis-à-vis de toutes les possibilités
que les approches de dynamique moléculaire ab initio en base locale laissent présager,
mais nous ne doutons pas des futures avancées dans ce domaine.
Annexe I
Systèmes à Échelles de Temps
Multiples
Le choix du pas de temps d'intégration est déterminé par la nature des forces
agissant sur le système. Dans les systèmes moléculaires, les forces résultent de diérentes
classes d'interactions entre particules et génèrent des mouvements dont les échelles de
temps caractéristiques sont disparates. Toutefois, le pas de temps doit être choisi tel que
le mouvement le plus rapide du système puisse être intégré de façon stable et précise.
Ceci conduit à des procédures numériques inecaces puisque les forces associées aux
mouvements les plus lents sont recalculées sur de petites échelles de temps sans qu'elles
soient notablement modiées. Ce type de problème à échelles de temps multiples est
quasiment toujours présent au cours des simulations de dynamique moléculaire, soit
par la nature même du système ou bien par la présence d'un  réservoir  de chaleur
tel qu'un thermostat.

Les méthodes permettant d'accélérer l'intégration des équations du mouvement sont

basées sur la dénition d'un  sous-système  rapide. Les équations du mouvement
associées au sous-système rapide sont intégrées sur des petits pas de temps alors que
les équations du mouvement associées aux degrés de libertés  lents  du système sont
quant à elles intégrées sur des pas de temps plus grands. Ces approches peuvent être
développées dans le cadre de la factorisation de Trotter de l'opérateur de Liouville.
202 ANNEXES

Considérons le cas où les degrés de liberté d'un système peuvent être séparés en deux
groupes, l'un  rapide  et l'autre  lent , notés respectivement x et y . Le système est
donc constitué d'éléments  lourds  (le sous-système y ) et  légers  (le sous-système
x). Nous noterons de manière abrégée x et px les positions et les impulsions associées
au sous-système x et de la même façon les variables du sous-système y . Décomposons
l'opérateur de Liouville :
iL = iLx + iLy (I.1)

où
∂ ∂
iLx = ẋ + Fx (x, y) (I.2)
∂x ∂px
∂ ∂
iLy = ẏ + Fy (x, y) (I.3)
∂y ∂py

À l'aide de la factorisation de Trotter, un propagateur discret en temps peut être déni :

µ ¶ µ ¶
∆t ∆t
Gxyx (∆t) = exp iLx exp (iLy ∆t) exp iLx (I.4)
2 2

Cette factorisation fait apparaître la première et la dernière composante du propaga-

teur, qui sont associées au mouvement rapide du sous-système x alors que la composante
du milieu est associée au mouvement lent du sous-système y . Le propagateur  rapide 
peut être lui aussi factorisé :
µ ¶ · µ ¶ µ ¶ µ ¶¸n/2
∆t δt ∂ ∂ δt ∂
exp iLx = exp Fx exp δtẋ exp Fx (I.5)
2 2 ∂px ∂x 2 ∂px

où δt = ∆t/n représente donc le petit pas de temps associé au mouvement rapide du

sous-système x. Le propagateur du sous-système lent y peut lui aussi être factorisé :
µ ¶ µ ¶ µ ¶ µ ¶
∆t ∆t ∂ ∂ ∆t ∂
exp iLy = exp Fy exp ∆tẏ exp Fy (I.6)
2 2 ∂py ∂y 2 ∂py

Si ces factorisations sont utilisées dans la relation (I.4), alors les degrés de liberté
rapides et leurs impulsions associées {x, px } seront déterminés numériquement à l'aide
de l'algorithme Verlet aux vitesses pour n pas de temps δt, alors que les degrés de
liberté lents seront déterminés à l'aide du même propagateur sur un seul pas de temps
 large  ∆t. Finalement, l'application de Gxyx peut se décomposer en trois étapes :
ANNEXES 203

 à partir de l'état initial {x(0), y(0), px (0), py (0)}, l'intégrateur Verlet aux vitesses
est utilisé n/2 fois avec un pas de temps δt = ∆t/n pour générer l'état au temps
∆t/2 des variables du sous-système x,
 les variables du sous-système y sont déterminées au temps (∆t) à l'aide d'un seul
pas de temps ∆t et du propagateur Verlet aux vitesses,
 le propagateur Verlet aux vitesses est utilisé n/2 fois avec un pas de temps δt =
∆t/n pour déterminer les variables du sous-système x à la date ∆t.
Les forces des coordonnées lentes ne sont donc calculées qu'une seule fois par pas de
propagation ∆t = nδt. Si la dimensionnalité du sous-système rapide est petite vis à
vis de la taille totale du système entier, les forces agissant sur les variables lentes sont
donc recalculées moins souvent que dans le cas des méthodes habituelles.
Annexe II

Intégrateurs d'Ordres Supérieurs

À partir de propagateurs d'un certain ordre, par exemple l'algorithme de Verlet

aux vitesses, des propagateurs d'ordre supérieur peuvent être générés sans pour autant

nécessiter l'évaluation des dérivées des forces par rapport aux positions. Le schéma

d'intégration de Yoshida-Suzuki [104, 105] est notamment très utile pour intégrer les

équations du mouvement des chaînes de thermostats de Nosé-Hoover [44, 91].

Supposons que l'opérateur de Liouville puisse s'écrire sous la forme générale sui-

vante :
M
X
iL = iLk (II.1)

k=1

Le développement symétrique selon la factorisation de Trotter de l'opérateur d'évolu-

tion conduit au propagateur suivant :

µ ¶ µ ¶
(1) δt δt
U (δt) = exp iL1 × ... × exp iLM −1
2 2
× exp (iLM δt) (II.2)

µ ¶ µ ¶
δt δt
× exp iLM −1 × ... × exp iL1
2 2

En particulier, Suzuki a montré que tous les propagateurs de la forme U (2m) (δt) satis-

font la relation de récurrence suivante :

U (2m) (δt) = U (2m−1) (δt)

(II.3)
£ ¤2 £ ¤2
= U (2m−3) (pm δt) × U (2m−3) ((1 − 4pm )δt) × U (2m−3) (pm δt)
206 ANNEXES

où les pm représentent les coecients de Yoshida-Suzuki, dénis par :

1
pm = (II.4)
4− 41/(2m−1)

Ainsi, à partir de la relation (II.2), l'équation (II.3) peut être utilisée pour générer des
factorisations d'ordre supérieur correspondant à des propagateurs d'ordre plus élevé.
La relation (II.3) révèle de plus que les propagateurs d'ordre pair et impair sont équi-
valents, et en particulier que U (2) (δt) = U (1) (δt), expliquant ainsi pourquoi le déve-
loppement (II.2) est précis en O ((δt)3 ). En prenant par exemple m = 3, on obtient un
algorithme d'intégration précis au cinquième, i.e. sixième, ordre.
On peut nalement combiner le schéma d'intégration de Yoshida-Suzuki et un dé-
veloppement en pas de temps multiples, ceci conduit à la relation suivante :
µ ¶ n
"m µ ¶#
δt Y Y wj δt
exp iLk = exp iLk (II.5)
2 i=1 j=1
2n

où les wj sont les paramètres d'intégration de Yoshida-Suzuki  simpliés , dépendant

de m. Le tableau 1 suivant regroupe quelques valeurs de ces coecients en fonction de
m.

m wj
1 1
3 w1 = w3 = 1/(2 − 21/3 ), w2 = 1 − 2w1
7 w1 = w7 = −1.17767998417887
w2 = w6 = 0.235573213359357
w3 = w5 = 0.78451361047756
w4 = 1 − 2(w1 + w2 + w3 )

Tab. 1  Paramètres d'intégration de Yoshida-Suzuki.

Annexe III
Force Thermodynamique et
Contraintes Multiples
Dans cette annexe, nous précisons l'expression de la force thermodynamique dans le
cas de contraintes multiples. Cette approche permet d'évaluer des diérences d'énergie
libre dans le cas où la coordonnée de réaction est  décrite  par plusieurs contraintes.
L'énergie libre F restreinte aux paramètres {ξ1 , . . . , ξp } s'exprime en fonction de la
distribution associée :

F (ξ1 , . . . , ξp ) = −kB T ln Ω(ξ1 , . . . , ξp ) (III.1)

Coordonnées généralisées
Nous supposons qu'un jeu de M − p fonctions (q1 , . . . , qM −p ) puisse être déni de
telle manière que le jeu (ξ1 , . . . , ξp , q1 , . . . , qM −p ) forme un jeu complet de coordonnées
généralisées pour le système moléculaire à N particules (M = 3N ). La dérivée partielle
par rapport à ξi de l'énergie libre selon la relation (III.1) s'écrit selon :

∂F kB T ∂Ω
=− (III.2)
∂ξi Ω ∂ξi
208 ANNEXES

La distribution Ω pour un ensemble canonique, restreinte à ces paramètres, s'exprime

selon :
Z Y µ ¶Y
H
Ω(ξ1⋆ , . . . , ξp⋆ ) ∝ ~ ~
dRI dPI exp − ~ I }) − ξ ⋆ )
δ(ξi ({R i (III.3)
I
kB T i

où les ξi⋆ représentent les valeurs auxquelles on souhaite se restreindre. Le jacobien

associé à la transformation des coordonnées cartésiennes en coordonnées généralisées
est déni selon :  
Jξ
J ≡  (III.4)
Jq
où Jξ désigne les p premières lignes du jacobien, et Jq les suivantes. On peut dénir la
matrice Z :
Z ≡ J M −1 J t (III.5)

où M représente le tenseur des masses. Cette matrice peut se décomposer en sous-blocs :

 
Zξ Zξq
Z=  (III.6)
Zqξ Zq

où Zξ est une matrice p × p, Zξq est de dimension p × (M − P ) et Zq est de dimension

(M − p) × (M − p). La matrice inverse de Z est notée G, et de la même manière :
 
Gξ Gξq
G=  (III.7)
Gqξ Gq

À l'aide de ces coordonnées généralisées, l'hamiltonien du système prend la forme sui-

vante :
1 1
H (ξ, q, pξ , pq ) = ptξ Zξ pξ + ptq Zq pq + ptξ Zξq pq + Ep (ξ, q) (III.8)
2 2
La combinaison des relations (III.2) et (III.3) avec ces coordonnées généralisées conduit
à l'expression suivante :
R ³ ´
dqdpq dpξ ∂H
∂ξi
exp − kH
BT
∇ξi F = R ³ ´ (III.9)
dqdpq dpξ exp − kH
BT
ANNEXES 209

Ainsi, la relation (III.2) peut être ré-écrite selon :

¿ À
∂F ∂H
= (III.10)
∂ξi ∂ξi ξ

Après diérenciation de la relation (III.8) et intégration sur les impulsions pξ et pq , on

peut obtenir une nouvelle expression pour la relation (III.9) :

∇ξ F = h∇ξ Ep + kB T ∇ξ ln |J |iξ (III.11)

La dérivée de l'énergie libre peut alors être vue comme le résultat de deux contributions :
la force mécanique agissant sur ξ et les variations de l'élément de volume associé au
système de coordonnées généralisées.

Force thermodynamique
Par la suite, nous utiliserons le fait que, pour un jeu donné {ξi⋆ }, il soit possible de
choisir une base q telle que :
Zqξ (ξ ⋆ , q) = 0 ∀q (III.12)
La relation (III.10) dépend explicitement du choix de toutes les coordonnées générali-
sées, incluant notamment q . Nous allons ici modier cette relation pour se libérer de
cette dépendance et la rendre indépendante du choix de q . Ceci peut être eectué en
intégrant analytiquement le plus de termes possibles de la relation (III.10).
Commençons par introduire les notations simpliées suivantes :
p 1 ∂
x′i ≡ ~I
MI R ∇′i ≡ √
~I
MI ∂ R
(III.13)
P~I
≡√p′xi
MI
L'équation d'évolution pour la variable pξi est la suivante :
dpξi ∂H
=− (III.14)
dt ∂ξi
Le moment conjugué pξ est déni comme la dérivée du lagrangien par rapport à ξ˙,
soit :
X dξj X dqk
pξi ≡ [Gξ ]ij + [Gξq ]ik (III.15)
j
dt k
dt
210 ANNEXES

Nous pouvons diérencier la relation (III.15) par rapport au temps et utiliser la rela-
tion (III.14) pour obtenir une expression pour ∂H
∂ξi
:

∂H dpξ X d [Gξ ]ij dξj d2 ξj

=− i =− − [Gξ ]ij 2
∂ξi dt dt dt dt
j
X d [Gξq ] dqk
(III.16)
ik d2 qk
− − [Gξq ]ik 2
k
dt dt dt

Puisque le jeu {q} vérie la relation (III.12), le dernier terme de la relation (III.16) est
égal à zéro. Cette relation peut ensuite se simplier en utilisant la règle de la chaîne
pour obtenir :
X£ · ¸
∂H ¤ d2 ξj X £ −1 ¤ ∂Zξ ′ X d [Gξq ] dqk
=− −1
Zξ ij 2 + Zξ ij ′
· p x p ξk
− ik
(III.17)
∂ξi j
dt jk
∂x jk k
dt dt

En exprimant p′x en fonction de pξ et de pq , on peut démontrer que le deuxième terme

de la relation (III.17) est égal à :
X£ · ¸ X£
−1
¤ ∂Zξ ′ ¤ ∂ [Zξ ]jk £ ′ ¤
Zξ ij ′
· p x p ξk
= Zξ−1 ij ′
Jξ rl pξr pξk
jk
∂x jk jklr
∂x l

X£
(III.18)
¤ ∂ [Zξ ]jk £ ′ ¤
+ Zξ−1 ij ′
Jξ r+p,l pqr pξk
jklr
∂x l

où J ′ est déni de manière analogue à J mais avec x′ à la place de x. Dans cette

dernière relation, nous avons séparé la partie droite en contributions paire et impaire
de pξ et pq .
Rappelons que nous voulons calculer :
Z µ ¶
H ∂H
dpq dpξ exp − (III.19)
kB T ∂ξi
³ ´
Puisque la base q choisie vérie la relation (III.12), la fonction exp − kHB T est paire
en pξ et pq . Ainsi dans la relation (III.18), toutes les contributions impaires en pξ et pq
s'annulent. Ceci conduit après intégration sur les impulsions pξk à l'expression suivante :
* · ¸ + * +
X£ ¤ ∂Z X £ ¤ ∂ [Z ξ ] £ ¤ ∂ξ
Zξ−1 ij
ξ
′
· p′x pξk = kB T Zξ−1 ij ′
jk
Zξ−1 kr ′
r
(III.20)
jk
∂x jk jklr
∂x l ∂x l
ξ ξ
ANNEXES 211

1 ∂qr
De plus les vecteurs vérient la relation (III.12), ils sont donc orthogonaux aux
Ms ∂xs

vecteurs ∇ξ1 , . . . , ∇ξp . Par conséquent le troisième terme de la partie droite de la rela-
tion (III.17) ne contribue pas à h∇ξ H iξ . Finalement, l'insertion de la relation (III.20)

dans l'équation (III.17) conduit à la relation suivante :

* + ¿ À
X 1 2
−1 −1 −1 d ξ
h∇ξ H iξ = kB T Z · ∂l Zξ · Zξ · ∇ξ − Zξ (III.21)
l
Ml ξ dt2 ξ
ξ

∂Zξ
où nous notons ∂l Zξ = ∂xl
. En notant λ le multiplicateur de Lagrange associé de la

procédure SHAKE, on obtient l'expression suivante :

* +
X 1
h∇ξi H iξ = −λξi + kB T Zξ−1 · ∂l Zξ · Zξ−1 · ∂l ξ (III.22)
l
M l
ξ

À l'aide de la relation de l'ensemble blue-moon, on peut nalement exprimer la dérivée

de l'énergie libre selon :

D ³ P £ −1 ¤ ´Econs
kB T
|Zξ |1/2 −λξi + 2 l Zξ il
(∇ ′
ξl · ∇ ′
ln |Zξ |)
ξ
∇ξi F = cons (III.23)
h|Zξ |−1/2 iξ

Cette relation a été démontrée par W.K. den Otter [127, 128]. Cette expression générale

permet de retrouver la formule présentée dans le cas d'une seule contrainte [125, 126] :

­ −1/2 ®cons
Z [−λ + kB T G] ξ
∇ξ F = cons (III.24)
hZ −1/2 iξ

où nalement Z et G ont pour expression :

X 1 µ ∂ξ ¶2
Z= (III.25)
~I
MI ∂ R
I

1 X 1 ∂ξ ∂2ξ ∂ξ
G= · · (III.26)
2 ~ I ∂R
Z I,J MI MJ ∂ R ~ J ∂R
~IR ~J
Annexe IV
Dérivées de Contraintes Holonomes
L'utilisation de contraintes géométriques holonomes σ en dynamique moléculaire

impose la connaissance du gradient de celle-ci par rapport aux positions atomiques. De

plus, lorsque cette contrainte est utilisée comme coordonnée de réaction dans l'ensemble

blue-moon, la connaissance des dérivées secondes est aussi requise pour déterminer la

 force  associée, i.e. la dérivée de l'énergie libre par rapport à cette contrainte. Nous

préciserons ici l'expression des termes correctifs Z et G, dénis selon :

X 1 ∂σ ∂σ
Z= · (IV.1)
~ I ∂R
MI ∂ R ~I
I

1 X 1 ∂σ ∂ 2σ ∂σ
G= 2 · · (IV.2)
~ ~ ~ ~J
Z I,J MI MJ ∂ RI ∂ RI ∂ RJ ∂ R

Dans cette annexe les dérivées premières et secondes des contraintes implémentées lors

de cette thèse sont présentées. Il faut noter que chacune des contraintes présentées

par la suite a été implémentée et généralisée pour les barycentres de masse de groupes

d'atomes. En eet, le passage d'une dérivée par rapport à un barycentre à la dérivée par

rapport à un atome est trivial ; soit ~G

R le barycentre de masse d'un ensemble d'atomes

~I}
{R et en notant MG la masse du barycentre :

∂ MI ∂
= (IV.3)
~I
∂R ~G
MG ∂ R
214 ANNEXES

Contrainte de distance
La distance entre deux centres A et B est la contrainte la plus simple que l'on puisse
~ A et R
imaginer. Les positions des centres A et B sont respectivement notées R ~ B . Cette

A B
d

Fig. 7  Contrainte de distance entre deux centres

contrainte peut ainsi s'écrire :

³ ´2
σ= R~A − R
~ B − d2 (IV.4)

où d désigne la distance imposée. Les dérivées premières de cette contrainte sont im-
médiates et s'expriment selon :

∂σ ³ ´
=2 R ~B
~A − R (IV.5)
∂R~A
∂σ ³ ´
= −2 R ~A − R
~B (IV.6)
∂R~B

De la même manière, les dérivées secondes sont triviales :

∂ 2σ ∂ 2σ
= =2 (IV.7)
~2
∂R ~2
∂R
A B

Finalement, l'expression des termes correctifs Z et G pour le calcul des énergies libres
est :
µ ¶2 µ ¶2
1 ∂σ 1 ∂σ
Z= +
MA ∂ R ~A MB ∂ R ~B
(IV.8)
MA + MB 4d2
= 4d2 =
MA · MB µ
où µ désigne la masse réduite de A et B .

1
G= (IV.9)
2d2
ANNEXES 215

Contrainte d'angle

Une contrainte d'angle θ est dénie par trois centres A, B et C . Nous adopterons
la dénition de la gure 8.

d1

θ
B d2 C

Fig. 8  Contrainte d'angle dénie par trois centres

La contrainte peut donc être dénie selon :

à !
~A − R
(R ~ B ) · (R
~C − R
~ B)
σ = arccos −θ (IV.10)
d1 d2

où θ désigne la valeur d'angle imposée. En dénissant les vecteurs normés suivants :

~ ~
~ BA = RA − RB
E (IV.11)
~A − R
kR ~ Bk
~ ~
~ BC = RC − RB
E (IV.12)
~C − R
kR ~ Bk

~ A conduit à :
La dérivée première de la contrainte par rapport à R

∂σ ~ BA − E
cos θE ~ BC
= (IV.13)
∂R~A d1 sin θ

Les centres A et C étant interchangeables, on obtient l'expression analogue pour la

~C :
dérivée première par rapport à R

∂σ ~ BC − E
cos θE ~ BA
= (IV.14)
∂R~C d2 sin θ

~ B s'exprime quant à elle :

La dérivée par rapport à R

∂σ ~ BA + (d2 − d1 cos θ)E

(d1 − d2 cos θ)E ~ BC
= (IV.15)
∂R~B d1 d2 sin θ
216 ANNEXES

Les dérivées secondes sont beaucoup plus compliquées à exprimer, elles ne sont donc
pas détaillées ici. On donne néanmoins le résultat pour les termes correctifs Z et G .

1 d21 + d22 − 2d1 d2 cos θ 1

Z= + + (IV.16)
MA d21 2 2
MB d1 d2 MC d22

On peut exprimer G simplement en fonction de Z :

4 sin θ
G=Z (IV.17)
MB d1 d2

Contrainte d'angle dièdre

De la même manière que précédemment, un angle dièdre est dénie à partir de

quatre centres. Nous prendrons les notations décrites par la gure 9.

A A

θ1 d1
d1 B
C
d2 Φ

B d3
C d3
θ2 D
D

Fig. 9  Contrainte d'angle dièdre dénie par quatre centres

De la même manière que précédemment, on peut dénir les vecteurs normés suivant :

~ ~
~ BA = RA − RB
E (IV.18)
~A − R
kR ~ Bk
~ ~
~ BC = RC − RB
E (IV.19)
~C − R
kR ~ Bk
~ CB = −E
E ~ BC (IV.20)
~ ~
~ CD = RD − RC
E (IV.21)
~D − R
kR ~ Ck
ANNEXES 217

La contrainte peut donc s'exprimer selon :

³ ´ ³ ´
~ BA × E
E ~ BC · E ~ CB × E
~ CD
σ = arccos  −Φ (IV.22)
sin θ1 sin θ2

où Φ désigne la valeur d'angle dièdre imposée. Après quelques développements ma-

thématiques assez fastidieux, on obtient les expressions suivantes pour les dérivées
premières :
∂σ ~ BA × E
E ~ CB
= (IV.23)
∂R~A d1 sin2 θ1
∂σ d2 − d1 cos θ1 ³ ~ ´
~ BC + cos θ2
³ ´
= E BA × E E ~ CB
~ CD × E (IV.24)
∂R~B d1 d2 sin2 θ1 d2 sin2 θ2
∂σ d2 − d3 cos θ2 ³ ~ ´ cos θ1 ³ ~ ´
= ~
ECD × ECB + ~
EBA × EBC (IV.25)
∂R~C d3 d2 sin2 θ2 d2 sin2 θ1

∂σ ~ CD × E
E ~ BC
= 2 (IV.26)
∂R~D d3 sin θ2
Les dérivées secondes ne seront pas développées ici, leurs expressions tenant sur plu-
sieurs pages ! Toutefois, pour  simplier  l'écriture des termes Z et G , dénissons :

d2 − d1 cos θ1
X1 = (IV.27)
d1 d2 sin θ1
d2 − d3 cos θ2
X2 = (IV.28)
d2 d3 sin θ2

et :

cos θ1
Y1 = (IV.29)
d2 sin θ1
cos θ2
Y2 = (IV.30)
d2 sin θ2

L'expression de Z est alors :

µ ¶2 µ ¶2
1 1 1 1
Z= +
MA d1 sin θ1 MD d3 sin θ2
1 ¡ ¢
+ X12 + Y22 − 2X1 Y2 cos Φ (IV.31)
MB
1 ¡ 2 ¢
+ X2 + Y12 − 2X2 Y1 cos Φ
MC
218 ANNEXES

L'expression de G est quant à elle :

µ ¶
X1 Y2 X2 Y1
G =4Z + sin Φ
MB MC
2 £ ¤
+ 2 2 2X1 Y2 − (X12 + Y22 ) cos Φ sin Φ
d2 MB
2 £ ¤
+ 2 2 2X2 Y1 − (X22 + Y12 ) cos Φ sin Φ
d2 MC
2 £¡ 2 ¢ ¤
− X1 + Y22 − 2X1 Y2 cos Φ X2 Y1 sin Φ
M B MC
2 £¡ 2 ¢ ¤
− X2 + Y12 − 2X2 Y1 cos Φ X1 Y2 sin Φ (IV.32)
M B MC
2 d1 − d2 cos θ1
− [X2 (X1 − Y2 cos Φ)] sin Φ
MB MC d1 d22 sin2 θ1
2 d1 − d2 cos θ1
− [Y2 (Y1 − X2 cos Φ)] sin Φ
MB MC d1 d22 sin2 θ1
2 d3 − d2 cos θ2
− [X1 (X2 − Y1 cos Φ)] sin Φ
MB MC d3 d22 sin2 θ2
2 d3 − d2 cos θ2
− [Y1 (Y2 − X1 cos Φ)] sin Φ
MB MC d3 d22 sin2 θ2

Contrainte de  projection 
Une contrainte de  projection  est une contrainte dénie à l'aide de trois centres.

d1

θ
B d2 C

Fig. 10  Contrainte de projection dénie par trois centres

À partir de la dénition de la gure 10, cette contrainte revient à xer la projection

de la distance entre A et B sur la distance entre B et C . Cette contrainte a donc pour
expression :
d1 cos θ
σ= −λ (IV.33)
d2
ANNEXES 219

où λ désigne la valeur imposée de cette projection. Les dérivées premières de cette

contrainte ont pour expressions :

∂σ 1 ³~ ´
= 2 R ~
C − RB (IV.34)
∂R~A d2

∂σ 2d1 ³ ´ ³ ´
= 3 cos θ R ~ B + 1 2R
~C − R ~B − R ~A
~ CR (IV.35)
∂R~B d2 d22
∂σ 2d1 ³ ´ ³ ´
= 3 cos θ R~B − R~C + 1 R ~A − R
~B (IV.36)
∂R~C d2 d22
Ces dérivées premières conduisent à l'expression de Z suivante :
· ¸
1 d21 1 d21 d1 d21
Z= 2 + + 2
1 + 2 + 6 cos θ + 8 2 cos θ (IV.37)
d2 MA d42 MB d22 MC d2 d2 d2

Les expressions des dérivées secondes pour cette contrainte sont comme précédemment
assez complexes. Le facteur G peut néanmoins se  simplier  et conduire à l'expres-
sion :
2 d1 cos θ 1 2d31 cos θ
G= +
MA MC d52 MC2 d72
2 1 h³ ~ ~
´ ³
~ ~
´i
− R C − R B · R C − R A
MA MB d62
(IV.38)
1 2 ³~ ~
´3 ³
~C − R ~B
´
− 2 8 R A − RC · R
MB d2
2 d1 cos θ ³ ~ ´ ³
~A · R ~B + R ~ C − 2R
~A
´
+ R C − R
MB MC d72
Bibliographie
[1] E. Schrödinger. Phys. Rev. 28 (1926), 1049.

[2] G. N. Lewis. J. Am. Chem. Soc. 38 (1916), 762.

[3] R. J. Gillespie. Quaterly Rev. 11 (1957), 339.

[4] R. Hoffmann. Angew. Chem., Int. Ed. Engl., Nobel Lecture 21 (1982), 711.

[5] K. Fukui. Angew. Chem., Int. Ed. Engl., Nobel Lecture 21 (1982), 801.

[6] M. H. Beck, A. Jäckle, G. A. Worth et H.-D. Meyer. Physics Reports

324 (2000), 1.

[7] R. Car et M. Parrinello. Phys. Rev. Lett. 55, 22 (1985), 2471.

[8] D. Marx et J. Hutter. Modern methods and algorithms of quantum chemistry,

vol. 1 de NIC Series. J. Grotendorst, John von Neumann Institute for Computing,
Jülich, 2000.

[9] M. J. Field. Chem. Phys. Lett. 172 (1990), 83.

[10] G. Lippert, J. Hutter et M. Parrinello. Theor. Chem. Acc. 103 (1999),

124.

[11] S. S. Iyengar, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scu-

seria et M. J. Frisch. J. Chem. Phys. 115, 22 (2001), 10291.

[12] D. Frenkel et B. Smit. Understanding Molecular Simulation  From Algo-

rithms to Applications. Academic Press, San Diego, 1996.

[13] W. Kolos. Adv. Quant. Chem. 5 (1970), 99.

[14] W. Kutzelnigg. Mol. Phys. 90 (1997), 909.

222 BIBLIOGRAPHIE

[15] P. A. M. Dirac. Proc. Cambridge Phil. Soc. 26 (1930), 376.

[16] E. Deumens, A. Diz, R. Longo et Y Öhrn. Mod. Phys. 66 (1994), 917.

[17] P. A. M. Dirac. The Principles of Quantum Mechanics, 3e éd. Oxford Univer-

sity Press, Oxford, 1947.

[18] P. Ehrenfest. Z. Phys. 45 (1927), 455.

[19] S. Klein, M. J. Bearpark, B. R. Smith, M. A. Robb, M. Olivucci et

F. Bernardi. Chem. Phys. Lett. 292 (1998), 259.

[20] M. Garavelli, F. Bernardi, M. Olivucci, M. J. Bearpark, S. Klein et

M. A. Robb. J. Phys. Chem. A 105 (2001), 11496.

[21] I. S. Y. Wang et M. Karplus. J. Am. Chem. Soc. 95 (1973), 8160.

[22] A. Warshel et M. Karplus. Chem. Phys. Lett. 32 (1975), 11.

[23] C. Leforestier. J. Chem. Phys. 68 (1978), 4406.

[24] M. J. Field. J. Phys. Chem. 95 (1991), 5104.

[25] B. Hartke et E. A. Carter. Chem. Phys. Lett. 189, 45 (1992), 358.

[26] B. Hartke et E. A. Carter. J. Chem. Phys. 97 (1992), 6569.

[27] B. Hartke, D. A. Gibson et E. A. Carter. Int. J. Quant. Chem. 45 (1993),

59.

[28] D. A. Gibson et E. A. Carter. J. Phys. Chem. 97 (1993), 13429.

[29] Z. Liu, L. E. Carter et E. A. Carter. J. Phys. Chem. 13 (1995), 4355.

[30] D. A. Gibson, I. V. Ionova et E. A. Carter. Chem. Phys. Lett. 240 (1995),

261.

[31] D. A. Gibson et E. A. Carter. Mol. Phys. 89 (1996), 1265.

[32] D. A. Gibson et E. A. Carter. Chem. Phys. Lett. 271 (1997), 266.

[33] H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Da-

niels, G. E. Scuseria et M. J. Frisch. J. Chem. Phys. 114, 22 (2001),

9758.
BIBLIOGRAPHIE 223

[34] L. Verlet . Phys. Rev. 159, 1 (1967), 98.

[35] M. P. Allen et D. J. Tildesley . Computer Simulation of Liquids. Oxford
University Press, Oxford, 1989.
[36] C. W. Gear . Numerical Initial Value Problems in Ordinary Dierential Equa-
tions. Prentice Hall, Englewood Clis, 1971.
[37] J. Stoer et R. Burlisch . Introduction to Numerical Analysis. Springer, New
York, 1980.
[38] W. H. Press, B. P. Flannery, S. A. Teukolsky et W. T. Vetterlin .
Numerical Recipes. Cambridge University Press, Cambridge, 1989.
[39] J.M. Millam, V. Bakken, W. Chen, W.L. Hase et H.B. Schlegel . J.
Chem. Phys. 111 (1999), 3800.
[40] G. J. Martyna et M. E. Tuckerman . J. Chem. Phys. 102, 20 (1995), 8071.
[41] R. Fortrie . Mémoire de DEA, Université Paul Sabatier, Toulouse, 2001.
[42] D. A. McQuarrie . Statistical Mechanics. Harper Collins, New York, 1976.
[43] M. E. Tuckerman, G. J. Martyna et B. J. Berne . J. Chem. Phys. 97
(1992), 90.
[44] G. J. Martyna, M. E. Tuckerman, D. Tobias et M. Klein . Mol. Phys.
87 (1996), 1117.
[45] H. F. Trotter . Proc. Am. Math. Soc. 10 (1959), 545.
[46] J. C. Slater . Phys. Rev. 34, 10 (1929), 1293.
[47] J. C. Slater . Phys. Rev. 38 (1931), 1109.
[48] N. W. Ashcroft et N. D. Mermin . Solid State Physics. Saunders College,
Philadelphia, 1976.
[49] E. Hückel . Z. Physik 60 (1930), 423.
[50] E. Hückel . Z. Physik 70 (1931), 204.
[51] S. F. Boys . Proc. Roy. Soc. London A200 (1950), 542.
224 BIBLIOGRAPHIE

[52] G. G. Hall. Proc. Roy. Soc. London A205 (1951), 541.

[53] C. C. J. Roothan. Rev. Mod. Phys. 31 (1960), 179.

[54] P.-O. Löwdin. Phys. Rev. 97 (1955), 1474.

[55] B. O. Roos. Adv. Chem. Phys. 69 (1987), 399.

[56] C. Møller et M. S. Plesset. Phys. Rev. 46 (1934), 618.

[57] B. Huron, J.-P. Malrieu et P. Rancurel. J. Chem. Phys. 58 (1973), 5745.

[58] R. Cimiraglia. J. Chem. Phys. 83 (1985), 1746.

[59] R. Caballol et J.-P. Malrieu. Chem. Phys. Lett. 188 (1992), 543.

[60] K. Andersson, P. Å Malmqvist, B. O. Roos, A. J. Sadlej et K. Wo-

linski. J. Chem. Phys. 94 (1990), 5483.

[61] F. Coester. Nucl. Phys. 7 (1958), 421.

[62] F. Coester et H. Hümmel. Nucl. Phys. 17 (1960), 477.

[63] J. Cizek et J. Paldus. Physica Scripta 21 (1980), 251.

[64] R. J. Bartlett. J. Chem. Phys. 93 (1989), 1697.

[65] P. Hohenberg et W. Kohn. Phys. Rev. B 136 (1964), 864.

[66] W. Kohn et L. J. Sham. Phys. Rev. A 140 (1965), 1133.

[67] A. D. Becke. Phys. Rev. A 38 (1988), 3098.

[68] C. Lee, W. Yang et R. G. Parr. Phys. Rev. B 37 (1988), 785.

[69] A. D. Becke. J. Chem. Phys. 98 (1993), 5648.

[70] H. Hellman. J. Chem. Phys. 3 (1935), 61.

[71] P. Gombas. Z. Phys. 40 (1935), 117.

[72] Ph. Durand et J.-C. Barthelat. Theor. Chem. Acc. 38 (1975), 293.

[73] P. Durand et J.-P. Malrieu. Advances in Chemical Physics : ab initio me-

thods in quantum chemistry, vol. 117. Wiley, 1987, p. 321.

BIBLIOGRAPHIE 225

[74] R. Poteau, I. Ortega, F. Alary, A. Ramìrez-Solìs, J.-C. Barthelat

et J.-P. Daudey. J. Phys. Chem. A 105 (2001), 198.

[75] R. Poteau, , F. Alary, H. A. E. Makarim, J. L. Heully, J.-C. Barthe-

lat et J.-P. Daudey. J. Phys. Chem. A 105 (2001), 206.

[76] L. Maron, O. Eisenstein, F. Alary et R. Poteau. J. Phys. Chem. A 106

(2002), 1797.

[77] F. Bessac, F. Alary, R. Poteau J.-L. Heully et J.-P. Daudey. J. Phys.

Chem. A 107 (2003), 9393.

[78] F. Bessac, F. Alary, Y. Carissan, J.-L. Heully, J.-P. Daudey et R. Po-

teau. J. Mol. Struct. (Theochem) 632 (2003), 43.

[79] P. Ordejón. Comp. Math. Sci. B 12 (1998), 157.

[80] G. E. Scuseria. J. Phys. Chem. A 103 (1999), 4782.

[81] A. Warshel et L. Levitt. J. Mol. Biol. 103 (1976), 227.

[82] V. Théry, D. Rinaldi, J.-L. Rivail, B. Maigret et G. Ferenczy. J.

Comp. Chem. 15 (1994), 269.

[83] M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber et

K. Morokuma. J. Chem. Phys. 100 (1996), 19357.

[84] S. Dapprich, K. S. Komaromi, K. S. Byun et K. Morokuma. J. Mol.

Struct. (Theochem) 461 (1999), 1.

[85] F. Maseras et K. Morokuma. J. Comp. Chem. 9 (1995), 1170.

[86] S. Humbel, S. Sieber et K. Morokuma. J. Chem. Phys. 105 (1996), 1959.

[87] P. Pulay. Mol. Phys. 17 (1969), 197204.

[88] H. B. Schlegel. Theor. Chem. Acc. 103 (2000), 294296.

[89] M. Head-Gordon et J. A. Pople. J. Phys. Chem. 92 (1988), 30633069.

[90] A. Szabo et N. S. Ostlund. Modern Quantum Chemistry : Introduction to

Advanced Electronic Structure Theory, 1e , révisée éd. Dover, New York, 1996.
226 BIBLIOGRAPHIE

[91] M. E. Tuckerman et M. Parrinello. J. Chem. Phys. 101, 2 (1994), 1302.

[92] J. P. Ryckaert, G. Cicotti et H. J. Berendsen. J. Comp. Phys. 23 (1977),

327.

[93] H. C. Andersen. J. Comput. Phys. 52 (1983), 24.

[94] R. Mac Weeny. Rev. Mod. Phys. 32 (1960), 335.

[95] H. B. Schlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E.

Scuseria et M. J. Frisch. J. Chem. Phys. 117, 19 (2002), 8694.

[96] P. C. Hariharan et J. A. Pople. Theor. Chem. Acc. 28 (1973), 213.

[97] A. Bergner, M. Dolg, W. Kuechle, H. Stoll et H. Preuÿ. J. Mol.

Phys. 80 (1993), 1431.

[98] H. C. Andersen. J. Chem. Phys. 72 (1980), 2384.

[99] S. Nosé et M. Klein. Mol. Phys. 50 (1983), 1055.

[100] S. Nosé. J. Chem. Phys. 81 (1984), 511.

[101] W. G. Hoover. Phys. Rev. A 31 (1984), 1695.

[102] M. E. Tuckerman, C. J. Mundy et G. J. Martyna. Europhys. Lett. 45

(1999), 149155.

[103] G. J. Martyna, M. Klein et M. E. Tuckerman. J. Chem. Phys. 97 (1992),

2635.

[104] H. Yoshida. Phys. Lett. A 150 (1990), 262.

[105] M. Suzuki. J. Math. Phys. 32 (1991), 400.

[106] R. G. Gordon. J. Chem. Phys. 43 (1965), 1307.

[107] D. W. Noid, M. L. Koszykowski et R. A. Marcus. J. Chem. Phys. 67

(1977), 404.

[108] R. Ramírez, T. López-Ciudad, P. Kumar et D. Marx. J. Chem. Phys.

121 (2004), 39733983.

[109] M. Aida et M. Dupuis. J. Mol. Struct. (Theochem) 633 (2003), 247.

BIBLIOGRAPHIE 227

[110] K. Wolinski, J. F. Hilton et P. Pulay. J. Am. Chem. Soc. 112 (1990),

8251.

[111] J. R. Cheeseman, G. W. Trucks, T. A. Keith et M. J. Frisch. J. Chem.

Phys. 104 (1996), 14.

[112] T. Helgaker, M. Watson et N. C. Handy. J. Chem. Phys. 113 (2000),

9402.

[113] R. S. Berry. J. Chem. Phys. 32 (1960), 933.

[114] R. S. Berry. Rev. Mod. Phys. 32 (1960), 447.

[115] M. J. S. Dewar, E. G. Zoebisch, E. F. Healy et J. J. P. Stewart. J.

Am. Chem. Soc. 107 (1985), 3902.

[116] G. M. Torrie et J.-P. Valleau. Chem. Phys. Lett. 28 (1974), 578.

[117] G. M. Torrie et J.-P. Valleau. J. Comput. Phys. 23 (1977), 187.

[118] E. A. Carter, G. Ciccotti, J. T. Hynes et R. Kapral. Chem. Phys. Lett.

156 (1989), 472.

[119] E. Paci, G. Ciccotti, G. Ferrario et R. Kapral. Chem. Phys. Lett. 176

(1991), 581.

[120] L. Rosso et M. E. Tuckerman. Molec. Simul. 28 (2002), 91.

[121] L. Rosso, P. Mináry, Z. Zhu et M. E. Tuckerman. J. Chem. Phys. 116

(2002), 4389.

[122] A. Laio et M. Parrinello. Proc. Natl. Acad. Sci. 99 (2002), 12562.

[123] M. Iannuzzi, A. Laio et M. Parrinello. Phys. Rev. Lett. 90 (2003), 238302.

[124] C. Dellago, P. G. Bolhuis et P. L. Geissler. Advances in Chemical

Physics , vol. 123. Wiley, 2002, pp. 178.

[125] W. K. den Otter et W. J. Briels. J. Chem. Phys. 109 (1998), 4139.

[126] M. Sprik et G. Ciccotti. J. Chem. Phys. 109, 18 (1998), 7737.

[127] W. K. den Otter et W. J. Briels. Mol. Phys. 98 (2000), 773.

228 BIBLIOGRAPHIE

[128] W. K. den Otter. J. Chem. Phys. 112 (2000), 7283.

[129] R. Kutteh et T. P. Straatsma. Reviews in Computational Chemistry,

vol. 12. Wiley-VCH, 1998, pp. 75136.

[130] V. Barone. J. Chem. Phys. 120 (2004), 3059.

[131] E. A. Solov'ev. Sov. Phys. JETP 48 (1978), 635.

[132] J. Huang, J. J. Valentini et J. T. Muckerman. J. Chem. Phys. 102 (1994),

5695.

[133] B. R. Johnson. J. Chem. Phys. 83 (1985), 1204.

[134] R. T. Skodje, F. Borondo et W. P. Reinhardt. J. Chem. Phys. 82 (1985),

4611.

[135] Q. Sun, J. M. Bowman et B. Gazdy. J. Chem. Phys. 89 (1988), 3124.

[136] D. Maynau, S. Evangelisti, N. Guihéry, C. J. Calzado et J.-P. Mal-

rieu. J. Chem. Phys. 116 (2002), 10060.

[137] F. Bessac, S. Hoyau et D. Maynau. Communication privée .

[138] J. Pitarch-Ruiz, S. Evangelisti et D. Maynau. Int. J. Quant. Chem. 101

(2005), 325.

[139] A. Migani et M. Olivucci. Conical intersections and organic reaction me-

chanisms, vol. 15. World Scientic, 2004, pp. 271322.

[140] N. Ferré et M. Olivucci. J. Am. Chem. Soc. 125 (2003), 6868.

 2
Lassata, necdum satiata

2 Fatiguée, mais non encore rassasiée. Mots empruntés à un vers de Juvénal (VI, 130), dans la

peinture énergique qu'il trace des débordements nocturnes de Messaline. Le vers complet est : Et

lassata viris, necdum satiata recessit.

 Résumé

Ce mémoire traite de la dynamique des molécules envisageant une description classique

pour les noyaux et quantique pour la structure électronique. Les approches  Born-

Oppenheimer  et  Car-Parrinello  sont discutées et un nouvel algorithme est présenté

puis validé de par la bonne conservation de l'énergie totale au cours du temps. Il est

ensuite étendu pour simuler l'ensemble canonique et appliqué à la détermination de

caractéristiques spectroscopiques de systèmes moléculaires. L'estimation de quantités,

telle l'énergie libre, est considérée à l'aide de la théorie de l'ensemble blue-moon.

Cette méthode est appliquée à deux réactions chimiques, mettant en évidence la mise

en défaut de l'approche traditionnelle basée sur l'approximation harmonique. Enn,

l'estimation de l'énergie de point zéro au delà de cette approximation est abordée.

Mots-clefs

Dynamique moléculaire ab initio  Car-Parrinello  Spectroscopie  Réactivité

chimique  Méthodologie

Abstract

This thesis concerns the dynamics of molecules where the nuclei and electrons are

treated by classical and quantum mechanics, respectively. Both Born-Oppenheimer and

Car-Parrinello approaches are discussed, a novel algorithm is presented and validated by

considering the good conservation of the total energy. This algorithm is then extended

to simulate the canonic ensemble and employed for the determination of spectroscopic

characteristics of molecular systems. By means of blue-moon ensemble, the free

energy estimation is considered. This approach is used for two chemical reactions and

shows the failure of the usual approach which is based on the harmonic approximation.

Finally the estimation of zero point energy beyond this approximation is treated.

Keywords

Ab initio molecular dynamic  Car-Parrinello  Spectroscopy  Chemical reactivity 

Methodology