Les surfaces en metaHurgie Ann. Chim. Fr., 1986, 11, pp.

29-36

ETUDE DES SURFACES PAR LA DYNAMIQUE MOLECULAIRE

V. PONTIKIS

Section de Recherches de Metallurgie Physique (SRMP), Centre d'Etudes Nucleaires

de Saclay, 91191 Gif-sur-Yvette Cedex, France

Summary : Molecular dynamics study of surfaces. - The classical molecular dyna-

mics method allows for the study of the structure of solid surfaces at the ato-
mic scale and the computation of their thermodynamical properties. We present
the basic features and the most important Limitations of the method. To illus-
trate the interest of studying surface properties by molecular dynamics, we re-
view the results of recent work on the structural and dynamical surface proper-
ties.

Resume : La dynamique moleculaire, dans L'approximation de La mecanique classi-

que, est une rnethode qui permet L'etude de La structure des surfaces a L'echelle
atomique et le calcul de Leurs proprietes thermodynamiques. Nous presentons Le
principe de La methode et ses Limites. Pour illustrer L'interet que presente
cette methode de simulation a L'ordinateur, nous passons en revue certains re-
sultats recents concernant la structure et Les proprietes dynamiques des surfa-
ces.

1. INTRODUCTION
La simulation a L'ordinateur des phenomenes physiques constitue une appro-
che compLementaire aux etudes theoriques traditionnelles et experimentales, qui
s'est deveLoppee durant Les trente dernieres annees. La dynamique moleculaire
(OM) est une methode de simulation a l'ordinateur mise au point par Alder et
Wainwright (1) pour des systemes de spheres dures et par Rahman (2) et Verlet
(3) dans Le cas des potentiels d'interaction continus. De nombreuses etudes par
DM ant ete consacrees a l'etude des proprietes d'equilibre (statiques) des liqui-
des. De ce point de vue elle est equivaLente a la methode de Monte Carlo intro-
duite par Metropolis et al. (4). Cependant, L'etude des proprietes dynamiques
des Liquides ou des solides ne peut etre effectuee que par la DM car elle seule
permet La determination de la trajectoire du systeme etudie dans L'espace des
phases. La methode de Monte Carlo permet le calcul des proprietes thermodynami-
ques en realisant un echantillonnage de l'espace des configurations du systeme,
mais ne donne aucune information sur l'evolution temporelle de celui-ci.
Tires a
part : V. PONTIKIS, SRMP, Centre d 1 Etudes Nucleaires de Saclay,
91191 GIF-SUR-YVETTE CEDEX, FRANCE

29
 La DM conserve toute la complexite du systeme etudie, li~a la nature du
modele, des forces d'interaction et de la structure, a L'oppose de L'approche
theorique traditionnelle dont Le souci permanent est la reduction des variables
afin de pouvoir traiter le probleme d'interet. C'est pourquoi La DM peut servir
de test aux modeles theoriques et permet L'etude des consequences de Leurs hypo-
theses de depart. Les proprietes d'equilibre des solides cristallins calculees
par DM (5) ont ete comparees, dans cette optique, aux resultats de la theorie
thermodynamique de perturbation (6) et le facteur de structure dynamique S(k,w)
aux predictions de la theorie de phonons self-consistente (5). Un autre aspect
interessant de la DM est qu 1 elle permet l'etude de la structure et des mecanis-
mes de diffusion a l'echelle atomique, difficilement ou pas accessibles par
l'experience. Finalement, sous reserve que l'on dispose d'un modele satisfaisant
pour decrire Les forces interatomiques, l'etude des solides et des liquides aux
fortes pressions et densites et aux temperatures elevees peut etre realisee par
La DM.
Les surfaces et plus generalement Les defauts des solides cristallins per-
turbent fortement la densite atomique et electronique a leur voisinage et ren-
dent, de ce fait, l'etude de leur structure et de leurs proprietes difficile.
De nombreuses etudes theoriques et Les calculs de OM reposent sur L'utilisation
de potentiels de paires effectifs qui constituent une approximation simplifica-
trice des forces d'interaction reelles. Si pour L'etude de certaines proprietes
des surfaces cette approximation est mauvaise, d'autres, tels Les phenomenes de
transport, permettent de s'en contenter. Selan La demarche usuelle, Les surfaces
sont etudiees par la DM a l'aide de potentiels de paires. La comparaison des re-
sultats obtenus avec ceux issus de l'experience indique alors si cette approxi-
mation est satisfaisante ou bien si i l faut affiner Le modele d'interactions
utilise, pour obtenir un meilleur accord entre simulation et experience.
Dans ce qui suit nous decrivons brievement la methode de DM dans l'approxi-
mation de la mecanique classique et nous discutons ses principales Limites. Son
application pour l'etude des surfaces est ensuite illustree par la presentation
des resultats de travaux recents concernant respectivement : a) Les proprietes
de vibration des surfaces, b) L'evolution de La structure atomique en fonction
de la temperature, c) Les mecanismes de diffusion superficiels et la transition
commensurable-incommensurable du krypton adsorbe sur le graphite.

