One-electron properties

 

Introduction

One-electron properties and wave function analysis can be computed from the SCF wave function by running the program properties. The selection of the property to be computed is driven by keywords (last keyword END)

At the end of the SCF process, the program crystal writes information on the crystalline system and its wave function as unformatted sequential data in Fortran unit 9, and as formatted data in Fortran unit 98. The tutorial refers to data stored in binary form ( Fortran unit 9).

The data set in fortran unit 9 (or 98) is read when running properties, and cannot be modified.

The data include:

Hamiltonian eigenvectors are not stored, they must be recomputed. when necessary (keyword NEWK).

The purpose of this tutorial is to teach how to extract from the wave function some basic one-electron properties:

  1. Atoms and bond populations (Mulliken Population Analysis) (keyword PPAN)
  2. Charge Electron Density (keyword ECHG)
  3. Band Structure ( keyword BAND)
  4. Density of States (keyword DOSS - preliminary computation of eigenvectors )

To compute one electron properties of a system, a "good" wave function for that system must be computed first. Convergence should be checked by restarting SCF from the solution obtained when convergence tools have been applied (FMIXING, LEVSHIFT, SPINLOCK, etc).

One-electron properties

Mulliken population analysis

A Mulliken population analysis can be performed by inserting the keyword PPAN. For each non-equivalent atom, information on the bond populations of first six neighbors is printed.

Mulliken population analysis output block follows: 

 MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM
 ALPHA+BETA ELECTRONS
 MULLIKEN POPULATION ANALYSIS - NO. OF ELECTRONS 20.000000

 ATOM Z CHARGE A.O. POPULATION

 1 MG 12 10.021 2.000 1.999 1.995 1.995 1.995 0.006 0.010 0.010
 0.010
 2 O 8 9.979 2.007 1.176 1.146 1.146 1.146 0.821 0.846 0.846
 0.846

 ATOM Z CHARGE SHELL POPULATION

 1 MG 12 10.021 2.000 7.986 0.035
 2 O 8 9.979 2.007 4.614 3.358

 OVERLAP POPULATION CONDENSED TO ATOMS FOR FIRST 6 NEIGHBORS

 ATOM A 1 MG ATOM B CELL R(AB)/AU R(AB)/ANG OVPOP(AB)
 2 O ( -1 0 0) 3.978 2.105 -0.011
 1 MG ( -1 0 0) 5.626 2.977 0.000
 2 O ( 0 0 0) 6.890 3.646 0.000
 1 MG ( -1 -1 1) 7.956 4.210 0.000
 2 O ( 1 0 -1) 8.895 4.707 0.000
 1 MG ( -2 0 1) 9.744 5.156 0.000

 ATOM A 2 O ATOM B CELL R(AB)/AU R(AB)/ANG OVPOP(AB)
 1 MG ( 1 0 0) 3.978 2.105 -0.011
 2 O ( -1 0 0) 5.626 2.977 -0.013
 1 MG ( 0 0 0) 6.890 3.646 0.000
 2 O ( -1 -1 1) 7.956 4.210 0.000
 1 MG ( -1 0 1) 8.895 4.707 0.000
 2 O ( -2 0 1) 9.744 5.156 0.000

Total atomic charges, atomic orbital and shell populations are reported as well as the overlap populations. Atomic charges are very close to a model of MgO as Mg++ O--. The overlap population between nearest neighbors is close to zero, as expected for an ionic crystal.

Compare the HF results with different levels of theory: LDA (LDA-PZ), GGA (PWGGA) and with a hybrid (B3LYP) method.

Here are the total atomic charges obtained at different level of theory:
HF LDA GGA B3LYP
Mg 10.021 10.123 10.104 10.091
O 9.979 9.877 9.896 9.909

Note that Mulliken population analysis is an arbitrary scheme for partitioning total electron charge in atom and bond contributions. Indeed, all such schemes are ultimately arbitrary. Atomic charges, unlike the electron density, are not a quantum mechanical observable, and are not unambiguously predictable from first principles.

Electron Charge Density

