Make your own free website on

A short description of my Ph.D. thesis:


"Atomic scale simulation of the surface processes. Hydrogen interaction with the magnesium oxide (001) surface"

(Download WORD 6.0 copy of this page)

 The quick development of the computational techniques (computational power and methods), together with the progress of the visualization (graphics, 2d and 3d representation, animation) and of the strategies to handle and process the data create new opportunities in the development of nature phenomena simulation. Thus, the numerical simulation replaces the production and testing of physical prototypes, which involves a reduction of the study time and an important saving of time, materials, even of staff, facilitating a faster market launching of the products.

One of the fields that take advantage from this opportunity is material science and engineering, where the management principles involve both quantitative measurement and phenomenological behaviours. This joined information requires combinatorial experience and a good understanding of the theory, that make possible to extract, from a great amount of data, the information to be used in material projecting.

The modern quantum theory is the theoretical frame that allows the characterization of a particle system. Quantum mechanics ensure a theoretical frame (values, mathematical equations) for the characterization and analysis of the behaviour of some material systems. It is impossible to solve the involved mathematical equations in absence of a very powerful ``brain''. The present levels of soft and hard development make the quantum mechanics equations and the analytical and numerical methods of solving them become an instrument (numerical simulation of matter) to predict the properties of a material system.

Almost all the physical and chemical properties of a material system are related to the total energy or to differences of total energy. Therefore it is possible to use the capacity of the quantum simulation to predict the total energy of a electrons-nuclei system (chemical system) in order to exploit a wide range of phenomena. For example, the lattice equilibrium constant of a crystal system is the lattice constant for which the total energy is minimized and the solid surfaces and defects adopt structures that minimize their total energy. Some total energy techniques can be successfully applied in the accurate prediction of the lattice equilibrium constants, the volume modules, the phononic spectra, the piezoelectric constants and the phase tranzition, function of pressure and temperature. This instrument is useful in the exploitation of a wide range of researching fields (chemistry, material sciences, biology, geophysics).

There are many examples on the relation between the atomic and the macroscopic level. A example is the application of the quantum simulation in the material sciences, namely in the defect description of a material. The defect comes from the atomic scale, where a bond is disturbed, which leads to the character modification or even to its breaking off. Therefore, it is clear that the site and the way that a defect appears in a material can be determined by the help of quantum simulation. Perhaps one of the most important prediction in the quantum calculation is the prediction of the existence of fulerenes before making it evident by experiments.

Having in mind these possibilities, I tried to approach by numerical simulations the MgO bulk and (001) surface phenomena, the system that are very important in catalysis and in material sciences. The doctor's degree paper resulted from my work comprises 196 pages and is structured on four chapters.

In the first chapter (Crystal Surface) I presented the possibilities of approaching, at an atomic scale, the surface local or collective phenomena, with a stress on the procedures to be followed in such an analysis. Thus, there are first described the surface characteristic values and their origin in the volume values, then the possibilities to simulate the clean surface (jellium, cluster, slab and supercell methods) and the phenomena related to molecular adsorption on the crystal surfaces for small coverages (adsorption, dissociation, absorption, diffusion, recombination, desorption) or for large coverages (formation of adsorption phase domains, nature of the domain walls, the phenomenon of hysteresis, phase transitions).

Following the introduction chapter, in the second chapter (Volume properties of MgO) the volume properties of the magnesium oxide crystal are studied. In this study I used a solid physics programme, CRYSTAL 95, based on the SCF-LCAO-CO method, for which HFT and DFT hamiltonians can be used. The first part of this chapter is dedicated to the choice of the combination between the basis set (of ten basis sets of different qualities), the hamiltonian (HFT exchange potential including correction of post-SCF correlation and local DFT or including-GGA-gradient-corrections exchange and correlation potentials) and the parameters of calculation tolerance for the integrals and series in the direct (TOLINTEG) or the reciprocal space (RECIPR), as well as the establishment of the calculation schemes. See this influences of integration/sumation parameters in Figure 1.
The selection is made on the basis of accurately reproducing the wave function and the dependence of the total energy E of the elementary cell on the lattice constant value a within the range 4.23 - 4.28 Å. The accuracy of reproducing the wave function is characterized by the approaching of the virial coefficient to 1, and the accuracy of the E(a) dependence by positioning its minimum a0 as close as the experimental equilibrium value a0exp=4.210 (assessed by means of the relative error value RE=(a0exp-a)/a0exp ) and E(a0) value and by the curve form defined by the volume modulus B% % EMBED Equation.2 and its derivative against pressure B\% % EMBED Equation.2 . Following the parabolic fitting E(a) I estimated RE and computational effort for each set of basis (Figure 2). I reduced RE by controling and adjusting of the tolerance parameters for evaluation of integrals and infinite series. I identified two basis sets that involve a reduced computational effort and ER:
-All Electrons sets B5 (Mg: 86-1G, O: 8-51G)
    -for the HFT hamiltonian: ER=-0.044% (for TOLINTEG:           ,  MULTIPOLE=  , RECIPR: )
    -for the post HF-SCF correction with LYP correlation potential in a hybrid scheme: RE=-0.020% (for TOLINTEG:           ,  MULTIPOLE=  , RECIPR: )
