Introductory Tutorials L. Valenzano

This tutorial aims at giving a quick overview of the basic features of the CRYSTAL code from simple calculations (total energy, geometry optimization, ...) to some advanced applications in the modelling of surfaces and defective systems. The tutorials was prepared for the MSSC2007 School by gathering parts of the already existing basic tutorials.
Basic information necessary to start usage of the CRYSTAL package are also presented in this tutorial.

See:

 How to start

Basic knowledge of Unix/Linux commands (create, edit, copy, move files, definition of environment variables)  is required.

• Install and test the CRYSTAL package following the instructions. Tutorials exercises run with sequential binary (crystal)
Check the commands runcry runprop work correctly.
A free evaluation copy for Unix/Linux systems (maximum of 4 atoms per cell) can be downloaded
• Install at least one of the visualization tools described in the documentation
• Check a postscript viewer be installed in your  computer.

Background requirements:

• basic concepts of crystallography: direct lattice, reciprocal lattice, space groups;
• basic quantum chemistry: Hartree-Fock and DFT methods;
• basic solid state physics: Bloch theorem, Bloch functions, band theory, Fermi energy surface;

A review article is devoted to a basic introduction to CRYSTAL as a tool for quantum simulation in solid state chemistry :

Roberto Dovesi, Bartolomeo Civalleri, Roberto Orlando, Carla Roetti, and Victor R. Saunders
Ab Initio Quantum Simulation in Solid State Chemistry
Reviews in Computational Chemistry, Chapter 1, Volume 21
Kenny B. Lipkowitz (Editor), Raima Larter (Editor), Thomas R. Cundari (Editor)
John Wiley & Sons, Inc, New York, 2005

 CRYSTAL input overview

The input deck for wave function calculation, an ASCII text file, is read by the program crystal.

The CRYSTAL input to crystal includes a title and three sections (referred to as "blocks").

Each block consists of keywords (case insensitive, written left justified) and numerical parameters (free format).

Each block ends with the keyword END (mandatory: 3 characters only are interpreted, any ending is allowed, ENDgeom, ENDbas, etc etc) or STOP. The latter will cause immediate termination of execution.

In order to test each block, the keyword TEST can be used. Input is processed up to that block, and the program then stops.

Optional keywords can be present in each section.
Extended information on input features are in "CRYSTAL User's Manual".

The input deck has the following structure (mandatory data):

 Title input block 1 Geometry input (see tutorials: geometry editing and geometry optimization)         standard geometry input         optional geometry optimization and editing keywords END input block 2 Basis set input (see tutorial: basis set)          standard basis set input         optional basis set related keywords END input block 3 Single particle Hamiltonian (default: RHF)  and SCF control (see tutorial: Hamiltonian, computational parameters and SCF)          SHRINK         sampling in reciprocal space (for 1D-2D-3D systems only)         optional general information and SCF related keywords END

Examples of complete CRYSTAL inputs are reported in the following table:

 Molecule-0D Polymer-1D Slab-2D Crystal-3D Periodicity METHANE POLYMER SN GRAPHITE SILICON BULK Title MOLECULE POLYMER SLAB CRYSTAL Structure 0 0 0 Code for 3D - See CRYSTAL User's Manual 44 4 77 227 Space group number 4.431 2.47 5.42 Lattice parameter(s) 2 2 1 1 Number of irreducible atoms 6     0.     0.     0. 16     0.     -0.844969     0. 6     -0.333     0.333     0. 14     .125     .125     .125 Atomic number and coordinates 1     .629     .629     .629 7     .1416     .667077     -.00093 Atomic number and coordinates keyword keyword keyword keyword optional geometry editing END END END END end of geometry input --------------------------------- --------------------------------- --------------------------------- --------------------------------- --------------------------------- 6   2 16   4 6   2 14   3 basis set input 1   0   3   2.   0. 1   0   3   8.   0. 1   0   3   2.   0. 1   0   3   2.   0. 1   1   3   4.   0. 1   1   3   8.   0. 1   1   3   4.   0. 1   1   3   8.   0. 1   1 1   1   3   6.   0. 99   0 1   1   3   4.   0. 1   0   3   1.   0. 7   2 99   0 99   0 1   0   3   2.   0. 1   1   3   5.   0. 99   0 keyword keyword keyword keyword optional geometry basis set editing END END END END end of basis set input --------------------------------- --------------------------------- --------------------------------- --------------------------------- --------------------------------- SHRINK SHRINK SHRINK shrinking factor keyword (mandatory) 8 8 8 8 8 8 keyword keyword keyword keyword optional SCF editing END END END END end of SCF input

 CRYSTAL input:

## Block 1: Geometry

Basically, a crystal structure is determined by:
1) space group;
2) cell parameters;
3) coordinates  of the atoms in the asymmetric unit.

This information can be found in crystallographic databases, for instance in the ICSD database.

Searching for the MgO structure, the output from the ICSD data base is the following, where the information needed to define the crystal structure is highlighted.

 COL  ICSD Collection Code 9863  DATE Recorded Jan 1, 1980; updated Jan 19, 1999 NAME Magnesium oxide               (Title) MINR Periclase  FORM Mg O       = Mg O  TITL X-ray determination of electron-density distributions in oxides,        Mg O, Mn O, Co O, and Ni O, and atomic scattering factors of their        constituent atoms  REF  Proceedings of the Japan Academy       PJACA 55 (1979) 43-48  AUT  Sasaki S, FujinoK, TakeuchiY  SYM  x, y, z               y, z, x                z, x, y        x, z, y               y, x, z                z, y, x        x, -y, -z             y, -z, -x              z, -x, -y        x, -z, -y             y, -x, -z              z, -y, -x        -x, y, -z             -y, z, -x              -z, x, -y        -x, z, -y             -y, x, -z              -z, y, -x        -x, -y, z             -y, -z, x              -z, -x, y        -x, -z, y             -y, -x, z              -z, -y, x        -x, -y, -z            -y, -z, -x             -z, -x, -y        -x, -z, -y            -y, -x, -z             -z, -y, -x        -x, y, z              -y, z, x               -z, x, y        -x, z, y              -y, x, z               -z, y, x        x, -y, z              y, -z, x               z, -x, y        x, -z, y              y, -x, z               z, -y, x        x, y, -z              y, z, -x               z, x, -y        x, z, -y              y, x, -z               z, y, -x CELL a=4.217(1) b=4.217(1) c=4.217(1) alpha=90.0 beta=90.0 gamma=90. (Cell parameters) V=75.0 D=3.56 Z=4            (Cell Volume and Density in gr/cm3) SGR  F m -3 m      (225) - cubic     (Hermann-Mauguin symbol and Space group number) CLAS m-3m  (Hermann-Mauguin) - Oh  (Schoenflies)  PRS  cF8  ANX  AX  PARM Atom__No OxStat  Wyck -----X----- -----Y----- -----Z-----  -SOF- Mg     1  2.000    4a  0.          0.          0.   (Atom fractional coordinates)       O      1 -2.000    4b  0.5         0.5         0.5 WYCK b a  ITF  Mg  1   B=0.312  ITF  O   1   B=0.362  REM  M PDF 43-1022  RVAL 0.013

The CRYSTAL geometry input for MgO, is then

 MGO BULK - GEOMETRY TEST CRYSTAL 0 0 0 225 4.217 2 12 0. 0. 0. 8 0.5 0.5 0.5  Title of the job Dimensionality of the system (CRYSTAL => 3D) Crystallographic information (3D only) Space Group number Lattice parameters (minimal set of data) Number of atoms in asymmetric unit Atomic position specification in fractionary coordinates TEST Optional keyword to stop execution after geometry input END End of the geometry input section

### Exercise: Ionic crystal

1) Run a TEST on the geometry input of MgO as reported previously (runcry mgo_geom, if the filename of the input defined above is mgo_geom.d12)

2) Visualize the output (filename mgo_geom.out). You should get:

  ******************************************************************************* LATTICE PARAMETERS (ANGSTROMS AND DEGREES) - BOHR = 0.5291772083 ANGSTROM PRIMITIVE CELL - CENTRING CODE 5/0 VOLUME= 18.747822 - DENSITY 3.541 g/cm^3 A B C ALPHA BETA GAMMA 2.98186930 2.98186930 2.98186930 60.000000 60.000000 60.000000 ******************************************************************************* ATOMS IN THE ASYMMETRIC UNIT 2 - ATOMS IN THE UNIT CELL: 2 ATOM X/A Y/B Z/C ******************************************************************************* 1 T 12 MG 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 2 T 8 O -5.000000000000E-01 -5.000000000000E-01 -5.000000000000E-01 TRANSFORMATION MATRIX PRIMITIVE-CRYSTALLOGRAPHIC CELL -1.0000 1.0000 1.0000 1.0000 -1.0000 1.0000 1.0000 1.0000 -1.0000 ******************************************************************************* CRYSTALLOGRAPHIC CELL (VOLUME= 74.99128631) A B C ALPHA BETA GAMMA 4.21700000 4.21700000 4.21700000 90.000000 90.000000 90.000000 COORDINATES IN THE CRYSTALLOGRAPHIC CELL ATOM X/A Y/B Z/C ******************************************************************************* 1 T 12 MG 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 2 T 8 O -5.000000000000E-01 -5.000000000000E-01 -5.000000000000E-01 

3) Check the volume and the density of your system with those reported in the ICSD table above. (Note that CRYSTAL computes the density with respect to the mass of the most abundant isotope and not with respect to the atomic weight).

4) Visualize the crystal structure with moldraw.

The crystallographic or conventional cell is used as standard option in input. It may be non-primitive, which means that it may not coincide with the cell of minimum volume (primitive cell) which contains just one lattice point. Note that, for maximum calculation efficiency, CRYSTAL works on the primitive cell. Hence, the conventional cell is transformed into the primitive cell (1/4 of the conventional cell for face centred cubic lattices).

All the structural information following is referred to the primitive cell.

In the output file, the lattice parameters of the primitive cell and the corresponding atomic positions (in fractionary units) are reported. For MgO two atoms only build up the primitive cell.

In figure, the conventional cell and the primitive cell, enclosed in the conventional cell, of MgO are shown:

 CONVENTIONAL CELL PRIMITIVE CELL

### Exercise: Molecular crystal

1) Look at the following table (data taken from the CSD database) that contains the basic information you need to write the input file urea_geom.d12

 Space group P421m,   No. 113 (tetragonal) Lattice parameters a=5.565A,   c=4.684A C:   x= 0.0,   y=0.5,   z=0.3260 O:   x=0.0,   y=0.5,   z=0.5953 Atomic positions N:   x=0.1459,   y=0.6459,   z=0.1766 H(1):   x=0.2575,   y=0.7575,   z=0.2827 H(2):   x=0.1441,   y=0.6441,   z=-0.0380 Cell volume 145.06  A3 Density 1.37 g/cm3

2) In order to check the geometry of your system, modify the following input (urea_geom.d12) to run a TEST on the geometry of urea (runcry urea_geom).

 UREA BULK - GEOMETRY TEST CRYSTAL 0 0 0 113 5.565 4.684 5 6 0.0 0.5 0.3260 position of multiplicity 2 8 0.0 0.5 0.5953 position of multiplicity 2 7 0.1459 0.6459 0.1766 position of multiplicity 4 1 0.2575 0.7575 0.2827 position of multiplicity 4 1 0.1441 0.6441 -0.0380 position of multiplicity 4 END 