2. DESCRIPTION DE LA METHODE
La premiere etape d'une simulation par la DM est la definition du modele
geometrique. On calcule Les coordonnees des N atomes du modele arranges sur un
reseau compatible avec la structure choisie (cubique, hexagonale, etc •• ) et l'on
introduit eventuellement des defauts,i.e..: defauts ponctuels, dislocations OU
joints de grains. Le nombre maximum de particules, N, qu'un modele de OM peut
prendre en compte depend de La capacite memoire et de La vitesse de l'ordinateur.
Sa valeur pour Les ordinateurs Les plus performants actuellement disponibles est
N ~ 100 000 particules. Pour Le systeme Le plus grand que L'on peut done simuler,
Les effets de surface sont loin d'etre negligeables. C'est pourquoi des condi-
tions aux limites periodiques sont generalement utilisees. Lorsque la condition
de periodicite aux limites est levee suivant une direction, nous obtenons un
cristallite d'extension "infinie" delimite par deux surfaces libres normales a
celle-ci. Un tel systeme sera utilise pour la simulation par la OM de La struc-
ture et des proprietes de transport superficielles.
Dans une deuxieme etape on doit preciser la nature des interactions atomi-
ques. Generalement des potentiels effectifs de paires sont utilises,mais d'au-
tres modeles plus complexes peuvent l'etre egalement.
Finalement, Les systemes des 3N equations newtoniennes de mouvement
d2x.
10(
m.
1 d t2
- F.
1 O<
(i = 1, N et a = 1, 3) [ 1]

30
Etude des surfaces par la dynamique moleculaire
31

dote de 3N conditions initiales est integre numeriquement et permet La determi-

nation de La trajectoire du systeme etudie dans l'espace des phases. Celle-ci
determine son evolution temporelle. Par sa definition un tel systeme est isole
ce qui implique que L'energie totale est conservee.
Grace a L'hypothese ergodique, Les valeurs des grandeurs thermodynamiques
sont obtenues en calculant Les valeurs moyennes temporelles des observables mi-
croscopiques <7, 8). La temperature T, et la chaleur specifique a volume cons-
tant Cv, sont determinees a partir des expressions :
1 2
T - I: m; <v. > [2]
i=1 , N ia.
a=1,3

3Nk 2T2 3k8

B
<1 - 2C ) [3]
2
v

ou m. est la masse de la particule i, N le nombre total de particules, k La

constante de Boltzmann, v. La a-ieme composante de la vitesse de La par~icule i,
et K represente La valeur 1 ~e L'energie cinetique totale.
La pression, dont L'expression est donnee par le theoreme du viriel, Les
fonctions de correlation, Les facteurs de structure et le coefficient de diffu-
sion peuvent etre calcules a partir de la trajectoire d'equilibre du systeme (8).
La conservation de L'energie totale est le critere de convergence et de
stabilite du schema des differences finies utilise pour integrer Les equations
de mouvement. La valeur de La chaLeur specifique calculee a partir des differen-
tes portions de La trajectoire permet de savoir quand Le systeme a atteint l'e-
quilibre thermodynamique.