-Valence Electrons set B9 with effective core potentials of Durand type (Mg: D-1G, O: D-11G) and DFT calculation with correlation and exchange potentials: ER=-0.080% (for TOLINTEG:           ,  MULTIPOLE=  , RECIPR: ).

The possibility to transfer the HFT-optimized effective core potentials in DFT is also analyzed.

Once the optimum work conditions identified, in the second part of the chapter the electronic structure (electronic bands - Figure 3, total and projected state densities - Figure 4, Mulliken population) and the structural properties (the lattice constant, the volume modulus, the derivative of the volume modulus against pressure and system stability) are characterized. A stress is laid on identifying the electronic bands and studying the influence of the electronic correlation. The main result is the reducing of the band gap in agreement with the experimental value by lowering of the conduction band.

In the third part, the parabolic fitting and the fitting against the state Murnaghan equation are studied comparatively, perfect correspondences being found. By solving the system made up of the two energy dependencies on the volume of the primitive cell, the volume dependence on the pressure outside the crystal is deduced (Figure 5).

The third chapter (Properties of MgO (001) surface) analyzes the properties of the ideal MgO (001) surface, the phenomena of relaxation and rumpling of MgO (001) surface, as well as their influence upon the electronic properties of the analyzed surface. At the end of this chapter, I approached the Li-doping of MgO (001) surface by replacing one cation in the first layer following the direct collision or the occupation of a cation vacancy.
In the section dedicated to the possibilities of modeling the MgO (001) surface two methods of modelling are discussed, semiempirical SCF-LCAO-MO methods (Tight Binding -BIGICON and MNDO-MOPAC 6.0) and ab initio SC-LCAO-CO methods (CRYSTAL 95). By SLAB (CRYSTAL 95) calculations it is established that for modelling the MgO (001) surface a number of three layers is enough to obtain a correct description of the electronic charge distribution, while 12-15 layers are enough for the convergence of the energy contribution to the total energy of each layer of the surface in the SLAB structure (Figure 6) or for the convergence of the upper boundary of the valence band (Figure 7). These convergences are analyzed for two AE bases (B5: 86-1G/8-51G and B7: 85-11G/8-411G). After this standardization of the SLAB method I studied formation of the surface bands and the surface share in the upper side of the valence band (Figures 8 and 9). Relaxation (-0.009 Å) and rumpling (0.021 Å) confirm the theoretical and experimental predictions: the cations of the first layer move inside the surface (0.015 Å), and the anions move outside the surface (0.0195 Å). Relaxation and rumpling alter slightly the surface bands and very slightly the charges on the first two atomic layers (Figure 10).

The study of Li-doping of the surface is preceded by the choice of the surface supercell parameters, so that the direct and indirect interaction between two defects from neighbouring supercells should be as weak as possible. It is established that a 2x2 supercell is big enough for the direct interaction to be annulled and the indirect one to contribute with less than 10-5 hatree to the total energy. I discussed MgO(001) surface doping with lithium an various reaction channels relating to the ionization level of the lithium atom and the substituted Mg atom and to the level of electric charge of the surface. Two substitution mechanisms are considered: the substitution of a magnesium ion by direct collision and the occupation of a magnesium vacancy. It was found that energetically favored is the occupation with a lithium atom of the vacancy V0 formed by removing a neutral Mg0, so that the surface remains electrically neutral. Since the heights of the potential barriers were not established, it was not possible to predict the most probable doping channel. For each type of defect the modifications of the state densities and migration of the electrons in the first two surface layers are separately studied by Mulliken analysis of the electronic charges of the ions (Figure 11).