3) Visualize the output (urea_geom.out). At the end of it, you should get

  LATTICE PARAMETERS (ANGSTROMS AND DEGREES) - PRIMITIVE CELL A B C ALPHA BETA GAMMA VOLUME 5.56500 5.56500 4.68400 90.00000 90.00000 90.00000 145.059850 COORDINATES OF THE EQUIVALENT ATOMS (FRACTIONARY UNITS) N. ATOM EQUIV AT. N. X Y Z 1 1 1 6 C 0.00000000000E+00 -5.00000000000E-01 3.26000000000E-01 2 1 2 6 C -5.00000000000E-01 0.00000000000E+00 -3.26000000000E-01 3 2 1 8 O 0.00000000000E+00 -5.00000000000E-01 -4.04700000000E-01 4 2 2 8 O -5.00000000000E-01 0.00000000000E+00 4.04700000000E-01 5 3 1 7 N 1.45900000000E-01 -3.54100000000E-01 1.76600000000E-01 6 3 2 7 N -1.45900000000E-01 3.54100000000E-01 1.76600000000E-01 7 3 3 7 N -3.54100000000E-01 -1.45900000000E-01 -1.76600000000E-01 8 3 4 7 N 3.54100000000E-01 1.45900000000E-01 -1.76600000000E-01 9 4 1 1 H 2.57500000000E-01 -2.42500000000E-01 2.82700000000E-01 10 4 2 1 H -2.57500000000E-01 2.42500000000E-01 2.82700000000E-01 11 4 3 1 H -2.42500000000E-01 -2.57500000000E-01 -2.82700000000E-01 12 4 4 1 H 2.42500000000E-01 2.57500000000E-01 -2.82700000000E-01 13 5 1 1 H 1.44100000000E-01 -3.55900000000E-01 -3.80000000000E-02 14 5 2 1 H -1.44100000000E-01 3.55900000000E-01 -3.80000000000E-02 15 5 3 1 H -3.55900000000E-01 -1.44100000000E-01 3.80000000000E-02 16 5 4 1 H 3.55900000000E-01 1.44100000000E-01 3.80000000000E-02 

Note that from the 5 irreducible atoms (given in input) the others have been generated by symmetry.

4) Check the volume and the density of your system with those reported in the table above.

5) Visualize the crystal structure with moldraw.

The origin of the lattice cell coincides with the origin of the symmetry elements in this representation. This particular choice implies the presence of fragments of molecules inside the unit cell.

6) Copy your urea_geom.d12 into urea_mol.d12, insert the keyword ATOMORDE, and run the new input: the atoms are reordered following chemical connectivity, so allowing identification of molecular fragments in the structure.

Can you identify molecular fragments?
How many of them do you get?
How are the atoms numbered now?

  ******************************************************************************* * SEARCHING THE MOLECULES OF THE CRYSTAL ******************************************************************************* MOLECULAR FRAGMENT N. 1 REFERENCE ATOM: 1 IN CELL 0 0 0 COORDINATES 0.000 -2.783 1.527 FRAGMENT COMPOSITION: SEQ. NO. NEW OLD CELL AT. NO. COORD.(A) 1 1 0 0 0 6 0.0000 -2.7825 1.5270 2 3 0 0 1 8 0.0000 -2.7825 2.7884 3 5 0 0 0 7 0.8119 -1.9706 0.8272 4 6 0 -1 0 7 -0.8119 -3.5944 0.8272 5 13 0 0 0 1 0.8019 -1.9806 -0.1780 6 9 0 0 0 1 1.4330 -1.3495 1.3242 7 14 0 -1 0 1 -0.8019 -3.5844 -0.1780 8 10 0 -1 0 1 -1.4330 -4.2155 1.3242 MOLECULAR FRAGMENT N. 2 REFERENCE ATOM: 2 IN CELL 0 0 0 COORDINATES -2.783 0.000 -1.527 FRAGMENT COMPOSITION: SEQ. NO. NEW OLD CELL AT. NO. COORD.(A) 9 2 0 0 0 6 -2.7825 0.0000 -1.5270 10 4 0 0 -1 8 -2.7825 0.0000 -2.7884 11 7 0 0 0 7 -1.9706 -0.8119 -0.8272 12 8 -1 0 0 7 -3.5944 0.8119 -0.8272 13 15 0 0 0 1 -1.9806 -0.8019 0.1780 14 11 0 0 0 1 -1.3495 -1.4330 -1.3242 15 16 -1 0 0 1 -3.5844 0.8019 0.1780 16 12 -1 0 0 1 -4.2155 1.4330 -1.3242 ******************************************************************************* 

Each molecule has 4 atoms in the reference cell and 4 (grey background) in neighbouring cells.

For molecular crystals an important CRYSTAL tool is the extraction of one (or more) molecules from the bulk structure on the basis of chemical connectivity, defined by the sum of covalent radii.
In CRYSTAL this task is accomplished by the MOLECULE keyword. See CRYSTAL User's Manual.

### Exercise: Silicon

1) Consider the following geometry input for silicon where the standard CRYSTAL origin setting (second choice according to the "International tables for Crystallography") has been chosen (see CRYSTAL User's Manual index, "origin setting"):

 SILICON BULK CRYSTAL 0 0 0 227 5.42 1 14 .125 .125 .125 TEST END 

2) As done in the previous cases, run a TEST and visualize the structure (2 atoms per cell).

### Exercise: Berillium

1) Consider the following geometry input for beryllium:

 BERILLIUM BULK CRYSTAL 0 0 0 194 2.29 3.59 1 4 0.333333333333 0.666666666667 .25 END 

2) As done in the previous cases, run a TEST and visualize the structure.

## 2. Block 2: Basis set

CRYSTAL performs ab initio calculations on periodic systems within the linear combination of Bloch functions, defined in terms of local functions, hereafter indicated as atomic orbitals (AO). These local functions are expressed as linear combination of a certain number of Gaussian type functions (GTF). They are characterized by the same centre with fixed coefficients and exponents defined in the input. The AOs belonging to a given atom are grouped into  shells : s, sp, p, d, f shells are allowed..

 shell type code shell type AO AO order max shell charge 0 s 1 s 2 1 sp 4 s, x, y, z 8 (s2 p6) 2 p 3 x, y, z 6 3 d 5 2z2-x2-y2,  xz,  yz,  x2-y2, xy 10 4 f 7 (2z2-3x2-3y2)z, (4z2-x2-y2)x, (4z2-x2-y2)y, (x2-y2)z, xyz, (x2-3y2)x, (3 x2-y2)y 0 - polarization only

For each atom (as many blocks as different types of atoms in the crystal structure) it must be specified:
the conventional atomic number  and the number of shells ns of the atomic basis set

for each shell (ns blocks of records): type of basis set (0-1-2), type of shell (0-1-2-3-4),
number of primitives GTF ng, shell electronic charge, scale factor

for each primitive (ng records - optional - to define for  basis set type 0 only):
exponent, contraction coefficient, [contraction coefficient]

The number of electrons assigned to each atom is the sum of the electrons in each shell, as given in input. This has two implications:

1. a maximum number of electrons can be attributed to each shell (see table above);
2. the atomic configuration (neutral atoms or ions, as in MgO input example) in the initial SCF guess (default choice: density matrix as superposition of atomic densities) is defined by the atomic shells population.

CRYSTAL can use either general basis sets, including s, p, d, f (polarization only) functions or standard Pople basis sets (internally stored).

 BS input code description 0 general basis set; gaussians exponent and contraction coefficients defined in input 1 Pople STO-nG type basis set - internally stored 2 Pople 3(6)-21G type basis set- internally stored

To better understand the input structure of the basis set in CRYSTAL input, we report different Basis Sets (BS) for silicon in the following table:

 STO-3G 6-21G modified 3-21G modified+polarization free basis set 14 3 1 0 3 2. 0. 1 1 3 8. 0. 1 1 3 4. 0.  14 4 2 0 6 2. 1. 2 1 6 8. 1. 2 1 2 4. 1. 0 1 1 0. 1. 0.16 1. 1.  14 5 2 0 3 2. 1. 2 1 3 8. 1. 2 1 2 4. 1. 0 1 1 0. 1. 0.16 1. 1. 0 3 1 0. 1. 0.5 1.  14 4 0 0 6 2. 1. 16115.9 0.00195948 2425.58 0.0149288 553.867 0.0728478 156.340 0.24613 50.0683 0.485914 17.0178 0.325002 0 1 6 8. 1. 292.718 -0.00278094 0.00443826 69.8731 -0.0357146 0.0326679 22.3363 -0.114985 0.134721 8.15039 0.0935634 0.328678 3.13458 0.603017 0.449640 1.22543 0.418959 0.261372 0 1 2 4. 1. 1.07913 -0.376108 0.0671030 0.302422 1.25165 0.956883 0 1 1 0. 1. 0.123 1. 1.  Number of shells: 3 Number of shells: 4 Number of shells: 5 Number of shells: 4 Number of AO: 9 Number of AO: 13 Number of AO: 18 Number of AO: 13 # electrons:14 # electrons:14 # electrons:14 # electrons:14
• STO-3G - standard Pople BS - . Minimal BS - the value of the scale factor "0." implies standard  STO-nG scale factor.
• 6-21G modified - Pople split valence 6-21G BS - the outer exponent of sp shell is modified ( 9.334E-02 for standard molecular calculations, 0.16 for crystal calculation)
• 3-21G modified+polarization - Pople split valence 3-21G BS - outer sp shell exponent modified as above - 1 d polarization function added
• free basis set - free split valence BS

The definition of the basis set for all the atoms ends with the records (mandatory):

 99 0 END 

### Exercise: Silicon

1) Write the input for silicon (geometry + basis set) with a basis set taken from the CRYSTAL website (Si_88-31G*_nada_1996).

2) Run a TEST based on this input.

3) How many functions (AO) do you get with the proposed basis set?

### Exercise: MgO

1) Write the input for MgO (geometry + basis set) with the following basis sets:

 Magnesium - ionic configuration, Mg 2+: 1s(2) 2sp(8) 12 3 0 0 8 2. 1. 68371.875 0.0002226 9699.34009 0.0018982 2041.176786 0.0110451 529.862906 0.0500627 159.186000 0.169123 54.6848 0.367031 21.2357 0.400410 8.74604 0.14987 0 1 6 8. 1. 156.795 -0.00624 0.00772 31.0339 -0.07882 0.06427 9.6453 -0.07992 0.2104 3.7109 0.29063 0.34314 1.61164 0.57164 0.3735 0.64294 0.30664 0.23286 0 1 1 0. 1. 3sp shell charge 0. 0.4 1. 1. Oxygen - ionic configuration, O2-: 1s(2) 2sp(8) 8 3 0 0 8 2. 1. 4000. 0.00144 1355.58 0.00764 248.545 0.05370 69.5339 0.16818 23.8868 0.36039 9.27593 0.38612 3.82034 0.14712 1.23514 0.07105 0 1 5 8. 1. 2sp shell charge 8. 52.1878 -0.00873 0.00922 10.3293 -0.08979 0.07068 3.21034 -0.04079 0.2043 1.23514 0.37666 0.34958 0.536420 0.42248 0.27774 0 1 1 0. 1. 0.210000 1. 1.   Magnesium - neutral configuration: 1s(2) 2sp(8) 3sp(2) 12 3 0 0 8 2. 1. 68371.875 0.0002226 9699.34009 0.0018982 2041.176786 0.0110451 529.862906 0.0500627 159.186000 0.169123 54.6848 0.367031 21.2357 0.400410 8.74604 0.14987 0 1 6 8. 1. 156.795 -0.00624 0.00772 31.0339 -0.07882 0.06427 9.6453 -0.07992 0.2104 3.7109 0.29063 0.34314 1.61164 0.57164 0.3735 0.64294 0.30664 0.23286 0 1 1 2. 1. 3sp shell charge 2. 0.4 1. 1. Oxygen neutral configuration: 1s(2) 2sp(6) 8 3 0 0 8 2. 1. 4000. 0.00144 1355.58 0.00764 248.545 0.05370 69.5339 0.16818 23.8868 0.36039 9.27593 0.38612 3.82034 0.14712 1.23514 0.07105 0 1 5 6. 1. 2sp shell charge 6. 52.1878 -0.00873 0.00922 10.3293 -0.08979 0.07068 3.21034 -0.04079 0.20433 1.23514 0.37666 0.34958 0.536420 0.42248 0.27774 0 1 1 0. 1. 0.210000 1. 1. 