3. LES LIMITES DE LA METHODE

3.1. TaiLle du modele

Comme il a ete indique precedemment, Le nombre maximum de particules que
l'on peut prendre en compte dans une simulation par la DM est N ~ 100 000. Dans
la pratique Les modeles utilises contiennent rarement plus de 5 000 particules.
-+
Les proprietes du systeme dans Le domaine des petits vecteurs d'onde k ne sont
done pas accessibles par la DM. Cette limite conduit egalement a des difficultes
lorsque, par exemple, nous souhaitons etudier Les defauts ponctuels. La presence
d'une Lacune dans un modele de 10 3 a 10 5 particules correspond a une valeur de
concentration de defauts cc~ 10- 3 a 10- 5 ) qui est tres eloignee de la realite.
De telles concentrations ne sont atteintes dans Les solides cristallins au mieux
qu'au point de fusion.
Pour un modele cornportant des surfaces libres, on doit faire un compromis
entre l'aire des surfaces et le volume. Un rapport surface sur volume (S/V) trop
important, i.e. une valeur elevee de L'energie de surface par atome, affecte la
stabilite thermodynamique du systeme et fait decroitre La valeur de La tempera-
ture de fusion par rapport a celle du meme systeme sans surfaces libres (9).

3.2. Duree de La simulation

La duree rnaximale d' une traj ectoi re de DM, Liee aux performances actuelles
des ordinateurs, est de L'ordre de T ~ 10- 0 s. Par consequent il n'est pas possi-
ble d'etudier Le domaine des basses frequences v < 10 8 s. Le spectre des pho-
nons ne pourra pas etre calcule dans ce domaine et L'etude des evenements ther-
miquement actives de frequence faible n'est pas envisageable par La DM.

3.3. Le potentiel
L'etude des surfaces par la OM requiert un modele de forces d'interaction
31
satisfaisant a la fois pour le volume et au voisinage des surfaces dont la pre-
sence perturbe la ·densite atomique et electronique. De tels modeles n'existent
pas actuellement. De ce fait Les etudes de simulation font appel a des poten-
tiels de paires semi-empiriques ou empiriques dont Les parametres sont determi-
nes a partir des valeurs des constantes physiques des materiaux. Les resultats
de calculs de DM peuvent etre tres eloignes de l'experience :

a) - Le calcul conduit a une augmentation de la distance interplanaire entre

la premiere et la deuxieme couches de surface par rapport aux distances inter-
planai res en volume alors que Les experiences de diffraction d'electrons Lents
indiquent une contraction pour la plupart des metaux.

b) - On n'observe pas de reconstruction de la surface du modele de DM alors

que celle-ci est observee experimentalement dans de nombreux cas.
Des progres restent a accompLir au niveau des modeles des forces thermody-
namiques pour permettre La description des systemes reels par simulation.
Malgre ses limites importantes, la DM peut etre utilisee avec succes dans
de nombreux cas et notamment aux temperatures proches du point de fusion.

4. APPLICATIONS DE LA OM A L'ETUDE DES SURFACES

4.1. Proprietes de vibration

