    Total energy calculation B. Civalleri

 Overview and goals

The simplest type of calculation in CRYSTAL is a single point energy calculation. It consists in the calculation of the wave function and energy for a given crystalline system with a well-specified geometric structure, i.e. at a single, fixed point on the potential energy surface. The computed energy is the total energy, sum of the electronic energy and nuclear repulsion energy.
To perform a single point energy calculation a well-defined level of calculation must be specified. A level of calculation is uniquely defined by the combination of a theoretical method with a basis set.
The analysis of the computed wave function, an important step in predicting properties, is presented in the tutorial "One electron properties".

Many aspects of the crystalline behavior can be understood when the total energy of the system is known, or more often, when the energy difference between two or more electronic or nuclear configurations is known as a function of some parameters. Other important properties are related to the derivatives of the total energy. For instance, first derivatives with respect to the nuclear atomic coordinates give the forces acting on the atoms and allow to obtain the equilibrium geometry when a geometry optimization is performed. Other important properties are related to second derivatives of the total energy of the system with respect to, for instance: volume, parameters of the unit cell, and the fractional coordinates of the atoms. The related observables are the bulk modulus, the elastic constants and vibrational properties, respectively.

Some of the energy differences of interest are listed below:

 Computed energy System 1 System2 Example Cohesive Bulk Isolated atoms Ionic, covalent Interaction Bulk Isolated molecules Molecular crystals Stability Bulk Bulk Polymorphism Surface Bulk Slab MgO(100) Surface stability Slab Slab MgO(100) vs MgO(110) Adsorption Slab + molecule Slab, molecule CO on MgO(100) Substitution Defective solid Perfect system, defect C in Si Super exchange AFM bulk FM bulk NiO Solid state reaction Bulk1,Blk2 Bulk12 MgO + Al2O3 MgAl2O4

We will focus mainly on the total energy calculation and a few related energy differences for  test systems: MgO, Si, Be, and Urea. They are simple examples of the most important categories of solid compounds: ionic, covalent, metallic and molecular, respectively.

For a list of references on CRYSTAL applications to solid-state systems and processes see the CRYSTAL web site: http://www.crystal.unito.it

This tutorial consists of a worked example and three exercises. The example concerns the calculation of the energy of formation and a few other related properties of an ionic crystal, MgO. The same strategy as for the example holds for the proposed exercises, in which a covalent system, Silicon, a metallic solid, Be, and a molecular crystal, urea, are considered. Since example and exercises cover different kinds of chemical bond a direct comparison among the computed properties can be performed showing how CRYSTAL can be a powerful tool in computational solid-state chemistry.

 Example 1. A simple inorganic system: MgO

MgO has been selected as a case system for ionic compounds.

Technical notes
Input decks and some output files for the proposed example and exercises can be found in mgo.zip.

Example 1.1: Total energy calculation

Hartree-Fock (HF) calculation

The total energy of the system, computed at HF level, is given by these lines of the output:

 ``` === SCF ENDED -CONVERGENCE ON ENERGY E(AU) -2.7466415145499E+02 CYCLES 10 TOTAL ENERGY(HF)(AU)( 10) -2.7466415145499E+02 DE-2.6E-06 DP 3.0E-04 PX 3.2E-03 ```

Total energy is in hartree. The number of SCF cycles as well as SCF convergence criteria are printed

Density functional theory (DFT) calculation

In CRYSTAL09, the total energy of the system, at DFT level, is computed at each cycle of the SCF iteration. At the end, the DFT total energy is given by these lines of the output:

 ``` === SCF ENDED -CONVERGENCE ON ENERGY E(AU) -2.7402721224954E+02 CYCLES 7 ENERGY EXPRESSION=HARTREE+FOCK EXCH*0.00000+(LDA EXCH)*1.00000+PZ CORR TOTAL ENERGY(DFT)(AU)( 7) -2.7402721224954E+02 DE-4.5E-07 DP 1.5E-04 ```

The value is in hartree. The adopted exchange-correlation functional is also reported.