In the last chapter (Hydrogen interaction with MgO (001) surface) I analyze the possibilities of (cluster, cluster+PC, supercell and embed) modelling the atomic and molecular hydrogen adsorption on the clean or Li-doped MgO (001) surface, as well as the molecular hydrogen dissociation on these surfaces.

Hydrogen adsorption and volume penetration are analyzed by Tight Binding method in two adsorption sites, atop a cation (Figure 12) and an anion (Figure 13). The surface was modeled by a 5x5x5 cluster made up of five atomic layers, with 13 atoms of one type and 12 atoms of the other type on each layer, alternatively arranged in a 5x5 lattice, depending on what type of atom the adsorption occurs. The adsorption height and the absorption depth under the first layer are analyzed in each type of study. It is thus established that in the phase of adsorption hydrogen is fixed at a height of 2.041 Å from the surface and 0.812 Å in the adsorption on Mg and oxygen, respectively, and in the phase of absorption hydrogen is linked to oxygen at a bond distance of 0.794Å and 0712Å, respectively, corresponding to the two types of adsorbtion. The charge migration in the system, during adsorption and absorption, is also studied. It is thus established that although the coloumbian field is not well reproduced, the system charges obtained from the Tight Binding calculation are close to the ionic values (±1.500 electrons) and that during adsorption the charge is perturbed only on hydrogen, adsorption centre and its first neighbours, while during absorption the perturbation extends to the seconds neighbours too.

The H2 adsorption atop Mg+2 was also studied by the free cluster method, representing the MgO (001) surface by a 5x5x2 free cluster to establish if the chosen (9,1) cluster permits the relaxation of the charge distribution within the cluster during adsorption. The Tight Binding results were confirmed, i.e. the electron exchange between surface and adsorbate is made mainly by the contribution of the first neighbours of the target- Mg.

Next I discuss the methods of improving the reproduction of the surface behaviour in presence of the adsorbate (cluster +PC, embed and supercell methods), starting from the necessity of an adequate description of DOS of the sensitive surface to the interaction with the atomic and molecular hydrogen.

The programme GAMESS used in cluster+PC calculations does not allow the introduction of the long-distance coloumbian field produced by the point charges ±2e that surround the cluster. Therefore, I imagined a trick, starting from the fact that the programme permits the definition of the core ECP parameters. The +2 charge was introduced as follows. The core Mg (Z=12) is described by a ECP of 10 electrons (thus removing 10 electrons and 10 protons from the core), with null effects that impose null values to all the parameters of the potential. Thus, Mg is represented by 2 electrons and two protons. By eliminating the atomic basis, the two electrons are delocalized from the atom in the entire system. The charge -2 is built up in a similar way, describing the oxygen by a null ECP potential of 10 electrons and removing the atomic basis. Therefore oxygen will have a -2 charge (a lack of two protons) and a lack of two electrons which will be annulled by the two excess electrons of a ±2 point charge. The adsorption region of the surface is modeled by a (9,1) cluster in the adsorption atop a cation (Mg+2 or Li+1) or by a cubic (4,4) cluster in the centre adsorption, embedded in a point charge lattice 9x9x4 and 10x10x4, respectively.

The SLAB-SUPERCELL modelling (programme CRYSTAL 95) is such standardized as the direct and indirect interactions of two hydrogen molecules in two neighbouring supercells are eliminated. As far as it was possible for the power of the available computer, it is established that a 2x2 supercell introduces an error of 10-6 hartree regarding the direct interaction H2-H2 and 10-5 hartree regarding the indirect interaction, mediated by the surface. In these calculations I used the acceptance parameters 6 6 6 7 14/8 8 8 and the basis B5. The SLAB structure is made up of 5 atomic layers.

The EMBED CLUSTER modelling (programme EMBED 96) is used only for the adsorption atop a cation (Mg+2 and Li+1). In this study I used the same (9,1) cluster. Preliminary calculations are performed to determine the electrochemical potential of the free cluster (0.5964 hartree/electron) and the charge transfer from the remainder of the surface to the free cluster (0.132 electrons). The lattice ionic charges of the free surface determined in the SLAB-SUPERCELL study are correctly reproduced, both for clean and Li-doped surfaces.