10 electrons are assigned to O (1s2, (2sp)8), 10 electrons to Mg (1s2, (2sp)8); the SCF guess is a superposition of ionic densities, Mg2+ O2-
To obtain a superposition of atomic densities, the shell charges should be: O (1s2, (2sp)6) and Mg (1s2, (2sp)8(3sp)2)

2) Run a TEST based on this input.

3) How many functions (AO) do you get with the proposed basis set?

In the next section we will use these inputs to run single point calculations on MgO.

### Exercise: Urea

1) Write the input for urea (geometry + basis set) with a basis set taken from the CRYSTAL website (use basis sets by Gatti - C, N, O, H).

2) Run a TEST on this input.

3) How many functions (AO) do you get running the program with the proposed basis set?

## Block 3:  Hamiltonian and SCF computational parameters

The Hamiltonian

The hamiltonian specification is in the third input block. The hamiltonians implemented in CRYSTAL are HF (Hartree-Fock) and DFT (Density Functional Theory). The default choice is RHF (Hartree Fock - Closed shell). All keywords in the third input block are optional , except SHRINK (not for molecules) and END keyword, that must be inserted to close the section.

The basic choices for the hamiltonian are reported in the following table:

 Wavefunction HF DFT Characteristics Closed default DFT one determinant; Open UHF SPIN two determinants; one for alpha and one for beta electrons

DFT (Density Functional Theory) Hamiltonian

The keyword DFT selects a DFT Hamiltonian. Exchange-correlation functionals are separated in an exchange component (keyword EXCHANGE) and a correlation component (keyword CORRELAT). If the correlation functional is not set, an exchange-only functional is used in the hamiltonian. If the exchange functional is not set, the Hartree-Fock exchange potential is used.

Numerical integration is used (default integration parameters can be modified by the keywords LGRID, XLGRID, XXLGRID) ) .

The DFT input block must terminate with the keyword END[DFT].

Among DFT methods, various functionals are available in CRYSTAL:

• Local (-density approximation, LDA): functionals of the electron density only;
• meta-GGA: as for GGA, but also including either the kinetic energy density or the Laplacian of the electron density.
• Global-Hybrid (GH): the exchange functional is a linear combination of Hartree-Fock and DFT (LDA or GGA or mGGA) exchange term.
• Range-Separated Hybrid (RSH): hybrid functionals in which the amount of HF exchange included depends on the distance between electrons. Short-, long- and middle-range hybdrid functional are available
• Double-Hybrid (DH): Along with the hybridization of the exchange term, it combines also a DFT correlation functional with an a-posteriori PT2 correlation contribution (see the CRYSCOR tutorial for further details).

See CRYSTAL User's Manual for a complete list of available exchange-correlation functionals.

Hybrid functionals such as B3LYP and B3PW (A.D.Becke, J. Chem. Phys. 98, 5648, (1993)) are global keywords that define both exchange and correlation functionals.

Execution techniques

CRYSTAL can run in different ways:

1) traditional SCF: default, mono- and two-electron integrals are stored on disk and read-in at each SCF iteration;
2) direct SCF: keyword SCFDIR, the calculation is repeated at each iteration of SCF cycle.

Before running the complete calculation, evaluation of the disk space to be allocated for your run is recommended (keyword TESTRUN in 3rd input block), to choose the type of SCF, traditional (default) or direct (keyword SCFDIR).  The estimated number of mono- and bi-electron integrals to be stored is printed at the bottom of a testrun output (1 word = 8 byte).

SCF (Self Consistent Field) computational parameters

For periodic systems (1D, 2D, 3D) the mandatory input information is SHRINK, followed by the shrinking factors IS and ISP.
IS is used to generate a commensurate grid of k points in the reciprocal space, according to Pack-Monkhorst method. The hamiltonian matrix computed in the direct space is Fourier transformed for each k value and diagonalized to obtain eigenvectors and eigenvalues.

A second shrinking factor, ISP, defines the sampling of k points used for the determination of Fermi energy. For conducting systems, ISP should be at least 2*IS

Important keywords are:

 TOLDEE ite keyword SCF stops when absolute value of energy difference is less than 10-ite [default value is 6]
 MAXCYCLE imax keyword SCF stops when the number of cycles exceeds imax [default value is 50].

A number of convergence tools are available. See http://www.crystal.unito.it/tutorials.

 Single point calculations

The initial SCF guess is the density matrix of the isolated atoms in the cell. The atomic configuration can be ionic or neutral.

Three checks have to be done before proceeding with the calculation:

1. Check the geometry of the system (atom-atom distance, density, volume of the cell)
2. Check the basis set (avoid linear dependence - see EIGS);
3. Estimate the disk space to store integrals (TESTRUN in the 3rd input block) and choose traditional SCF (default) or direct SCF (keyword SCFDIR)
Run crystal (remove any TEST from your input deck), and check that SCF ends with the string:
 == SCF ENDED - CONVERGENCE ON ENERGY      E(AU) -x.xxxxxxxxxxxxE+yy CYCLES  nn


Verify the convergence of the SCF, restarting the calculation (without any convergence tools) from the density matrix (keyword GUESSP) obtained previously. If the program stops after 2 SCF cycles, convergence has been reached properly.

### Exercise: MgO bulk: check of computational parameters: k-points net, energy tolerance, integration grid (DFT hamiltonian only)

1) Run a single point calculation on MgO with the basis sets for Mg and O used before (see Basis set exercise: MgO). Use the following specifications: HF hamiltonian (default choice),

SHRINK
4 4
FMIXING
30


2) Analyze your output: how many SCF cycles did your system need to reach convergence? How many k points do you have in the reciprocal space? How ionic is the system (look at the electronic charge of all atoms printed at the last SCF cycle)? Which is the energy of the system?

3) Repeat the calculation increasing the shrinking factor as follows: from 4 to 6; from 6 to 8; from 8 to 12. Fill in the following table with the collected data and try to comment it.

 Shrinking factor Number of k points SCF energy (hartree)

4) Repeat the calculation using the following specifications: HF hamiltonian, SHRINK 8 8, FMIXING 30, TOLDEE 7 (default value 6). How many SCF cycles  to reach convergence? How has the energy of the system been affected by increasing the tolerance on the energy?

5) Repeat the previous calculation with TOLDEE 8. How many SCF cycles ar now needed to reach convergence? How is the energy of the system affected by increasing the tolerance on the energy?

6) Collect the results from the last 3 exercises (HF calculations with shrinking factor 8 8) and fill in the following table:

 TOLDEE SCF energy (hartree) 6 (default)

7) Run a new single point calculation using the following specifications: B3LYP hamiltonian, SHRINK 8 8, FMIXING 30.

8) Analize your output: how many SCF cycles were necessary to reach convergence? Which is the energy of the system? How much does it differ from point no.3? How ionic is the system now with respect to point no.3 (look at the electronic charges of the atoms in the cell printed at the last SCF cycle)?

9) With the same specifications used at point 7), run a single point calculation with slightly larger DFT integration grid by means of the keyword LGRID.

10) Repeat the previous exercise using a larger grid (XLGRID keyword).

11) Collect the results from the last 3 steps (DFT calculations with different integration grids) and fill in the following table:

 DFT integration grid Size of grid SCF energy (hartree) default

A few hints for the solutions are reported below

 SCF block ex. 1) SHRINK 4 4 FMIXING 30 END 

 SCF block ex. 3) DFT B3LYP END SHRINK 8 8 FMIXING 30 END 

 HF case (SHRINK 8 8): TOTAL ATOMIC CHARGES: 10.0220272 9.9779728 B3LYP case (SHRINK 8 8): TOTAL ATOMIC CHARGES: 10.0919882 9.9080118 

 Geometry optimization

The determination of equilibrium structure is of primary importance in the modelling of chemical systems.

CRYSTAL can compute analytical gradients of the energy with respect to both cell parameters and atom coordinates for 0, 1, 2, and 3-dimensional systems (insulators).

Geometry optimization can be performed in either symmetrized fractional coordinates or redundant internal coordinates.

The keyword OPTGEOM opens the structure optimization input (closed by END[OPT]). The keyword must be inserted in the geometry input section (block 1) as the last keyword of the section, before END[GEOM]

By default, the full relaxation of the atomic coordinates and lattice parameters is performed. When the position of all atoms is defined by the symmetry (like in MgO bulk), no optimization of atomic coordinates can be performed and only the lattice parameters are optimized.

Relaxation of the atomic positions at fixed lattice parameters is also possible with the keyword ATOMONLY. Symmetrized directions are generated according to the space group symmetry.

It is also possible to optimize the cell parameters at fixed atomic positions with the keyword CELLONLY. Symmetrized elastic distortions are generated according to the space group symmetry. They do not coincide with deformation of the crystallographic axes.

Optional choice (INTREDUN) is optimization in redundant internal coordinates of atomic positions and cell parameters.

Restart is possible for a job which is terminated abruptly (RESTART).

### Exercise: MgO bulk - optimization of cell parameter

1) Optimize the cell parameter of MgO bulk. Mg and O are in special symmetry points, their position can not be modified. For Mg and O use the basis sets taken from Basis set exercise: MgO. Use the following specifications: HF hamiltonian, SHRINK 8 8, FMIXING 30. The input deck is the following:

 MGO BULK - CELL OPT CRYSTAL 0 0 0 225 4.217 2 12 0. 0. 0. 8 0.5 0.5 0.5 OPTGEOM CELLONLY END END 12 3 0 0 8 2. 1. 68371.875 0.0002226 9699.34009 0.0018982 2041.176786 0.0110451 529.862906 0.0500627 159.186000 0.169123 54.6848 0.367031 21.2357 0.400410 8.74604 0.14987 0 1 6 8. 1. 156.795 -0.00624 0.00772 31.0339 -0.07882 0.06427 9.6453 -0.07992 0.2104 3.7109 0.29063 0.34314 1.61164 0.57164 0.3735 0.64294 0.30664 0.23286 0 1 1 0. 1. 0.4 1. 1. 8 3 0 0 8 2. 1. 4000. 0.00144 1355.58 0.00764 248.545 0.05370 69.5339 0.16818 23.8868 0.36039 9.27593 0.38612 3.82034 0.14712 1.23514 0.07105 0 1 5 8. 1. 52.1878 -0.00873 0.00922 10.3293 -0.08979 0.07068 3.21034 -0.04079 0.20433 1.23514 0.37666 0.34958 0.536420 0.42248 0.27774 0 1 1 0. 1. 0.210000 1. 1. 99 0 END SHRINK 8 8 FMIXING 30 END

2) How many optimization cycles did CRYSTAL perform?
Which is the optimized value for the cell parameter?
Which is the corresponding energy?