Run the MgO calculation at both HF and DFT level of theory and locate the lines above in the output.
Note: For DFT calculation use the LEVSHIFT and FMIXING options. For instance, for a LDA (LDA-PZ) calculation the DFT input is (see the CRYSTAL09 User's Manual for details):

 DFT DFT input block EXCHANGE LDA CORRELAT PZ ENDDFT Keyword to define the exchange functional Selected exchange functional Keyword to define the correlation functional Selected correlation functional End of the DFT input block END End of method input section LEVSHIFT 6 1 FMIXING 30 Keyword to apply an energy shift Energy shift (6*0.1 Hartrees), locking active Keyword to mix the Fock/KS matrix percent of Fock/KS matrix mixing END End of SCF input section

Example 1.2: A simple geometry optimisation

Since MgO has a cubic cell and all the atoms are in special positions, the only degree of freedom is the lattice parameter. Then, a simple geometry optimization can be devised.

Draw the total energy versus the lattice parameter (x-y chart) and locate the equilibrium lattice parameter. By starting from the experimental geometry, repeat it at different levels of theory: HF, LDA (LDA-PZ), GGA (PWGGA) and with a hybrid (B3LYP) method.

Hints: Span the lattice parameter in the region: 4.12 - 4.26 (15 pts) at HF, GGA and B3LYP level, and 4.07 - 4.21 (15 pts) at LDA level. Use the xmgr (or gnuplot) program to plot the curve.

Here are the equilibrium lattice parameters for the adopted level of theory:

 MgO HF LDA GGA B3LYP Exp. a 4.191 4.139 4.208 4.205 4.195

HF gives a lattice parameter in good agreement with the experimental value. DFT results show that LDA gives a lattice parameter which is underestimated with respect to the experimental data, whereas gradient-corrected and hybrid methods lattice parameters are in better agreement although slightly overestimated.

Example 1.3: Cohesive energy

The cohesive energy is defined as: where Ea is the atomic energy for each atom (or ion) belonging to the crystal unit cell and N is the number of atoms in the unit cell.
The choice of an atomic reference energy is a fundamental and controversial point in the evaluation of the cohesive energy.
As the basis set is not complete, two different basis sets are used for isolated atoms and for bulk. The same basis set is used for inner electrons, whereas extra functions must be added to the bulk basis set for an accurate description of the valence electron of the isolated atoms.

Three calculations are necessary to compute cohesive energy of  MgO: MgO bulk, oxygen atom, and magnesium atom.
At HF level, calculate the atomic energies of isolated atoms by crystal (ATOMHF keyword: atomic wave function only is computed, no periodic calculation). Test  one more sp function for better describing the atoms. Repeat this calculation with a 85-11G basis set for Mg and a 8-411G basis set for O. The required basis sets can be found at the CRYSTAL web site: http://www.crystal.unito.it/Basis_Sets/ptable.html

CRYSTAL does not allow calculation of atomic wave function with DFT Hamiltonian

Example 1.4: Super cell calculation and resource usage

The total energy per formula unit is very stable when evaluated with unit cells of increasing size. This aspect is very important for the evaluation of defect formation energies with the super cell scheme and when the relative stability of different magnetic phases is computed.

a) Modify the MgO input deck to create super cells containing 8, 16 and 32 atoms. Run each calculation at HF level and then compute the total energy per formula unit. Compare the results with the MgO bulk calculation. Hint: To speed up the calculation for the 32 atoms super cell reduce the value of the shrinking factor from 8 to 4.

Here is the HF total energy for the MgO super cell as well as the total energy per formula unit (hartree):

 2 8 16 32 Total energy (TE) -274.664 -1098.66 -2197.31 -4394.63 TE per formula unit -274.664 -274.664 -274.664 -274.664

Exercise: Repeat such calculations at LDA level. Is the energy stable also for DFT methods?

b) Have a look at the CPU time and at the elapsed time. The total CPU time should be separated into two contributions: on one hand, the CPU time to compute the integrals, on the other hand, the CPU time to perform the SCF iteration.
To locate the two contributions, digit the command:

and:

grep 'TTT END' *.out

The former is the CPU time for the integrals calculation, the latter is the total CPU time. To obtain the SCF CPU time, then, subtract the CPU time for the integrals calculation to the total CPU time.
How does the cost of the calculation increase with the size of the system?

Extensions

Many extensions of the proposed example are possible in which the knowledge of the total energy is a fundamental step.

For instance, in the field of surface science, different surfaces (slab models) can be cut out from the bulk (e.g. MgO(100) and MgO(110)). By comparing the total energy of bulk and slab models the surface energy can be easily computed. When the surface energy is available the different stability of surfaces can be obtained as well as the dependence of the surface energy on the thickness of the slab model. See "Surfaces and adsorption"

Another important application in materials science is the study of defective systems, in which a defect is present in the bulk structure. See "Defects".

 Exercise 1. A simple covalent system: Si

Silicon  is an important material which finds many application as semi-conductor in electronics.

Repeat the calculation previously proposed for MgO and compare the results with the MgO ones. Results can be interpreted in  the same way.

Technical notes
Input decks and some output files for the proposed exercise can be found in si.zip.

 Exercise 2. A metallic solid: Be

Be is a simple metallic solid. The Fermi energy must be computed by integration: the value of the Gilat shrinking factor, ISP (see "SCF") must be set to 2*IS (shrinking factor Pack-Monkhorst net).

Repeat the calculation proposed for MgO and Si and compare results with the previous ones. Note that Be bulk structure has two lattice parameters to be optimized. In this case, the geometry optimization of the cell may be performed by fixing the c/a ratio at its experimental value (1.5677) and by plotting the total energy vs the volume of the cell.

Technical notes
Input decks and some output files for the proposed exercise can be found in be.zip

 Exercise 3. A molecular crystal: Urea

The formation energy of a molecular crystal is obtained as difference between the energy of the bulk and the sum of the energies of isolated molecules.

Two calculations are required: urea crystal and urea molecule (there are two equivalent molecules in the crystallographic cell).

By starting at HF/3-21G level and adopting the experimental structure determined at 12 K (see  "Geometry inputfor the structural data):

1) Compute the total energy for the urea crystal;
2) For the molecular calculation, use the keyword MOLECULE to extract a molecule from the molecular crystal and run CRYSTAL

Then, calculate the formation energy of the urea crystal.

Check the computed formation energy for the basis set dependence. Repeat the HF calculation with a 6-21G(d,p), a 6-31G(d,p), and a 6-311G(d,p) basis set.

Compare the HF results with DFT calculations. Repeat the calculations above by using the following exchange-correlation functionals: SVWN, PWGGA, B3LYP.

To check the importance of the  Basis Set Superposition Error (BSSE)  use the MOLEBSSE keyword: a molecule is extracted from the bulk, surrounded by the basis functions of the neighboring molecule, up to a given distance from the outermost atoms.

Technical notes
Input decks and some output files for the proposed exercise can be found in urea.zip   