By comparing the DOS obtained with the four methods (free cluster, cluster+PC, embed cluster and supercell) for the free or Li-doped surface it is found that the DOS determined by the supercell method is properly reproduced by the DOS determined by the embed cluster method. It is also found that the positioning of the cluster in PC adjusts the DOS of the free cluster especially in the upper region of the VB and in the lower region of the CB, active energetic regions in H2 adsorption (Figures 14 and 15). The problems that remain are with anomal polarization of anion by +2 PC neighbors.

Then I analyzed, by the cluster+PC, the adsorption and the possibilities of dissociation of the hydrogen molecule in the adsorption sites atop the cation (Mg+2 or Li+1) or center, by building 4-dimensional PES E=E(d,h,theta,phi) function of the distance d between the hydrogen atoms, the height h of the mass centre of H2 molecule to the surface and the angles theta, phi that describe the orientation of the molecular axis, by minimizing the total energy depending on theta and phi.

For each type of adsorption I identified the molecule orientations and the energetically extreme points (d,h), as well as the dissociation possibilities of the molecule (Figures 16 17,  18). Thus, it is found that the molecule physically absorbs in 2-3 heights, characterized by the molecule orientation perpendicularly and then in parallels to the surface along the O-O bonds, with distances d of the order 0.9 Å. The molecule dissociation is prevented by potential barriers with a height of 1.20 eV (center adsorption) , 4.62 eV (adsorption on Mg+2), 2.014 eV (adsorption on Li+1). The centre dissociations and those on Mg+2 are prevented also because of the positioning of these potential barriers in the exit channel at large distances d (1.535 A and 2.140 A) in comparison with the adsorption on Li+1 (1.178 A). In the latter case, a translation excitation of the hydrogen molecule may be favourable to the direct dissociation, without remaining too long a time in the last adsorption state. The molecule may transfer its normal translation energy to the H-H vibration mode which is parallel to the surface. In the other cases the reduction in the height h leaves unchanged the distance d till the transitional state is attained. In this case the vibration excitation is possible only by accumulation of vibration energy on the H-H bond. This accumulation is facilitated by coupling the molecular vibration with the surface phonons.

Once the molecule dissociated, the hydrogen atoms occupy various sites in relation to the surface ions, but in all the situations the tendency of the hydrogens to form strong bonds with the surface anions is constant.

For the atop adsorptions the calculations in the extreme sites of PES are performed again by the other complementary methods, embed cluster and supercell. The computational effort involved by the embed method did not permit the study of all the dissociation adsorption sites. Both methods reproduce the tendency that a polarization appear in molecule in the adsorption phase through a charge transfer on the atom that is closer to the surface (0.009 for the embed calculations and 0.025 in the case of the cluster), with no charge transfer between the molecule and the surface. For both methods the energetic order of the extreme points in the phase of adsorption is maintained, but an exact characterization of the extremes implies new optimizations, that are not presented in this thesis because of the great computational effort involved. It was found that the CLUSTER+Point Charges modelling is qualitative and it does not reproduce the quantitative the charge transfer between surface and hydrogen. The set of basis STO-6G used in cluster study involves a large overlap population and clearly the Mulliken population have limited meaning. In the SLAB-SUPERCELL (program CRYSTAL 95) and EMBED CLUSTER (program EMBED 96) modellings are used sets of basis adapted at the MgO crystal (Chapter 2) in that this problem is eliminated since the diffuse GTO are eliminated. To avoid this situation an orbital projection must be used. Another reason is anomaly polarization of anion induced by +2 PC's neighbour. This problem can be eliminated by replacing of this PC's with pseudocations treated with Total Ion Potentials.


I thank to Dr. C. Pisani and co-authors and to Dr. M. Schmidt and co-autors of the CRYSTAL 95 and EMBED 96 and respectively GAMESS programs as they kindly made the programs available for me free of charge.

[ Institute of Physical Chemistry ] [ The Material Modelling Group ]
[ My Home Page ] [ My Publications ]

Developed by Viorel Chihaia

Last updated on : Mon Aug 1 11:08:41 1991