Aux temperatures superieures a la temperature de Debye C8 0 ), Les proprietes
de vibration des solides peuvent etre etudiees par la dynamique moleculaire.
Dans ce cas elle s'avere plus interessante que Les methodes de dynamique de re-
seau notamment parce que Les effets anharmoniques sont entierement pris en comp-
te sans aucune approximation. Dans le cas du krypton solide Les calculs de DM de
Dickey et Paskin (10) reproduisent ceux de Grindley et Howard obtenus par dyna-
mique de reseau dans L'approximation harmonique (11) et sont en accord avec
l'experience. De meme Shukla et Mountain (12) ont montre que la DM et la dynami-
que de reseau sont des methodes complementaires qui permettent une comparaison
detaiLlee des facteurs de Debye-Waller calcules avec Les valeurs experimentales.
Allen, de Wette et Rahman (13) ont calcule Les valeurs des deplacements
quadratiques moyens, <u 2 >, pour des atomes de surface et de volume dans Les so-
l ides de gaz rares. Ces auteurs ant montre que Les valeurs de <u 2 > obtenues par
DM pour Les atomes de surface sont significativement plus elevees que celles
obtenues par La dynamique de reseau, alors que pour le volume, L'accord est ex-
cellent entre Les deux methodes. Ce desaccord est attribue a l'anharmonicite,
plus forte en surface qu'en volume, et aux modifications des constantes de force
pour Les couches superficielLes dues aux relaxation et a une dilatation plus for-
te en surface qu'en volume lorsque la temperature augmente (dilatation differen-
tielle).
Tous ces effets sont pris en compte par la OM qui est done plus eff icace
que la dynamique de reseau au mains aux temperatures elevees CT> 8 0 ).
Allen et al. C13) ont egalement montre que lorsque la temperature augmente
Le rapport <u 2 f >/<u 2 > croit, de meme que le spectre des phonons s'e-
largit vers Le~u¥r~&oencesv~l~~~es. Ces resultats sont dus a l'infLuence conju-
guee de La dilatation differentielle et de L'anharmonicite.

4.2. Structure des surf aces

La structure des surfaces aux temperatures elevees a fait l'objet de plu-
sieurs travaux de simulation par La OM (14-16). Une des principales motivations
de ces etudes est d'aboutir a une meilleure comprehension des mecanismes de fu-
sion des solides cristallins. En effet, pour expliquer L'absence de toute meta-
stabilite a la fusion pour la grande majorite des materiaux, on a avance l'hy-
pothese que la fusion est un processus initie au voisinage des dislocations,

32
Etude des surfaces par la dynamique moteculaire
33

joints de grains et surfaces. Des arguments theoriques a partir de modeles sim-

ples (17-18) et des simulations par la DM (14-16) concluent en faveur d'une
instabilite des surfaces a une temperature inferieure au point de fusion.
Rosato et al. (9) et Pontikis et al. (20) ont etudie la structure de la face
<110) d'un cristal de structure cfc en fonction de la temperature par la OM. Ces
auteurs ont montre que pour des temperatures T > T = 0,7 Tf . la surface de-
vient instable du fait de la production de defautscponctuelsu~~ ~res grand nom-
0

bre. Le caractere cooperatif de cette instabilite conduit a l'identification de

cette instabilite avec une transition rugueuse (9, 20).
Cette interpretation est en accord avec celle donnee par Broughton et
Gilmer (15) et elle est compatible avec Les resultats experimentaux de Gorse et
Lapujoulade C20) obtenus pour la face (110) du cuivre par diffraction d'helium.
La 1igure 1 illustre l'aspect de la surface du modele etudie a basse temperature
et aux temperatures T > T •
Un parametre d'ordrecbien adapte a l'instabilite observee est le nombre des
liaisons coupees par atome, pour Les atomes de surface. La figure 2 montre la
variation de ce parametre avec la temperature calculee par La DM.

s
1
(a)

I
I
• I
l
I
( b) I
0.5 . I
I
I
I
I
I
(c)
I
I _!_fusionc-
cp
0.38
~
Fig.1 Fig. 2
l/"1

FIG. 1. - Instabilite d'une surface cfc (110) avant le point de fusion. Configu-
rations instantanees des deux premieres couches de surface. Potentiel Lennard-
Jones (12-6) pour l'argon. Ca) T = 40 K; Cb) T = 56 K; (c) T = 77 K (d'apres
La ref. C19)).

FIG. 2. - Variation du parametre d'ordre caracterisant l'instabilite de la sur-