Here follows the solution to this exercise:

  ******************************************************************************* CRYSTALLOGRAPHIC CELL (VOLUME= 73.60933348) A B C ALPHA BETA GAMMA 4.19093535 4.19093535 4.19093535 90.000000 90.000000 90.000000 COORDINATES IN THE CRYSTALLOGRAPHIC CELL ATOM X/A Y/B Z/C ******************************************************************************* 1 T 12 MG 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 2 T 8 O -5.000000000000E-01 5.000000000000E-01 5.000000000000E-01 

  ****************************************************************** * OPT END - CONVERGED * E(AU): -2.746642641310E+02 POINTS 5 * ****************************************************************** 

### Exercise: Urea bulk and cohesive energy

Before starting this exercise, let us recall how to compute the cohesive energy in a molecular crystal.

The cohesive energy of a molecular crystal is computed as :

E = E(bulk)/Z - E(mol)

where E(bulk) is the total energy of the unit cell and must be referred to the number, Z, of molecules in the unit cell, and E(mol) is the total energy of the isolated molecule in the gas phase. It corresponds to the packing energy due to the interaction among the molecules in the crystal. Note that both molecule and crystal must be in their equilibrium structure. The computed E value can be compared with experimental sublimation energies.
E can be written as the sum of two contributions:

E = E(cond) + E(conf)

where the first term refers to the condensation of molecules that keep the same conformation they have in the crystal and this term corresponds to the lattice energy; the second term is the difference between the energies of the isolated molecules in the bulk and gas phase conformation. E(conf) is, of course, negligible for rigid molecules, but can be significant for floppy molecules. As opposed to the C2v geometry of molecular urea in the bulk, the most stable structure has a C2 symmetry with an anti configuration (see figure below). Therefore, when computing the interaction energy, E(conf) between the C2 and the C2v conformers must be taken into account. Computed cohesive energies have also to be corrected for the BSSE through the counterpoise method as accomplished by using the MOLEBSSE keyword. This error affects E(cond) and must be computed on the optimized bulk structure.

By including the BSSE correction, the interaction energy then becomes:

E = E(cond) + E(conf) + BSSE

where, for crystalline urea, we have

E(cond) = E(bulk)/2 - E(mol - C2v)

E(conf) = E(mol - C2v) - E(mol - C2)

BSSE = E(mol - bulk) - E(mol - ghosts)

because there are two molecules in the cell.

In summary, one needs:

• E(bulk) = total energy of the unit cell in its equilibrium geometry
• E(mol - C2v) = total energy of the molecule in its optimized geometry with the same conformation as in the bulk. In this case with a planar C2v conformation
• E(mol - C2) = total energy of the molecule in its optimized geometry with the same conformation as in the gas phase. C2 conformer in the present case
• E(mol - bulk) = total energy of the molecule with the same geometry as in the bulk.
• E(mol - ghosts) = total energy of the molecule as obtained by augmenting the basis set with the ghost functions of surrounding atoms. This is obtained by using the MOLEBSSE option.

All the requested calculations will be performed with the HF Hamiltonian and with 3-21G basis sets for all the atoms involved.
1. Run a full optimization of urea bulk structure. Analyze the output to locate the equilibrium geometry. Which is the minimum energy? How many optimization steps did the program perform?
2. Run a single-point total energy calculation on the urea molecule as extracted from the optimized structure of the bulk.
3. Compute the BSSE correction for the urea molecule in the bulk. Which energy do you get?
4. Optimize the structure of the urea molecule in the C2v conformation. This is the planar conformation of molecules in bulk.
5. Optimize the structure of the urea molecule in the C2 conformation. This is the non-planar conformation of molecules in gas phase.
6. Collect all the results and fill in the following table. Then, calculate the cohesive energy of crystalline urea.
 System Label SCF energy (hartree) E(bulk) E1 E(mol - C2v) E2 E(mol - C2) E3 E(mol - bulk) E4 E(mol - ghosts) E5

The geometry input blocks of this exercise are the following:

 Geometry block ex. 1) UREA BULK - OPTIMIZATION CRYSTAL 0 0 0 113 5.565 4.684 5 6 0.0 0.5 0.326 8 0.0 0.5 0.5953 7 0.1459 0.6459 0.1766 1 0.2575 0.7575 0.2827 1 0.1441 0.6441 -0.0380 OPTGEOM FULLOPTG END END 

 Geometry block ex. 2) UREA BULK - SINGLE-POINT MOLECULE IN THE BULK CRYSTAL 0 0 0 113 5.59037346 4.67007538 5 6 0.000000000000E+00 -5.000000000000E-01 3.176227372521E-01 8 0.000000000000E+00 -5.000000000000E-01 -4.093385097982E-01 7 1.437784027012E-01 -3.562215972988E-01 1.678775482565E-01 1 2.564186903752E-01 -2.435813096248E-01 2.695980532036E-01 1 1.406087223437E-01 -3.593912776563E-01 -4.827103388750E-02 MOLECULE 1 1 0 0 0 END 

 Geometry block ex. 3) UREA BULK - BSSE CRYSTAL 0 0 0 113 5.59037346 4.67007538 5 6 0.000000000000E+00 -5.000000000000E-01 3.176227372521E-01 8 0.000000000000E+00 -5.000000000000E-01 -4.093385097982E-01 7 1.437784027012E-01 -3.562215972988E-01 1.678775482565E-01 1 2.564186903752E-01 -2.435813096248E-01 2.695980532036E-01 1 1.406087223437E-01 -3.593912776563E-01 -4.827103388750E-02 MOLEBSSE 1 1 0 0 0 1 4.00 END 

 Geometry block ex. 4) UREA MOLECULE C2v MOLECULE 15 5 6 4.010040108509E-18 4.010040108509E-18 1.511178418937E-01 8 3.627793330734E-17 3.627793330734E-17 1.367129216024E+00 7 0.0000000000000000 1.152517719977E+00 -6.029124102443E-01 1 0.0000000000000000 1.186574472949E+00 -1.524738700435E+00 1 0.0000000000000000 2.002684092087E+00 -6.757791827922E-02 OPTGEOM END END 

 Geometry block ex. 5) UREA MOLECULE C2 MOLECULE 5 5 6 4.010040108509E-18 4.010040108509E-18 1.511178418937E-01 8 3.627793330734E-17 3.627793330734E-17 1.367129216024E+00 7 1.475136395636E-01 1.152517719977E+00 -6.029124102443E-01 1 -2.591127703870E-01 1.186574472949E+00 -1.524738700435E+00 1 6.956803083756E-02 2.002684092087E+00 -6.757791827922E-02 OPTGEOM END END 

Exercise: Recalculate the cohesive energy of crystalline urea using the following data obtained with the 6-31G(d,p) basis set:

 System Label SCF energy (hartree) E(bulk) E1 -4.4804947630236E+02 E(mol - C2v) E2 -2.2398536052639E+02 E(mol - C2) E3 -2.2399715236476E+02 E(mol - bulk) E4 -2.2399397220548E+02 E(mol - ghosts) E5 -2.2400050884526E+02

Compare the result with the previous calculation with 3-21G basis.

### Exercise: Optimization of urea bulk in redundant coordinates

1) Run a full geometry optimization for urea bulk using redundant internal coordinates (INTREDUN).

2) Analyze your output: which is the equilibrium geometry xxxxxxxx? How many optimization steps are necessary to reach the minimum? Compare it with the optimization with symmetrized fractional coordinates.

3) Check that the minimum energy you get from this optimization is the same as the one you got for symmetrized fractional coordinates.

 Vibrational frequencies

The crystal lattice is never rigid. Atoms actually move around their equilibrium positions inside the crystalline structure. Theoretical calculation of atom vibrations gives access to a number of properties such as the calculation of vibrational spectra, the prediction of structure stability or instability of a crystalline phase and the computation of thermodynamical properties (e.g. heat capacity, entropy, ...).

Vibration frequency calculations in CRYSTAL are performed at the $$\Gamma$$ point. The dynamical matrix is computed by numerical evaluation of the first-derivative of the analytical atomic gradients. The point group symmetry of the system is fully exploited to reduce the number of points (SCF+gradient) to be considered. On each numerical step, the residual symmetry is preserved during the SCF process and the gradients calculation.

The keyword FREQCALC opens the vibrational frequencies calculation input (closed by END). The keyword must be inserted in the geometry input section (block 1) as the last keyword of the section. It is worth noting that meaningful vibrational frequencies are only computed when the crystalline structure is a stationary point on the potential energy surface, generally a minimum.

Geometry optimization and frequencies calculation must be computed in two different runs.

Note that the transverse optical (TO) modes are computed by default.

Restart is possible for a job which is terminated abruptly (RESTART).

### Exercise: MgO

1) Run a vibrational frequencies calculation on MgO. (add FREQCALC at the end of geometry input block)

2) Analyze your output and locate where the relevant information on the vibrational frequencies calculation are reported:

a) how many SCF+gradient calculations were performed?;

b) to which irreducible representations belong the vibration modes?;

c) which modes are IR active and which modes are Raman active?

d) look for the vibrational temperatures;

e) analyze the thermodynamic properties of your system.

The solutions to this exercise are reported below:

 Geometry block ex. 1) MGO BULK - FREQUENCIES CRYSTAL 0 0 0 225 4.191 2 12 0. 0. 0. 8 0.5 0.5 0.5 FREQCALC frequency calculation with default computational parameters END END 

From output:

Numerical calculation of second derivatives from analytical gradients:

  NUMBER OF SCF+GRADIENT CALCULATIONS 2 MAX ABS(DGRAD): MAXIMUM ABSOLUTE GRADIENT DIFFERENCE WITH RESPECT TO THE CENTRAL POINT DE: ENERGY DIFFERENCE WITH RESPECT TO THE CENTRAL POINT (DE IS EXPECTED TO BE POSITIVE FOR ALL DISPLACEMENTS) ATOM MAX ABS(DGRAD) TOTAL ENERGY (AU) N.CYC DE SYM CENTRAL POINT -2.746642754933E+02 15 0.0000E+00 48 2 O DX 4.4659E-04 -2.746642742333E+02 6 1.2600E-06 8 2 O DY GENERATED FROM A PREVIOUS LINE 2 O DZ GENERATED FROM A PREVIOUS LINE 

Symmetry analysis:

  CLASS | GROUP OPERATORS (SEE SYMMOPS KEYWORD) -------------------------------------------------------------------- C2 | 2; 3; 4; C3 | 5; 11; 7; 9; 6; 12; 10; 8; C2' | 13; 14; 17; 18; 23; 24; C4 | 15; 16; 19; 20; 21; 22; I | 25; SGH | 26; 27; 28; S6 | 29; 35; 31; 33; 30; 36; 34; 32; SGD | 37; 38; 41; 42; 47; 48; S4 | 39; 40; 43; 44; 45; 46; IRREP/CLA E C2 C3 C2' C4 I SGH S6 SGD S4 ------------------------------------------------------------------------------- MULTIP | 1 3 8 6 6 1 3 8 6 6 ------------------------------------------------------------------------------- Fu | 3.00 -1.00 0.00 -1.00 1.00 -3.00 1.00 0.00 1.00 -1.00 Fu -(3, 2); 

