Total energy calculation B. Civalleri

• Overview and goals
• Example 1. A simple inorganic system: MgO
• Example 1.1: Total energy calculation
• Example 1.2: A simple geometry optimization
• Example 1.3: Cohesive energy
• Example 1.4: Super cell calculation and resource usage
• Exercise 1. A simple covalent system: Si
• Exercise 2. A metallic solid: Be
• Exercise 3. A molecular crystal: Urea

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.

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. 

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.7466419192577E+02 CYCLES   6

 TOTAL ENERGY(HF)(AU)(   6) -2.7466419192577E+02 DE 4.7E-08 tst 1.5E-08 PX 6.5E-04 

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

Density functional theory (DFT) calculation

In CRYSTAL 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.7402726958443E+02 CYCLES   8

 ENERGY EXPRESSION=HARTREE+FOCK EXCH*0.00000+(LDA      EXCH)*1.00000+PZ      CORR

 TOTAL ENERGY(DFT)(AU)(  8) -2.7402726958443E+02 DE-5.1E-08 tester 3.2E-07 

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: From CRYSTAL14 the use of LEVSHIFT and FMIXING options is the default for DFT calculations. For instance, for a LDA (LDA-PZ) calculation the input lines in bold characters are implicit (see the CRYSTAL User's Manual for details):

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. Note that the purpose of this exercise is purely didactic: for the automatic geometry optimization procedure, see "Geometry optimization".

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 xmgrace or gnuplot program to plot the curve.

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

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.

where E_a 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). To check the importance of the  Basis Set Superposition Error (BSSE)  use the ATOMBSSE keyword: an atom is extracted from the bulk, surrounded by the basis functions of the neighboring atoms up to a given distance.
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 accept the keyword ATOMHF with a DFT Hamiltonian, but ATOMBSSE is allowed.

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):

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.  To locate the these two values, look at the end of the output for a line like the following one:

TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT END         TELAPSE       75.88 TCPU       75.88 

or, if you are running a Unix system, open a shell in the folder containing the output and digit the command:

grep 'TTT END' *.out

This command will return the calculation time of all outputs cointained within the folder.
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".

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.

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

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 input" for 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