face (110) en fonction de La temperature. Le point de fusion du modele est indi-
que (d'apres la ref. (19)).

4.3. Phenomenes de transport

La DM permet l'etude des mecanismes de diffusion superficielle a l'echelle
atomique et de ce fait la comparaison directe des resultats obtenus avec ceux
des experiences de microscopie a champ d'ions. Les etudes de DM ont montre l'e-
xistence d'une grande variete de mecanismes de diffusion, fortement dependants
de la structure atomique de la surface (21-23).
Le mecanisme de diffusion predominant sur Les faces <110) des solides cfc,
mis en evidence par la DM, est un mecanisme d'echange entre les adatomes et Les
33
atomes de surface. Ce mecanisme conduit a une anisotropie de diffusion faible
et independante de la temperature, malgre Le caractere fortement anisotrope de
cette surface. Wrigley et Ehrlich (24) ont mis en evidence le mecanisme d'echan-
ge, par microscopie a emission d'ions, en observant la migration d'atomes de
tungstene sur une face (110) de l'iridium. Des mesures du coefficient de diffu-
sion de l'argent sur la face (110) du cuivre, r~alis~es par Ghaleb et
Perraillon (25) conduisent a une faible anisotropie de diffusion independante de
la temperature indiquant ainsi que le mecanisme d'echange opere aussi dans ce
cas. La figure 3 est une illustration de ce mecanisme obtenue par la OM.

(a)

{ b)

( c)

FIG. 3. - Mecanismes d'echange pour la migration d'adatomes sur une face (110)
de la structure cfc. Positions atomiques instantanees, avant (a) et apres (c)
le saut, et configuration de col (b).

4.4. Phases adsorbees

De nombreux travaux experimentaux et theoriques ont ete consacres a L'etu-
de de la physisorption du krypton sur le graphite. Ce systeme constitue en effet
un modele simple permettant L'etude des proprietes de la transition entre Les
phases commensurable et incommensurable. Le krypton adsorbe sur le graphite peut
former une couche solide commensurable avec le substrat a des temperatures infe-
rieures au point de fusion. Dans ce cas un tiers seulement des sites d'adsorp-
tion est occupe. Pour des recouvrements plus forts on s'attend a l'apparition de
parois de phase incommensurable separant Les domaines commensurables. Les pre-
dictions theoriques sur La structure de ces parois sont en desaccord. Les unes
34
Etude des surfaces par la dynamique moleculaire 35

sont en faveur d'une organisation des parois selon une structure en "nid d'abeil-
le", Les autres prevoyant une structure en parois paralleles C26-27). L'expe-
rience ne permet pas malheureusement de trancher entre ces deux possibilites par
manque d'informations detaillees a l'echelle atomique.
Recemment Koch et al. (28) ont pu montrer par OM que la phase incommensura-
ble s'organise en parois formant un reseau "nid d'abeille". Cette structure res-
te stable lorsque La temperature augmente. L'epaisseur des parois est indepen-
dante de la temperature mais leur forme est diffuse aux temperatures elevees.
Ces resultats constituent un exemple illustrant la maniere dont la simula-
tion, en apportant des informations auxquelles l'experience n'a pas acces,
permet de trancher entre differentes predictions theoriques.

5. REMARQUES ET CONCLUSION
La simulation a l'ordinateur est un outil puissant pour L'etude des pheno-
menes physiques et l'elaboration de modeles, mais, neanmoins, un certain scepti-
cisme est souvent exprime a son encontre (29). Les experimentateurs doutent de
son interet, car "elle ne pourra pas remplacer l'experience" et, quant aux theo-
riciens, ils considerent que c'est une mauvaise approche pour L'etude des pheno-
menes physiques. Cette situation resulte probablement d'une meconnaissance de
L'approche et de ses potentialites et est deja en train d'evoluer rapidement.
Par certains aspects la simulation a l'ordinateur est tres proche de l'experien-
ce. Le systeme etudie est "prepare" et uobserve" pour en etudier Les proprietes.
C'est pourquoi certains designent la simulation par le terme "experience" a
l'ordinateur. La demarche suivie pour extraire l'information pertinente de la
grande quantite de donnees resultant de la simulation est quant a elle similaire
a l'approche theorique.
Finalement, le risque de faire de mauvaises ou d'inutiles simulations est
aussi grand que le risque de faire de mauvaises experiences ou d'elaborer des
theories sans interet.
IL est certain que Les developpements, presents et futurs, des ordinateurs
entraineront L'utilisation de plus en plus frequente des techniques de simula-
tion notamment en Sciences des Materiaux.