The ground-state electron charge density is an observable of primary importance. It may be expected that both HF and DFT hamiltonians reproduce the essential features of the density, such as the concentration of charge along covalent bonds, the nearly uniform distribution of the conduction electrons in metals, the core expansion or contraction. In an CO-BF-LCAO scheme the electron density can be expressed as follows:

$\displaystyle \rho (\hbox{\boldmath {$r$}})$ $\textstyle =$ $\displaystyle \sum_{\hbox{\boldmath {$g$}} \hbox{\boldmath {$l$}}} \sum_{\mu \n... ...x{\boldmath {$r$}}) \chi^{\hbox{\boldmath {$l$}}}_{\nu}(\hbox{\boldmath {$r$}})$


CRYSTAL computes the charge density (and its derivatives, if requested), in a grid of points defined in input. See keyword MAPNET in CRYSTAL User's Manual. Data are written in fort.25, and saved in input_file_name.f25.

The total electron density maps provide a pictorial representation of the total electronic distribution. However, more useful information is obtained by considering difference maps, given as a difference between the crystal electron density and a "reference" electron density. The "reference" density is a superposition of atomic (or ionic) charge distributions. In this case, two identical run of the ECHG option must be submitted (in the same input file) and the second must be preceded by the keyword PATO, to obtain a density matrix as superposition of atom densities). See keyword PATO in CRYSTAL User's Manual.

The isolines in the figures represent negative (dotted-line) positive (continuous) and zero value (dot-dashed) of the density in atomic units (electron/bohr3). In the total electron density maps the isolines span the range from -0.1 to 0.1 a.u with a step of 0.01 a.u.; in the difference maps the range is reduce from -0.01 to 0.01 with a step of 0.001 a.u..

MgO charge density maps.