Vibrational modes:

  MODES EIGV FREQUENCIES IRREP IR INTENS RAMAN (HARTREE**2) (CM**-1) (THZ) (KM/MOL) 1- 3 0.2379E-21 0.0000 0.0000 (Fu ) A ( 0.00) I 4- 6 0.4503E-05 465.7569 13.9630 (Fu ) A ( 0.00) I 

  VIBRATIONAL TEMPERATURES (K) [MODE NUMBER;IRREP] TO MODES 670.1 [ 4;Fu ] 670.1 [ 5;Fu ] 670.1 [ 6;Fu ] 

  HARMONIC VIBRATIONAL CONTRIBUTIONS TO THERMODYNAMIC FUNCTIONS AT GIVEN TEMPERATURE AND PRESSURE: (EL = ELECTRONIC ENERGY E0 = ZERO-POINT ENERGY ET = THERMAL CONTRIBUTION TO THE VIBRATIONAL ENERGY PV = PRESSURE * VOLUME TS = TEMPERATURE * ENTROPY) AU/CELL EV/CELL KJ/MOL EL : -274.664275493308 -7473.994906731630 -721130.95368190 E0 : 0.003183216619 0.086619727857 8.35753405 ******************************************************************************* THERMODYNAMIC FUNCTIONS WITH VIBRATIONAL CONTRIBUTIONS AT (T = 298.15 K, P = 0.10132500E+00 MPA): AU/CELL EV/CELL KJ/MOL ET : 0.000752096746 0.020465592918 1.97462973 PV : 0.000000427709 0.000011638560 0.00112295 TS : 0.001068385430 0.029072245567 2.80504555 ET+PV-TS : -0.000315860975 -0.008595014089 -0.82929287 EL+E0+ET+PV-TS: -274.661408137664 -7473.916882017863 -721123.42544072 OTHER THERMODYNAMIC FUNCTIONS: mHARTREE/(CELL*K) mEV/(CELL*K) J/(MOL*K) ENTROPY : 0.003583382292 0.097508789426 9.40816888 HEAT CAPACITY : 0.006339433012 0.172504742223 16.64417903 

### Exercise: calcite

1) Analyze the frequency output of calcite (CaCO3):

a) how many SCF+gradient calculations were performed?;

b) to which irreducible representations belong the vibration modes?;

c) which modes are IR active and which are Raman active?

d) look for the vibrational temperatures;

e) analyze the thermodynamic properties of your system.

The solutions to this exercise are reported below:

  MAX ABS(DGRAD): MAXIMUM ABSOLUTE GRADIENT DIFFERENCE WITH RESPECT TO THE CENTRAL POINT DE: ENERGY DIFFERENCE WITH RESPECT TO THE CENTRAL POINT (DE IS EXPECTED TO BE POSITIVE FOR ALL DISPLACEMENTS) ATOM MAX ABS(DGRAD) TOTAL ENERGY (AU) N.CYC DE SYM CENTRAL POINT -1.883021622778E+03 31 0.0000E+00 12 1 CA DX 3.6415E-04 -1.883021621785E+03 7 9.9308E-07 1 1 CA DY GENERATED FROM A PREVIOUS LINE 1 CA DZ 5.3227E-04 -1.883021621201E+03 7 1.5762E-06 3 3 C DX 4.0721E-03 -1.883021611246E+03 10 1.1531E-05 1 3 C DY GENERATED FROM A PREVIOUS LINE 3 C DZ 1.5595E-03 -1.883021618367E+03 9 4.4102E-06 3 5 O DX 2.4149E-03 -1.883021615750E+03 9 7.0277E-06 1 5 O DY GENERATED FROM A PREVIOUS LINE 5 O DZ 5.1208E-04 -1.883021621726E+03 9 1.0521E-06 1 

 LASS | GROUP OPERATORS (SEE SYMMOPS KEYWORD) -------------------------------------------------------------------- C3 | 2; 3; C2 | 4; 5; 6; I | 7; S6 | 8; 9; SGD | 10; 11; 12; IRREP/CLA E C3 C2 I S6 SGD --------------------------------------------------- MULTIP | 1 2 3 1 2 3 --------------------------------------------------- A1g | 1.00 1.00 1.00 1.00 1.00 1.00 A2g | 1.00 1.00 -1.00 1.00 1.00 -1.00 Eg | 2.00 -1.00 0.00 2.00 -1.00 0.00 A1u | 1.00 1.00 1.00 -1.00 -1.00 -1.00 A2u | 1.00 1.00 -1.00 -1.00 -1.00 1.00 Eu | 2.00 -1.00 0.00 -2.00 1.00 0.00 A1g-(1, 1); A2g-(1, 3); Eg -(2, 4); A1u-(1, 2); A2u-(1, 4); Eu -(2, 6); 

  MODES EIGV FREQUENCIES IRREP IR INTENS RAMAN (HARTREE**2) (CM**-1) (THZ) (KM/MOL) 1- 2 -0.2676E-21 0.0000 0.0000 (Eu ) A ( 0.00) I 3- 3 0.3271E-21 0.0000 0.0000 (A2u) A ( 0.00) I 4- 5 0.3172E-06 123.6060 3.7056 (Eu ) A ( 144.60) I 6- 6 0.3602E-06 131.7273 3.9491 (A2u) A ( 104.21) I 7- 8 0.4969E-06 154.7137 4.6382 (Eg ) I ( 0.00) A 9- 9 0.8120E-06 197.7660 5.9289 (A2g) I ( 0.00) I 10- 11 0.1004E-05 219.8920 6.5922 (Eu ) A ( 189.99) I 12- 13 0.1613E-05 278.7694 8.3573 (Eg ) I ( 0.00) A 14- 15 0.1676E-05 284.1535 8.5187 (Eu ) A ( 812.65) I 16- 16 0.1744E-05 289.8341 8.6890 (A1u) I ( 0.00) I 17- 17 0.1842E-05 297.8520 8.9294 (A2u) A ( 347.09) I 18- 18 0.2028E-05 312.5693 9.3706 (A2g) I ( 0.00) I 19- 20 0.1036E-04 706.3442 21.1757 (Eg ) I ( 0.00) A 21- 22 0.1037E-04 706.7574 21.1881 (Eu ) A ( 39.53) I 23- 23 0.1564E-04 867.9843 26.0215 (A2u) A ( 205.01) I 24- 24 0.1627E-04 885.2715 26.5398 (A2g) I ( 0.00) I 25- 25 0.2442E-04 1084.5219 32.5131 (A1g) I ( 0.00) A 26- 26 0.2444E-04 1084.9299 32.5254 (A1u) I ( 0.00) I 27- 28 0.4216E-04 1425.1279 42.7243 (Eu ) A ( 5579.85) I 29- 30 0.4408E-04 1457.1303 43.6837 (Eg ) I ( 0.00) A 

  VIBRATIONAL TEMPERATURES (K) [MODE NUMBER;IRREP] TO MODES 177.8 [ 4;Eu ] 177.8 [ 5;Eu ] 189.5 [ 6;A2u] 222.6 [ 7;Eg ] 222.6 [ 8;Eg ] 284.5 [ 9;A2g] 316.4 [ 10;Eu ] 316.4 [ 11;Eu ] 401.1 [ 12;Eg ] 401.1 [ 13;Eg ] 408.8 [ 14;Eu ] 408.8 [ 15;Eu ] 417.0 [ 16;A1u] 428.5 [ 17;A2u] 449.7 [ 18;A2g] 1016.3 [ 19;Eg ] 1016.3 [ 20;Eg ] 1016.9 [ 21;Eu ] 1016.9 [ 22;Eu ] 1248.8 [ 23;A2u] 1273.7 [ 24;A2g] 1560.4 [ 25;A1g] 1561.0 [ 26;A1u] 2050.4 [ 27;Eu ] 2050.4 [ 28;Eu ] 2096.5 [ 29;Eg ] 2096.5 [ 30;Eg ] 

  HARMONIC VIBRATIONAL CONTRIBUTIONS TO THERMODYNAMIC FUNCTIONS AT GIVEN TEMPERATURE AND PRESSURE: (EL = ELECTRONIC ENERGY E0 = ZERO-POINT ENERGY ET = THERMAL CONTRIBUTION TO THE VIBRATIONAL ENERGY PV = PRESSURE * VOLUME TS = TEMPERATURE * ENTROPY) AU/CELL EV/CELL KJ/MOL EL : -1883.021622777611 -51239.623327891757 -4943872.57388462 E0 : 0.036144143164 0.983532137296 94.89643450 ******************************************************************************* THERMODYNAMIC FUNCTIONS WITH VIBRATIONAL CONTRIBUTIONS AT (T = 298.15 K, P = 0.10132500E+00 MPA): AU/CELL EV/CELL KJ/MOL ET : 0.008623804538 0.234665651658 22.64179562 PV : 0.000002966432 0.000080720717 0.00778837 TS : 0.015372123347 0.418296742062 40.35950416 ET+PV-TS : -0.006745352377 -0.183550369687 -17.70992017 EL+E0+ET+PV-TS: -1882.992223986824 -51238.823346124140 -4943795.38737028 OTHER THERMODYNAMIC FUNCTIONS: mHARTREE/(CELL*K) mEV/(CELL*K) J/(MOL*K) ENTROPY : 0.051558354341 1.402974147449 135.36644025 HEAT CAPACITY : 0.051280212801 1.395405531372 134.63617974 

 One-electron properties

At the end of the SCF process, the program crystal writes information on the crystalline system and its wave function as unformatted sequential data on file fort.9, and as formatted data on file fort.98. Formatted data allow transfer from one platform to another: for instance,  optimization of a structure or frequency calculation on a powerful, parallel system (crystal or Pcrystal program), analysis of the wave function (properties program) and visualization on your PC.

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

The wave function data  are read (when fort.98 is used, the first keyword must be RDFMWF, to convert formatted data in binary data) when running properties, and are not modified.

The data include:

• Crystal structure, geometry and symmetry operators.
• Basis set.
• Reciprocal lattice k-points sampling information.
• Irreducible Fock/KS matrix in direct space.
• Irreducible density matrix in direct space.

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

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

• Atoms and bond populations (Mulliken Population Analysis) (keyword PPAN)
• Charge Electron Spin Density (keyword ECHG)
• Band Structure (keyword BAND)
• Density of States (keyword DOSS - preliminary computation of eigenvectors required)

A complete list of properties is presented in http://www/crystal.unito.it/features.html

To run a property calculation:

runprop inputfile wavefunction_filename

(where inputfile.d3 is the input to properties and wavefunction_filename.f9 contains the wave function of a preliminary and mandatory crystal calculation).

Compute the wave function for MgO, Si, Be, urea bulk:

runcry mgo
runcry si
runcry be
runcry urea

the script runcry saves wf data in file mgo.f9, si.f9, be.f9, urea.f9

### Exercise: Mulliken population analysis

A Mulliken population analysis can be performed by inserting the keyword PPAN. For each non-equivalent atom, information on the electronic charge of atoms, shells, AOs, and bond populations up to the first six neighbours are printed.

1) Perform Mulliken population analysis of the charge density of MgO bulk:

runprop mgo_mul_pop    mgo.

The properties input file, mgo_mul_pop.d3, looks like

 PPAN END 

  ALPHA+BETA ELECTRONS MULLIKEN POPULATION ANALYSIS - NO. OF ELECTRONS 20.000000 ATOM Z CHARGE A.O. POPULATION 1 MG 12 10.022 2.000 1.999 1.995 1.995 1.995 0.006 0.010 0.010 0.010 2 O 8 9.978 2.007 1.177 1.146 1.146 1.146 0.820 0.845 0.845 0.845 ATOM Z CHARGE SHELL POPULATION 1 MG 12 10.022 2.000 7.986 0.036 2 O 8 9.978 2.007 4.616 3.355 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.984 2.108 -0.010 1 MG ( -1 0 0) 5.635 2.982 0.000 2 O ( 0 0 0) 6.901 3.652 0.000 1 MG ( -1 -1 1) 7.969 4.217 0.000 2 O ( 1 0 -1) 8.910 4.715 0.000 1 MG ( -2 0 1) 9.760 5.165 0.000 ATOM A 2 O ATOM B CELL R(AB)/AU R(AB)/ANG OVPOP(AB) 1 MG ( 1 0 0) 3.984 2.108 -0.010 2 O ( -1 0 0) 5.635 2.982 -0.012 1 MG ( 0 0 0) 6.901 3.652 0.000 2 O ( -1 -1 1) 7.969 4.217 0.000 1 MG ( -1 0 1) 8.910 4.715 0.000 2 O ( -2 0 1) 9.760 5.165 0.000 