6. BIBLIOGRAPHIE
( 1) ALDER (B.J.), WAINWRIGHT (T.E.), J. Chem. Phys., 1959, 31, 459.
(2) RAHMAN (A.), Phys. Rev., 1964, 136, A104. -
(3) VERLET CL.), Phys. Rev., 1967, 159, 98.
(4) METROPOLIS (N.), ROSENBLUTH CA.W.), ROSENBLUTH (M.N.), TELLER CA.M.),
TELLER CE.), J. Chem. Phys., 1953, 21, 1087
( 5) HENSEN (J.P.), Phys. Rev., 1976, 815, 878
(6) WEIS (J.J.), Malec. Phys., 1974, 28, 187
(7) LEBOWITZ (J.L.), PERCUS (J.K.), VER°LET (L.), Phys. Rev., 1967, 153, 250
(8) HENSEN {J.P.), Mc DONALD (l.R.), Theory of simple liquids, Academic Press~
New York, 1976
(9) ROSATO CV.), CICCOTTI CG.), PONTIKIS CV.), Soumis a Phys. Rev. B, 1985
(10) DICKEY (J.M.), PASKIN (A.), Phys. Rev., 1969, 188, 1407
( 11) GRINDLEY CJ.), HOWARD CR.), Lattice Dynamics, R.F. Wallis ed., p. 129, 1965
{12) SHUKLA (R.C.), MOUNTAIN (R.D.), Phys. Rev., 1982, 825, 3649
(13) ALLEN CR.E.), DE WETTE CF.W.), RAHMAN CA.), Phys. Rev., 1969, 179, 873
(14) BROUGHTON (J.Q.), WOODCOCK (L.V.), J. Phys. C, 1979, 11, 2743
( 15) BROUGHTON (J.Q.), GILMER (G.H.), {a) J. Chem. Phys., 1983, 79, 5105. Cb) J.
Chem. Phys., 1983, 79, 5119 -
(16) GILMER (G.M.), BROUGHTON (J.Q.), J. Vac. Sci. Technol., 1983, 81, 298
(17) PIETRONERO CL.), TOSATTI (E.), Solid State Commun., 1979, 32, "2"5°5
(18) JAYANTHI CC.), TOSATTI CE.), PIETRONERO (L.), Phys. Rev., 1985, 831, 3456
(19) PONTIKIS CV.), ROSATO CV.), Surf. Sci., 1985, a paraitre.
C2C) GORSE CD.), LAPUJOULADE CJ.)> Surf. Sci., 1985, a paraitre
35
C21) DE LORENZI CG.), JACUCCI CG.), PONTIKIS CV.), Proc. IC55-4 et ECOSS-3,
Cannes, 1980, Vol. I, Ed. D.A. Degras et M. Costa, p. 54 et Surf. Sci.,
1982, 116, 391
C22) MRUZIK (M.R.), POUND (G.H.), J. Phys., 1981, F11, 1403
(23) GAROFALINI CS.M.), HALICIOGLU CT.), Surf. Sci., 1981, 104, 199
(24) WRIGLEY CJ.D.), EHRLICH CG.), Phys. Rev. Letters, 1980, 44, 661
(25) GHALEB CD.), PERRAILLON CB.), Surf. Sci., 1985, a paraitre
(26) VILLAIN CJ.), GORDON CM.B.), Surf. Sci., 1983, 125, 1
(27) COPPERSMITH (D.S.), ~ISHER (D.S.), HALPERIN (8.J.), LEE (P.A.),
BRINKMAN CW.F.), Phys. Rev. Letters, 1981, 44, 549. Phys. Rev., 1982, 825,
349 ~
(28) KOCH (S.W.), RUDGE CW.E.), ABRAHAM CF.F.), Surf. Sci., 1984, 145, 329
(29) ANDERSON CP.W.), La Recherche, 1980, 107, 98

• © Masson, Paris, 1986

36