ECHG keyword
0 order of the derivatives: if not 0, the charge density gradients are computed
65 number of point along the B-A segment (see CRYSTAL User's Manual)
COORDINA cartesian coordinates of points A.B,C defining the window in a plane
-4. $\mbox{ }$-4. $\mbox{ }$ 0.0 cartesian coordinates of point A (Angstrom)
4. $\mbox{ }$-4. $\mbox{ }$ 0.0 cartesian coordinates of point B (Angstrom)
4. $\mbox{ }$ 4. $\mbox{ }$ 0.0 cartesian coordinates of point C (Angstrom)
MARGINS margins are added to the window, in the order AB,CD,AD,BC
1.5 1.5 1.5 1.5 width of the margins, in the order AB,CD,AD,BC
END

Charge densities can be also calculated directly from XCrySDen see here

Exercises -

  1. Compute MgO wave function (input mgo.d12)
  2. Run properties with the given input file to obtain a total electron density map (see below).
  3. Add the definition of a new density matrix (using the keyword PATO) to the input, as superposition of isolated atoms (or ions) densities (according to the electron charge you attributed to the atoms when defining the basis set). Repeat the ECHG block.
    To plot the difference of the first two blocks of data in mgo.25 (first block, crystal density; second block, sum of atom/ion densities) edit the file mgo.maps:

    1 line 160 1 1 0

    Explain how the electronic structure is affected by the crystal field; which kind of bond is realized in the crystal?

Figure 2: MgO electron density: total and difference maps.
\begin{figure}{\psfig {figure=tabelle/mgodiff.eps,height=16cm}} \protect\end{figure}

 

Band Structure

Useful information can be extracted from single particle hamiltonian (HF/DFT) eigenvalues spectrum, band structure, resulting from SCF calculation. The topology of the occupied manifold and of the first conduction bands is usually correct. The shape of the Fermi surface in metals can therefore be predicted with fair accuracy.
In 2D systems the occurrence of surface states in the band gaps can be detected, as well as their dependence on surface relaxation or reconstruction; the presence of adsorbed species corresponds to new peaks in the DOSS energy spectra.
Furthermore, by considering the manifold of occupied states and by analyzing the composition of the corresponding crystalline orbitals, one can get information on the nature of the crystalline chemical bonds.

In the input file for the band structure (see also the CRYSTAL User's Manual, keyword BAND) two main information must be specified:

  1. the path in the Brillouin zone of the reciprocal space.
    The shape of the Brillouin zone depends on the reciprocal lattice, which is univocally defined by the direct lattice (that is the Bravais lattice of the system). A selection of convenient paths in the reciprocal space are reported for each Bravais lattice in the sheet documentation;
  2. how many bands have to be written in fort.25 for plotting. Suggestion: all. Computing bands can be very demanding for large system, avoid missing information.
  3. the range of energy/bands to be visualized. Sometimes it is important to concentrate the attention in the region close to the Fermi energy rather than display the whole set of bands.


MgO BAND structure

The first step is to define the number and position of the system bands. There are two atoms per unit cell, a Mg and an O, described by means of a 8-61G and 8-51G basis set, respectively, for a total of 18 atomic orbitals. Consequently, the number of bands is 18.

The second step consists in the designation of a proper path in the Brillouin zone. MgO Bravais lattice is face-centered cubic (see the sheet documentation for the Brillouin zone of a face-centered cubic lattice).

The third step is the selection of the bands to be plotted. The number of electrons per cell is 20 and the core electrons are 12. This means that the lower six bands, which correspond to the core atomic orbitals, are flat, low in energy and not sensitively affected by the formation of the crystalline network. To investigate the energy region involved in the bond formation we can select the band from the 7-th to the 14-th.

Here is reported and commented the input file to perform a band structure evaluation.

BAND keyword
MgO title
4 $\mbox{ }$ 8$\mbox{ }$ 60$\mbox{ }$ 1 18$\mbox{ }$ 1$\mbox{ }$ 0

 

 

 

 4: number of segments in the reciprocal space to explore;
8: shrinking factor in term of which the coordinates of the extreme of the segments are expressed;
60: total number of k points along the path;
1: first band to be saved;
18: last band to be saved;
1: plotting option (if 1, write data on fortran unit 25);
0: no printing option are activated
0 0 0 $\mbox{ }$ 4 0 4 path in the reciprocal space: $\Gamma \rightarrow$ X
4 0 4 $\mbox{ }$ 4 2 6 X $\rightarrow$ W
4 2 6 $\mbox{ }$ 4 4 4 W $\rightarrow$ L
4 4 4 $\mbox{ }$ 0 0 0 L $\rightarrow \Gamma$


Exercises:

Calculate the band structure of MgO using XCrySDen see here

 

Figure 6: MgO doss and band structure.
\begin{figure}{\psfig {figure=tabelle/mgoband2.eps,height=10cm}} \protect\end{figure}


 

Density of States

For a summary on the physical meaning of the density of state see CRYSTAL User's Manual).

 


The calculation of density of states requires the Fock/KS eigenvectors in the k points defined by Pack-Monkhorst net.
The keyword NEWK allows calculation of eigenvectors.

In the input file for the density of states ( DOSS) calculation (see also the CRYSTAL User's Manual) the following information must be specified:

In case of a projected density onto a subset of AOs, some other data must be specified in the input file:


MgO density of states

As in the case of the band structures, only the valence bands (from the 7-th) are included. The density of states is projected onto the whole set of AOs of the atom 1 (which is the Mg) and of the atom 2 (which is the O).

NEWK keyword
6$\mbox{ }$ 6$\mbox{ }$ 6: shrinking factor for reciprocal space Pack-Monkhorst net (see SCF input)
6: shrinking factor for reciprocal space Gilat net (see SCF input)
1 0 1: evaluation of the Fermi level with the new k-points net
0: no print options
DOSS keyword
2$\mbox{ }$ 200$\mbox{ }$ 7$\mbox{ }$ 14$\mbox{ }$ 1$\mbox{ }$ 12$\mbox{ }$ 0$\mbox{ }$

2: number of projections (the total DOS is always performed);
200: number of points along the energy axis in which the DOSS is calculated;
7: first band
14: last band
1: plot option (if 1, the program stores the data in fortran unit 25);
12: degree of the polynomial used for the DOSS expansion;
0: printing option
-1$\mbox{ }$ 1 projection onto all the AOs (-1) of Mg (1, first atom)
-1$\mbox{ }$ 2 projection onto all the AOs (-1) of O (2, first atom)


Exercises: Calculate the DOS of MgO using XCrySDen see here