MgO is an ionic system, 10.022 and 9.978  electrons are attributed to Mg and O, respectively.

2) Perform Mulliken population analysis of the charge density of Si bulk.

3) Perform Mulliken population analysis of the charge density of urea bulk.

4) Localize the 4 hydrogen bonds linking each urea molecule in the urea crystal looking at the bond population and to the X-H.....Y distance.

## Electron charge density

### Exercise: Electron charge density of MgO bulk

The total electron charge density maps (keyword ECHG) provide a pictorial representation of the total electronic distribution.
The charge density is computed in a grid of points, defined according to instructions given at the entry MAPNET of CRYSTAL User's Manual

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 may be 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 compute a density matrix as superposition of atom densities).

Charge density at grid points are saved in the file input_filename.f25.

In general, you have to follow the scheme reported in the table that refers to MgO:

 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. -4. 0.0 cartesian coordinates of point A (Angstrom) 4. -4. 0.0 cartesian coordinates of point B (Angstrom) 4. 4. 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

1) Insert in  properties input the keyword PATO, to compute a total density matrix as superposition of atomic densities, followed by the same ECHG block defined above. Your input file is

 ECHG 0 65 COORDINA -4. -4. 0.0 4. -4. 0.0 4. 4. 0.0 MARGINS 1.5 1.5 1.5 1.5 END PATO 0 0 ECHG 0 65 COORDINA -4. -4. 0.0 4. -4. 0.0 4. 4. 0.0 MARGINS 1.5 1.5 1.5 1.5 END END   A | | | | B___________________C The position of the three points is given by their cartesian coordinates 

2) In order to visualize the total electron charge density map and the difference electron density map, type:

maps  input_filename

where  input_filename.f25  is the file fort.25 saved by the script runprop at the end of properties calculation.
Input parameters default values to draw the maps are supplied.
Two maps are drawn in one postscript file, showing the crystal density and the superposition of ionic densities (as in the basis set input the configuration  defined is Mg2+ O2-)

See Crgra2006 user manual  to draw difference maps.

The total density and difference density maps look like:

### Exercise:  Electron charge density of Si bulk

In the case of Silicon the plane is defined with reference to three atoms in the crystal (keyword ATOMS of MAPNET). This choice is independent from the geometry of the system, but the length of AB and BC is modified, according to the positions of the atoms (after a geometry optimization the atoms may be moved from the original position):

 ECHG keyword 0 charge density only 65 ATOMS coordinates of point A,B,C given by the position of 3 atoms 1    1 1 -1 Point A - first atom in the cell (1 1 -1) 1    0 0 0 Point B - first atom in the reference crystal cell 1    0 0 2 Point C - first atom in the cell  (0 0 2) MARGINS definition of a larger window 2.0 2.0 2.0 2.0 Width of margins, in the order AB,CD,AD,BC END

1) Using the wave function for Silicon , run properties to obtain a total electron density map.

2) Perform the difference maps by appending the part relative to the atomic superposition to the previous input

3) Visualize and compare the electron charge density maps.

The solutions to this exercise are the following:

1. Electron charge density maps input file:

 ECHG 0 65 ATOMS 1 1 1 -1 1 0 0 0 1 0 0 2 MARGINS 2. 2. 2. 2. END PATO 0 0 ECHG 0 65 ATOMS 1 1 1 -1 1 0 0 0 1 0 0 2 MARGINS 2. 2. 2. 2. PPAN END 

### Exercise: Charge electron density of Urea bulk

In the case of urea bulk the plane is defined with reference to three atoms in the crystal (keyword ATOMS).

The definition of the coordinates of points A, B, C is:

 ECHG keyword 0 charge density only 65 ATOMS coordinates of point A,B,C given by the position of 3 atoms 2    0 0 1 Point A - first atom in the cell ( 1 1 -1) 2    0 0 0 Point B - first atom in the reference crystal cell 2    1 1 0 Point C - first atom in the cell  ( 0 0 2) RECTANGU The window is modified (larger) to be rectangular MARGINS definition of a larger window 2.0 2.0 2.0 2.0 Width of margins, in the order AB,CD,AD,BC PPAN Mulliken population analysis END

Note - The three points A,B,C define a parallelogram. The keyword RECTANGU defines the window as a rectangular inscribing the parallelogram.

1) Using the wave function for urea, run properties to obtain a total electron charge density map  in the plane defined by the second atom in the cell (100), the second atom in the reference cell and the second atom in the cell (110).

3) Recognize where the hydrogen bonds are located, with the help of Mulliken population analysis (value of overlap population between X--H------Y.

The solutions to this exercise are the following:

Electron density maps input file:

 ECHG 0 65 ATOMS 2 0 0 1 2 0 0 0 2 1 1 0 RECTANGU MARGINS 2. 2. 2. 2. PPAN END 

## Band Structure calculations

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 the input file for the band structure (see also the CRYSTAL User's Manual, keyword BAND) some 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).

2. The range of bands to be written in fort.25 for plotting.
Suggestion: all (the number of bands is the number of AO). Computing bands can be very demanding for large system: all eigenvalues are computed, avoid to miss information.

In order to compute the bands, the following input data are required:
1. The number of segments defining the path in the reciprocal space;
2. The shrinking factor in term of which the coordinates of the extreme of the segments are expressed (suggestion: 12);
3. The total number of k points along the path;
4. The first band to be saved;
5. The last band to be saved;
6. The plotting option (1 writes data on fortran unit 25);
7. The printing options (0 for no printing).

### Example: MgO bulk

Let us consider 5 segments along the following path (see figure below):

 Gamma(0,0,0) --> X(1/2,0,1/2) --> W(1/2,1/4,3/4) --> L(1/2,1/2,1/2) --> Gamma(0,0,0) --> W(1/2,1/4,3/4)

Using 12 as shrinking factor, we obtain the coordinates of the above special lattice points reported in the commented input deck below.
For example, the three integers 6 3 9 in the table correspond to point W with coordinates (6/12 3/12 9/12), or equivalently (1/2 1/4 3/4) in reciprocal space.

We want to plot ALL the bands, i.e. we need to consider ALL the AO (see SCF output obtained previously).

The paths in the Brillouin zone can be found e.g. in M. Lax, Symmetry Principles in Solid State and Molecular Physics (Dover Publications, 2001.

Here is reported and commented the input deck to compute MgO band structure:

 BAND keyword MgO title 5 12 30 1 18 1 0 5: number of segments in the reciprocal space to explore; 12: shrinking factor in term of which the coordinates of the extreme of the segments are expressed; 30: 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    6 0 6 segment in the reciprocal space: X 6 0 6    6 3 9 X W 6 3 9    6 6 6 W L 6 6 6    0 0 0 L 0 0 0    6 3 9

1) Run a properties calculation with the following input (mgo_band.d3):

 BAND MGO 5 12 30 1 18 1 0 0 0 0 6 0 6 6 0 6 6 3 9 6 3 9 6 6 6 6 6 6 0 0 0 0 0 0 6 3 9 END 

runprop mgo_band mgo

2) Check if file mgo_band_mgo.f25 has been created.

3) Visualize bands: band mgo_band_mgo. Default control file is supplied.

### Exercise: Silicon bulk

1. Write the input file for the Si band structure, si_band.d3, considering that:
• It is a cubic face centered lattice so its reciprocal lattice is body-centered cubic.
• A convenient choice of path in the reciprocal space could be = (0,0,0), X = (1/2, 0, 1/2), W = (1/2, 1/4, 3/4) and L = (1/2, 1/2, 1/2)
• Core electrons are 20 (because there are two Si in the unit cell), valence electrons 8; band range 1-18 will save eigenvalues for all occupied bands + 4 virtual
2. Run the properties calculation: runprop si_band si.
3. Check if the output is correct and if file si_band_si.f25) has been created.
4. Visualize bands: band si_band_si (the file si_band_si.band06.ps is the postscript file generated)

The solutions to this exercise are the following:

• Si band structure: input deck si_band.d3

 BAND Si 4 12 120 11 18 1 0 0 0 0 6 0 6 6 0 6 6 3 9 6 3 9 6 6 6 6 6 6 0 0 0 END 

• ## Density of States

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 from a given hamiltonian matrix.

In the input file for the density of states (DOSS) calculation (see main tutorials) the following information must be specified:

• The range of energy to compute DOSS. It can be defined by two bands (the range is from the bottom of the first band to the top of the second band), or by two energies - In case of very flat bands, better energy range.;
• The degree of the polynomials to be used in the DOS expansion.

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

• The label of the atom(s) (in the unit cell order) on which the projection has to be performed
• For each atom(s) the sequential numbers of the atomic orbitals to include in the projection (look at the output file of wave function calculation).
• To project the DOS onto the whole set of AOs of an atoms, write a -1 followed by the sequential number of the selected atom in the unit cell).

### Exercise: MgO bulk

The density of states is projected onto the whole set of AOs of the atom 1 (Mg) and of the atom 2 (O). File mgo_doss.d3:

 NEWK keyword 0 0 0: shrinking factor for reciprocal space Pack-Monkhorst net is taken from the SCF input 0: shrinking factor for reciprocal space Gilat net is taken from the SCF input 1 1 1: calculation of the density matrix from the new eigenvectors (it corresponds to one more SCF cycle) 1: print options 66 100 66: printing option number 100: max number of eigenvalues printed DOSS keyword 2 100 7 14 1 14 0 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 (first valence band) 14:    last band (fourth virtual band) 1:      plot option (if 1, the program stores the data in file fort.25); 14:    degree of the polynomial used for the DOSS expansion; 0:      printing option -1 1 projection onto all the AOs (-1) of Mg (1, first atom) -1 2 projection onto all the AOs (-1) of O (2, first atom)

1) Calculate the DOSS for MgO and process the data file fort.25, saved as mgo.f25 by runprop, with the command:

runprop   mgo_doss mgo                DOSS data saved in mgo_doss_mgo.f25
doss
mgo_doss_mgo                       DOSS plot postscript file mgo_doss_mgo.doss06.ps

 NEWK 6 6 1 1 66 100 DOSS 2 200 7 14 1 14 0 -1 1 -1 2 END 

2) The printed output,  mgo_doss_mgo.outp, looks like:

  ############################################################################### TOTAL AND PROJECTED DENSITY OF STATES - FOURIER LEGENDRE METHOD FROM BAND 7 TO BAND 14 ENERGY RANGE -0.11613E+01 0.11754E+01 ############################################################################### NUMBER OF LEGENDRE POLYNOMIALS 14 NUMBER OF SYMMETRIZED PWS FOR EXPANSION 16 NUMBER OF K POINTS OF SECONDARY NET 16 NUMBER OF PROJECTIONS 2 NUMBER OF ENERGY POINTS 198 *** PROJECTION 1 ATOMIC ORBITALS 1 2 3 4 5 6 7 8 9 *** PROJECTION 2 ATOMIC ORBITALS 10 11 12 13 14 15 16 17 18 BAND INTEGRATED DOSS PER PROJECTION AND TOTAL 1 2 T 7 0.00991 1.99009 2.00000 8 0.01277 1.98723 2.00000 9 0.01147 1.98853 2.00000 10 0.00491 1.99509 2.00000 11 0.47259 1.52741 2.00000 12 1.03565 0.96435 2.00000 13 1.24271 0.75729 2.00000 14 1.44140 0.55860 2.00000 

3) Visualize the file mgo_doss_mgo.doss06.ps:

Se Crgra2006 User manual to look at a different band energy range.

 Basic modeling of surfaces

Surfaces can be easily modelled with CRYSTAL by using the slab approach. This means that a thin-film (2D infinite system) perpendicular to the selected surface is cut out of the bulk structure. Note that the bulk structure must be carefully optimised before creating a surface.
Two aspects of general interest must be remarked:

1. Not all crystalline surfaces are physically stable or worthy of investigation
2. Need to have a slab thick enough to converge the property of interest

## Creating a slab model of a surface

All one needs to run a surface calculation with CRYSTAL is:

1. construct the bulk unit cell (see previous sections)
2. specify the Miller indices of the selected crystallographic face
3. specify the surface termination
4. specify the slab thickness

This is accomplished through the keyword SLABCUT. A new 2D unit cell is automatically found as well as the layer group operators. To gather all the information needed to the SLABCUT option, the keyword SLABINFO has to be used in a TEST calculation.

The keyword SLABINFO, followed (next record) by the crystallographic (Miller) indices of the surface, defines a new 3D cell, with two lattice vectors perpendicular to the [hkl] direction chosen.

Let us take MgO (rock salt structure) as a simple example and build a slab model of the (100) surface.

1) starting from the bulk structure optimized at the HF level, insert the SLABINFO keyword in the geometry input block and run a TEST calculation.
 MGO BULK - GEOMETRY TEST CRYSTAL 0 0 0 225 4.191 2 12 0.    0.    0. 8 0.5   0.5   0.5 MgO bulk structure SLABINFO 1 0 0 Keyword to perform a re-orientation of the unit cell perpendicular to a given (hkl) face TEST END Keyword to perform a test calculation End of the geometry input section

The atoms in the new reference cell are re-ordered according to their z coordinate, in order to recognize the layered structure, parallel to the (hkl) plane. The layers of atoms are numbered.
  ******************************************************************************* * CELL ROTATION ******************************************************************************* DIRECT-LATTICE FUNDAMENTAL VECTORS X Y Z A1 0.0000 2.0955 2.0955 A2 2.0955 0.0000 2.0955 A3 2.0955 2.0955 0.0000 ********************************************************************** * DEFINITION OF THE NEW LATTICE VECTORS * * TWO OF WHICH BELONG TO THE SELECTED PLANE * ********************************************************************** PLANE INDICES = 1 0 0 PRIMITIVE INDICES = 0 1 1 NEW FUNDAMENTAL DIRECT-LATTICE VECTORS B1,B2,B3 B1= 1 A1 0 A2 0 A3 B2= 0 A1 -1 A2 1 A3 B3= 1 A1 0 A2 -1 A3 LATTICE PARAMETERS (ANGSTROM AND DEGREES) A B C ALPHA BETA GAMMA 2.96348 2.96348 2.96348 120.00000 60.00000 90.00000 NUMBER OF PRIMITIVE CELLS CONTAINED IN THE NEW FUNDAMENTAL CELL : 1 *** THE CARTESIAN FRAME IS ROTATED SO THAT THE SELECTED PLANE IS IN THE X-Y PLANE NEW FUNDAMENTAL DIRECT-LATTICE VECTORS B1,B2,B3 X Y Z B1 2.9635 0.0000 0.0000 B2 0.0000 2.9635 0.0000 B3 1.4817 -1.4817 2.0955 VOLUME OF THE 3D CELL 18.403 AREA OF THE 2D CELL 8.782 B3 MODULUS 2.963 B3 Z PROJECTION 2.095 B3 X-Y PROJECTION 2.095 UNIT CELL ATOM COORDINATES : LAB AT.NO. CARTESIAN (ANG) CRYSTALLOGRAPHIC 1 12 0.000 0.000 0.000 0.000000 0.000000 0.000000 2 8 -1.482 1.482 0.000 -0.500000 0.500000 0.000000 1 12 0.000 0.000 0.000 0.000000 0.000000 0.000000 2 8 1.482 1.482 0.000 0.500000 0.500000 0.000000   ATOMS CLASSIFIED ACCORDING TO THE Z COORDINATE : LAYER 1 Z= 0.0000; LABEL AT.NO. X,Y,Z (ANG.) 1 12 0.00000 0.00000 0.00000 2 8 1.48174 1.48174 0.00000   *************************** CELL ROTATION COMPLETE *************************** 

In the case of MgO by specifying the (100) direction, a single atomic layer is singled out in which Mg and O are on the same plane. This allows one to check which the surface termination is and decide the number of atomic layers defining the slab model.
Note that for MgO, one atomic layer also represents the repeated unit that builds up the slab model.

2) with this information, the slab model can be created by specifying the SLABCUT keyword. For a 3-layer model of the (100) surface of MgO, the input looks like:
 MGO BULK - GEOMETRY TEST CRYSTAL 0 0 0 225 4.191 2 12 0.    0.    0. 8 0.5   0.5   0.5 MgO bulk structure SLABCUT 1 0 0 1 3 Keyword to cut out a slab from the bulk structure Miller indices of the selected face Surface termination and number of atomic layers TEST END Keyword to perform a test calculation End of the geometry input section

According to SLABINFO option, atomic layer 1 is chosen to be the surface termination and then 3 atomic layers are considered to create the slab model.
The corresponding output part follows
  ... ... LAYER 1 Z= 0.0000; LABEL AT.NO. X,Y,Z (ANG.) 1 12 0.00000 0.00000 0.00000 2 8 1.49093 1.49093 0.00000 *************************** CELL ROTATION COMPLETE ***************************   ******************************************************************************* * TWO DIMENSIONAL SLAB PARALLEL TO THE SELECTED PLANE ******************************************************************************* SURFACE LAYER NO. 1 TOTAL NUMBER OF LAYERS: 2 ORIGIN SHIFT TO MAXIMIZE USE OF SYMMETRY= -0.500 NUMBER OF SYMMOPS= 16 COORDINATES OF THE ATOMS BELONGING TO THE SLAB LAB AT.N. CARTESIAN (ANG) FRACTIONAL (3D) FRAC (2D) Z(ANG) 1 12 0.74547 -0.74547 1.05425 0.000 0.000 0.500 0.250 -0.250 1.05425 2 8 -0.74547 0.74547 1.05425 0.500 0.500 0.500 -0.250 0.250 1.05425 3 12 -0.74547 0.74547 -1.05425 0.000 0.000 -0.500 -0.250 0.250 -1.05425 4 8 0.74547 -0.74547 -1.05425 0.500 0.500 -0.500 0.250 -0.250 -1.05425 ******************************* SLAB GENERATED ******************************* 

A slab model of 3 atomic layer is then created from the 3D structure:

### Exercise 1:

Define a 3-layer slab model for the (110) face of MgO. Use the TEST keyword in the first block and then visualize the structure (solution)

## Computing the surface formation energy

The surface formation energy per unit area can be computed as

Ens = [E(n) - nEbulk] / 2A                    (1)

where E(n) is the energy of an n-layer slab per formula unit, Ebulk is the energy of the bulk structure per formula unit and A is the area of the primitive surface unit cell (8.782 Å2 as reported in the output). In this case, the energy of the bulk primitive unit cell has been obtained in the previous exercises to be: -274.6642392283 Hartree.

Thus, for the 3-layer model of the (100) surface of MgO, E3s = 3.331 mHa/ Å2 = 1.443 J/m2
Note: 1 Ha/ Å2 = 435.9748 J/m2.

As more layers are used in the calculation, Ens will converge to the surface energy if numerical errors in the calculation of E(n) and Ebulk cancel exactly. As CRYSTAL constructs the Hamiltonian matrix in direct space consistently in the SLAB (2D) and CRYSTAL (3D) calculations the main source of numerical differences between E(n) and Ebulk is due to k-space sampling.

### Exercise 2:

Compute E(n) and thus Ens for n=1, 2, 3, 4, 5, 6, 7 layers and confirm that for MgO(100) the surface energy converges to 1.444 Jm-2 (solutions).
 Layers SCF energy (hartree)

MgO is a wide band gap insulator and the convergence of Ebulk and E(n) with respect to k-space sampling is excellent. Ens converges rapidly and the surface energy is well defined. In metals or small band gap semiconductors this may not be the case. A more numerically stable definition of the surface energy is then:

Ens = {E(n) - n [E(n) - E(n-1)]} / 2A                    (2)

In this expression Ebulk is replaced with E(n)-E(n-1). If you consider each additional layer in the slab to be added as the central layer it is clear that E(n)-E(n-1) should converge to the energy of a single layer in the bulk crystal. However, in this case the bulk energy has effectively been computed under exactly the same numerical conditions as the surface layers leading to a cancellation of any systematic error.

### Exercise 2 bis:

Repeat the calculation of the surface formation energy for MgO(100) with eq. (2).

### Exercise 3:

Calculate the surface formation energy of a 3-layer slab model for the (110) face of MgO. Compare the result with E3s for MgO(100).
Which is the most stable surface? How does the surface energy converge with the slab thickness? Try 1, 2, 3, 4, 5, 6, 7 layers.

### Exercise 4:

Compare HF results with SVWN, PBE and B3LYP ones for a 3-layer slab model of the MgO(100) surface.
Remember that the bulk structure must be optimized before defining the slab model.

## Properties of the surface

Many of the interesting properties of surfaces (optical, chemical, electronic, magnetic) are related to the change in electronic structure which occurs in the near surface region. Here, we simply consider the Mulliken charges.

### Exercise 5

Check the convergence of the Mulliken charges with the slab thickness for MgO(100) and MgO(110). Compare the values for the two faces.
Compare the charges of Mg and O atoms at the surface layer and in the inner ones with the corresponding charges in the bulk structure.

Further information and more advanced exercises can be found in the Surface and Adsorption tutorial.

 Basic modelling of defects

Defects are created in crystals when the periodic structure is altered, for instance, by displacing, removing, adding or substituting one or more atoms. Such a local perturbation in the structure causes a modification of the properties of the system, which is often relevant.

When a point defect is created, the translational symmetry of the system is lost. Periodic models can nevertheless be applied to simulate a local defect by constructing a supercell, that is, a set of several unit cells inside which the defect is inserted.

In the 2-D example of generation of a supercell, the new lattice vectors length is 5 times the original unit cell ones. Translational symmetry is maintained, the area of the new cell is 25 times the original one.

A point defect is inserted at the centre of the cell: it could be anywhere, due to translational invariance. But this is the way to preserve as much symmetry as possible, because symmetry elements are conventionally referred to the origin of the reference system.
The distance between two point defects determines the entity of their interaction.

This basic tutorial presents the procedure followed to study a simple defect, Ca substitution in bulk MgO.
The first application of supercell model using the CRYSTAL code was published in J. Phys. Condens. Matter 5 (1993) 4793-4804.

No strong rearrangement in the electronic distribution of the crystal is expected after substitution of Ca for Mg.

The principal effect to be expected from the substitution of an Ca ion for Mg in bulk MgO is relaxation of the lattice, due to the different sizes of Mg2+ and Ca2+ cations (70 pm vs 100 pm). The experimental Mg-O distance is 210 pm, to be compared with 240 pm in CaO.

The substitutional defect formation energy refers to the ionic limits, and has to be evaluated according to the equation:

ECa = EMgO-Ca - N EMgO + EMg2+ - ECa2+

where:

ECa      defect formation energy
EMgO-Ca total energy per supercell of the defective structure
N
number of MgO units in the corresponding perfect lattice supercell
EMgO          energy of bulk MgO per MgO unit
EMg2+
total energy of isolated Mg2+ ion
ECa2+
total energy of isolated  Ca2+ ion

The defective structure is allowed to relax, in order to be in a minimum energy geometry.

Modelling of the system consists of two steps:

• definition of the perfect system, and choice of a supercell of suitable size;
• creation of the defect.

## Perfect system

### Defining the bulk structure

The space group of MgO is Fm3m , the asymmetric unit includes two atoms in special positions of multiplicity 1:

Mg (0,0,0)     O (1/2,1/2,1/2)

or exchanging the coordinates

O (0,0,0)     Mg (1/2,1/2,1/2)

The two choices are physically equivalent, but they are not equally convenient in our calculation, because of what is conventional in CRYSTAL. The defect consists of substitution of one Ca atom for Mg only: if the atom at the origin is substituted, the number of symmetry operators is not reduced.

The input to define the bulk structure geometry  of MgO is:
 MgO - bulk CRYSTAL 0 0 0 225 4.21 2 12  0.0 0.0 0.0  8  0.5 0.5 0.5 ENDG Title of the job Dimensionality of the system Crystallographic information (3D only) Space Group number Lattice parameter Number of non equivalent atoms Conventional atomic number,  fractionary coordinates End of the geometry input section

The geometry of the bulk MgO must be optimized, as shown in Geometry optimization.

### Supercell size

A supercell is defined by an expansion matrix, given in input, to multiply the lattice parameters matrix (see keyword SUPERCEL in CRYSTAL User's Manual). The expansion matrix to obtain a supercell of given volume must be chosen so as not to reduce the number of symmetry operators.
crystal checks the compatibility of the symmetry operators of the original cell with the new cell, and removes the symmetry operators not compatible with the structure. The number of symmetry operators can be reduced, not increased.

Techniques to define the expansion matrix are given in  "Defects"  tutorial. The supercell are named Snv, being nv the ratio between the volume of the supercell and the volume of the primitive one.

 Expansion matrix: from primitive S1 to crystallographic S4: SUPERCEL -1 1 1 1 -1 1 1 1 -1  Expansion matrix: from primitive S1 to S32: SUPERCEL -2 2 2 2 -2 2 2 2 -2 
 Expansion matrix: from primitive S1 to S8: SUPERCEL 2 0 0 0 2 0 0 0 2  Expansion matrix: from primitive S1 to S16: SUPERCEL -3 -1 -1 -1 3 -1 -1 -1 3 
 Expansion matrix: from primitive S1 to S64: SUPERCEL 4 0 0 0 4 0 0 0 4  Expansion matrix: from primitive S1 to S2: reduction of symmetry SUPERCEL 1 1 -1 0 0 1 1 -1 0 

The size of the reciprocal lattice cell reduces with increasing size of the direct lattice cell. The sampling in k space requires less and less points, at the limit 1 k point for very large supercell.
The Pack-Monkhorst net shrinking factor (see "SCF") ranges from 8 to 2.  The value of IS must be checked to guarantee stability of the energy value within the SCF convergence tolerance.

### Exercise 1

Prepare a geometry input (insert TEST) to generate the cells S1, S4 , S8, S16, S32 , S64 , and visualize them.

### Exercise 2

Determination of the best shrinking factor for S1 (2 atoms) S4 (8 atoms), S8 (16 atoms), S16 (32 atoms).
Look at bulk MgO input, to complete the geometry input of exercise #1 with basis set and SCF control

 S1 - 2 atoms - 18 AO - space group 225 - Fm-3m- 48 symmops Shrink n. k points n. SCF cycles cpu time (seconds) energy (hartree) 2 4 8 12

 S4 - 8 atoms - 72 AO - space group 221 - Pm-3m -  48 symmops Shrink n. k points n. SCF cycles cpu time (seconds) energy (hartree) 2 4 8 12

 S8 - 16 atoms - 144 AO - space group 225 - Fm-3m - 48 symmops Shrink n. k points n. SCF cycles cpu time (seconds) energy (hartree) 2 4 8 12

 S16 - 32 atoms - 288 AO - space group 229 - Im-3m -  48 symmops Shrink n. k points n. SCF cycles cpu time (seconds) energy (hartree) 2 4 8 12

 S32 - 64 atoms - 576 AO - space group 221 - Pm-3m - 48 symmops Shrink n. k points n. SCF cycles cpu time (seconds) energy (hartree) 2 4 8 12

 S64 - 128 atoms - 1152 AO - space group 225 - Fm-3m - 48 symmops Shrink n. k points n. SCF cycles cpu time (seconds) energy (hartree) 2 4 8 12

### Size-extensivity

The first essential requirement of supercell method is the size-extensivity. Each extensive property for an n-cells system must be equal to n times the corresponding value obtained for 1 cell system.

### Exercise 3

Looking at the outputs of exercise 2, verify that the energy per MgO pair is equal up to the SCF convergence threshold, for the supercells S8, S16 , S32, S64

 Type of cell Total energy  (hartree) Energy per formula unit (hartree) S1 (primitive cell) S8 S16 S32 S64

## Defective system

The program CRYSTAL allows creation of a substitutional defect using the keyword ATOMSUBS, inserted in geometry input, after the instructions to create the supercelll. The basis set for new atom must be defined in basis set input. From now on, work with S16

 ATOMSUBS 1 1 20 ATOMSYMM NEIGHBOR 11 keyword to substitute atoms number of atoms to be substituted sequence number of the atoms, new atomic number keyword to analyse local symmetry of atoms keyword to increase the range of neighbouring analysis number of neighbours to be printed

The Mg atom at the origin is substituted, to avoid reduction of symmetry.
(the keyword ATOMSYMM prints the symmetry relationship between atoms. The atom at the origin does not have any symmetry related). What is the distance between two Ca atoms? Check your answer looking at the neighbours analysis printout.

  N = NUMBER OF NEIGHBORS AT DISTANCE R ATOM N R/ANG R/AU NEIGHBORS (ATOM LABELS AND CELL INDICES) 1 CA 6 2.0955 3.9599 17 O 0 0 0 19 O 0 0 0 21 O 0 0 0 24 O 0 0 0 31 O 0 0 0 32 O 0 0 0 1 CA 12 2.9635 5.6002 2 MG 0 0 0 4 MG 0 0-1 5 MG 0 0 0 6 MG -1 0 0 7 MG 0 0 0 8 MG 0 0 0 10 MG 0 0 0 12 MG 0 0 0 13 MG 0-1 0 14 MG 0 0 0 15 MG 0 0 0 16 MG 0 0 0 1 CA 8 3.6295 6.8588 18 O 0 0 0 18 O 0 0-1 20 O -1-1-1 20 O 0 0 0 26 O 0 0 0 26 O -1 0 0 28 O 0 0 0 28 O 0-1 0 1 CA 6 4.1910 7.9198 3 MG 0 0 0 3 MG -1-1 0 9 MG 0 0 0 9 MG -1 0-1 11 MG 0 0 0 11 MG 0-1-1 1 CA 24 4.6857 8.8547 22 O 0 0 0 22 O 0-1 0 22 O 0 0 1 ............................................ 30 O 0 0 0 30 O 0-1-1 30 O 0-1 0 1 CA 24 5.1329 9.6998 2 MG -1-1-1 2 MG 0 0-1 4 MG 0 0 0 ............................................ 1 CA 12 5.9270 11.2003 3 MG -1-1-1 3 MG -1 0 0 3 MG 0-1 0 3 MG 0 0 1 9 MG -1-1-1 9 MG 0 0-1 9 MG -1 0 0 9 MG 0 1 0 11 MG -1-1-1 11 MG 0-1 0 11 MG 0 0-1 11 MG 1 0 0 1 CA 30 6.2865 11.8798 17 O 1 1 1 17 O 1 0 0 17 O 0 1 0 ............................................ 1 CA 24 6.6266 12.5224 2 MG -1 0-1 2 MG 0-1-1 4 MG 1 0 0 .......................................... 1 CA 24 6.9500 13.1336 18 O -1-1-1 18 O 1 0 0 18 O 0 1 0 ........................................... 1 CA 8 7.2590 13.7176 1 CA -1-1-1 1 CA 1 1 1 1 CA -1 0 0 1 CA 1 0 0 1 CA 0-1 0 1 CA 0 1 0 1 CA 0 0-1 1 CA 0 0 

  ATOM 1 ATOMIC NUMBER 20 - NUMBER OF SYMMOPS WHICH DO NOT MOVE THE ATOM 48 SEQUENCE NUMBERS OF THE SYMMOPS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 NO EQUIVALENT ATOMS ................................................................ ATOM 17 ATOMIC NUMBER 8 - NUMBER OF SYMMOPS WHICH DO NOT MOVE THE ATOM 8 SEQUENCE NUMBERS OF THE SYMMOPS 1 2 15 16 27 28 37 38 NUMBER OF ATOMS EQUIVALENT BY SYMMETRY 5 SEQUENCE NUMBERS OF THESE ATOMS 19 31 24 21 32 SEQUENCE NUMBERS OF THE CORRESPONDING SYMMOPS 3 5 6 7 10 ......................................................... ATOM 2 ATOMIC NUMBER 12 - NUMBER OF SYMMOPS WHICH DO NOT MOVE THE ATOM 4 SEQUENCE NUMBERS OF THE SYMMOPS 1 14 26 37 NUMBER OF ATOMS EQUIVALENT BY SYMMETRY 11 SEQUENCE NUMBERS OF THESE ATOMS 4 10 12 16 5 8 13 6 15 14 7 SEQUENCE NUMBERS OF THE CORRESPONDING SYMMOPS 2 3 4 5 6 7 8 9 10 11 12 

### Exercise 4

Compute the formation energy of the defect.

Relaxation of the structure

The cell parameter (the cell is cubic, there is only 1 degree of freedom) is fixed, corresponding to the one optimized for the perfect bulk structure.

Looking at the atom symmetry and the neighbour analysis, we see that there is a star of 6 Oxygen symmetry related at a distance of 2.0955 Angstrom around the defect, then  a star of 12 Mg at a distance of 2.9635 Angstrom, and out of it a star of 8 Oxygen at 3.6295 Angstrom. The atoms in each star are symmetry related, relaxation of the structure modifies the "radius" of each star.

 ATOMSUBS 1 1 20 ATOMSYMM NEIGHBOR 11 OPTGEOM END keyword to substitute atoms number of atoms to be substituted sequence number of the atoms, new atomic number keyword to analyse local symmetry of atoms keyword to increase the range of neighbouring analysis number of neighbours to be printed optimization - default atom coordinates only

 Unrelaxed Relaxed Energy (hartree) first star (Ox) (A) second star (Mg) (A) third star (Ox) (A) fourth star (A)

The SCF guess is a superposition of atomic densities. The electronic configuration given in input assigns 10 electron to both Mg and Ox, and  18 to Ca. The energy for the equivalent ions is printed at the beginning of SCF process:

 Ion Total energy (hartree) Ca2+ Mg2+ O2-

### Exercise 5

Verify the assumption "No strong rearrangement in the electronic distribution of the crystal is expected after substitution of Mg with Ca" by looking at the net atomic charges of the defective system, to be compared with the perfect MgO bulk ones.

 Atom Net atomic charges Ca O1 Mg1 O2 Mg2 Mg perfect bulk O perfect bulk