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.
Background requirements:
A review article is devoted to a basic introduction to CRYSTAL as a tool for quantum simulation in solid state chemistry :
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 1D2D3D systems only) optional general information and SCF related keywords END 
Examples of complete CRYSTAL inputs are reported in the following table:
Molecule0D  Polymer1D  Slab2D  Crystal3D  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: 
Basically, a crystal structure is determined by:
1) space group;
2) cell parameters;
3)
coordinates of the atoms in the asymmetric unit.
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 Xray determination of electrondensity 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) 4348 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 (HermannMauguin symbol and Space group number) 
CLAS
m3m (HermannMauguin)  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 431022 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 
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.000000000000E01 5.000000000000E01 5.000000000000E01 TRANSFORMATION MATRIX PRIMITIVECRYSTALLOGRAPHIC 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.000000000000E01 5.000000000000E01 5.000000000000E01 
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 nonprimitive, 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:


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 
P42_{1}m, 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 A^{3} 
Density  1.37 g/cm^{3} 
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.00000000000E01 3.26000000000E01 2 1 2 6 C 5.00000000000E01 0.00000000000E+00 3.26000000000E01 3 2 1 8 O 0.00000000000E+00 5.00000000000E01 4.04700000000E01 4 2 2 8 O 5.00000000000E01 0.00000000000E+00 4.04700000000E01 5 3 1 7 N 1.45900000000E01 3.54100000000E01 1.76600000000E01 6 3 2 7 N 1.45900000000E01 3.54100000000E01 1.76600000000E01 7 3 3 7 N 3.54100000000E01 1.45900000000E01 1.76600000000E01 8 3 4 7 N 3.54100000000E01 1.45900000000E01 1.76600000000E01 9 4 1 1 H 2.57500000000E01 2.42500000000E01 2.82700000000E01 10 4 2 1 H 2.57500000000E01 2.42500000000E01 2.82700000000E01 11 4 3 1 H 2.42500000000E01 2.57500000000E01 2.82700000000E01 12 4 4 1 H 2.42500000000E01 2.57500000000E01 2.82700000000E01 13 5 1 1 H 1.44100000000E01 3.55900000000E01 3.80000000000E02 14 5 2 1 H 1.44100000000E01 3.55900000000E01 3.80000000000E02 15 5 3 1 H 3.55900000000E01 1.44100000000E01 3.80000000000E02 16 5 4 1 H 3.55900000000E01 1.44100000000E01 3.80000000000E02 
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?
Your output should look like:
******************************************************************************* * 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.
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).
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.
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 (s^{2} p^{6}) 
2 
p  3  x, y, z  6 
3 
d  5  2z^{2}x^{2}y^{2}, xz, yz, x^{2}y^{2}, xy  10 
4 
f  7  (2z^{2}3x^{2}3y^{2})z, (4z^{2}x^{2}y^{2})x, (4z^{2}x^{2}y^{2})y, (x^{2}y^{2})z, xyz, (x^{2}3y^{2})x, (3 x^{2}y^{2})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 n_{s }of the atomic basis set
for
each shell (n_{s }blocks of records): type of basis
set (012), type of shell (01234),
number of primitives GTF n_{g}, shell electronic charge, scale factor
for each primitive (n_{g} 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:
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 STOnG 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:
STO3G 
621G modified  321G 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 
The definition of the basis set for all the atoms ends with the records (mandatory):
99 0 END 
1) Write the input for silicon (geometry + basis set) with a basis set taken from the CRYSTAL website (Si_8831G*_nada_1996).
2) Run a TEST based on this input.
3) How many functions (AO) do you get with the proposed basis set?
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 (1s^{2}, (2sp)^{8}), 10 electrons to Mg (1s^{2}, (2sp)^{8}); the SCF guess is a superposition of ionic densities, Mg^{2+} O^{2 }To obtain a superposition of atomic densities, the shell charges should be: O (1s^{2}, (2sp)^{6}) and Mg (1s^{2}, (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.
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?
The hamiltonian specification is in the third input block. The hamiltonians implemented in CRYSTAL are HF (HartreeFock) 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 
The keyword DFT selects a DFT Hamiltonian. Exchangecorrelation functionals are separated in an exchange component (keyword EXCHANGE) and a correlation component (keyword CORRELAT). If the correlation functional is not set, an exchangeonly functional is used in the hamiltonian. If the exchange functional is not set, the HartreeFock 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:
See CRYSTAL User's Manual for a complete list of available exchangecorrelation 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.
CRYSTAL can run in different ways:
1) traditional SCF: default, mono and twoelectron integrals are stored on disk and readin at each SCF iteration;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 bielectron integrals to be stored is printed at the bottom of a testrun output (1 word = 8 byte).
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 PackMonkhorst 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:
== 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.
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 3dimensional 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).
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.000000000000E01 5.000000000000E01 5.000000000000E01

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

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 C_{2v} geometry of molecular urea in the bulk, the most stable
structure has a C_{2} symmetry with an anti configuration (see figure below).
Therefore, when computing the interaction energy, E(conf) between the C_{2} and
the C_{2v} 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  C_{2v})
E(conf) = E(mol  C_{2v}) 
E(mol  C_{2})
BSSE = E(mol  bulk)  E(mol  ghosts)
because there are two molecules in the cell.
In summary, one needs:
All the requested calculations will be performed with the HF Hamiltonian and with 321G basis sets for all the atoms involved.
System 
Label  SCF energy (hartree) 
E(bulk)  E_{1}  
E(mol  C2v)  E_{2}  
E(mol  C2)  E_{3}  
E(mol  bulk)  E_{4}  
E(mol  ghosts)  E_{5} 
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  SINGLEPOINT MOLECULE IN THE BULK CRYSTAL 0 0 0 113 5.59037346 4.67007538 5 6 0.000000000000E+00 5.000000000000E01 3.176227372521E01 8 0.000000000000E+00 5.000000000000E01 4.093385097982E01 7 1.437784027012E01 3.562215972988E01 1.678775482565E01 1 2.564186903752E01 2.435813096248E01 2.695980532036E01 1 1.406087223437E01 3.593912776563E01 4.827103388750E02 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.000000000000E01 3.176227372521E01 8 0.000000000000E+00 5.000000000000E01 4.093385097982E01 7 1.437784027012E01 3.562215972988E01 1.678775482565E01 1 2.564186903752E01 2.435813096248E01 2.695980532036E01 1 1.406087223437E01 3.593912776563E01 4.827103388750E02 MOLEBSSE 1 1 0 0 0 1 4.00 END 
Geometry block ex. 4) UREA MOLECULE C2v MOLECULE 15 5 6 4.010040108509E18 4.010040108509E18 1.511178418937E01 8 3.627793330734E17 3.627793330734E17 1.367129216024E+00 7 0.0000000000000000 1.152517719977E+00 6.029124102443E01 1 0.0000000000000000 1.186574472949E+00 1.524738700435E+00 1 0.0000000000000000 2.002684092087E+00 6.757791827922E02 OPTGEOM END END 
Geometry block ex. 5) UREA MOLECULE C2 MOLECULE 5 5 6 4.010040108509E18 4.010040108509E18 1.511178418937E01 8 3.627793330734E17 3.627793330734E17 1.367129216024E+00 7 1.475136395636E01 1.152517719977E+00 6.029124102443E01 1 2.591127703870E01 1.186574472949E+00 1.524738700435E+00 1 6.956803083756E02 2.002684092087E+00 6.757791827922E02 OPTGEOM END END 
Exercise: Recalculate the cohesive energy of crystalline urea using the following data obtained with the 631G(d,p) basis set:
System 
Label  SCF energy (hartree) 
E(bulk)  E_{1}  4.4804947630236E+02 
E(mol  C2v)  E_{2}  2.2398536052639E+02 
E(mol  C2)  E_{3}  2.2399715236476E+02 
E(mol  bulk)  E_{4}  2.2399397220548E+02 
E(mol  ghosts)  E_{5}  2.2400050884526E+02 
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 firstderivative 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).
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.4659E04 2.746642742333E+02 6 1.2600E06 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.2379E21 0.0000 0.0000 (Fu ) A ( 0.00) I 4 6 0.4503E05 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 = ZEROPOINT 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+PVTS : 0.000315860975 0.008595014089 0.82929287 EL+E0+ET+PVTS: 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 
1) Analyze the frequency output of calcite (CaCO_{3}):
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.6415E04 1.883021621785E+03 7 9.9308E07 1 1 CA DY GENERATED FROM A PREVIOUS LINE 1 CA DZ 5.3227E04 1.883021621201E+03 7 1.5762E06 3 3 C DX 4.0721E03 1.883021611246E+03 10 1.1531E05 1 3 C DY GENERATED FROM A PREVIOUS LINE 3 C DZ 1.5595E03 1.883021618367E+03 9 4.4102E06 3 5 O DX 2.4149E03 1.883021615750E+03 9 7.0277E06 1 5 O DY GENERATED FROM A PREVIOUS LINE 5 O DZ 5.1208E04 1.883021621726E+03 9 1.0521E06 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.2676E21 0.0000 0.0000 (Eu ) A ( 0.00) I 3 3 0.3271E21 0.0000 0.0000 (A2u) A ( 0.00) I 4 5 0.3172E06 123.6060 3.7056 (Eu ) A ( 144.60) I 6 6 0.3602E06 131.7273 3.9491 (A2u) A ( 104.21) I 7 8 0.4969E06 154.7137 4.6382 (Eg ) I ( 0.00) A 9 9 0.8120E06 197.7660 5.9289 (A2g) I ( 0.00) I 10 11 0.1004E05 219.8920 6.5922 (Eu ) A ( 189.99) I 12 13 0.1613E05 278.7694 8.3573 (Eg ) I ( 0.00) A 14 15 0.1676E05 284.1535 8.5187 (Eu ) A ( 812.65) I 16 16 0.1744E05 289.8341 8.6890 (A1u) I ( 0.00) I 17 17 0.1842E05 297.8520 8.9294 (A2u) A ( 347.09) I 18 18 0.2028E05 312.5693 9.3706 (A2g) I ( 0.00) I 19 20 0.1036E04 706.3442 21.1757 (Eg ) I ( 0.00) A 21 22 0.1037E04 706.7574 21.1881 (Eu ) A ( 39.53) I 23 23 0.1564E04 867.9843 26.0215 (A2u) A ( 205.01) I 24 24 0.1627E04 885.2715 26.5398 (A2g) I ( 0.00) I 25 25 0.2442E04 1084.5219 32.5131 (A1g) I ( 0.00) A 26 26 0.2444E04 1084.9299 32.5254 (A1u) I ( 0.00) I 27 28 0.4216E04 1425.1279 42.7243 (Eu ) A ( 5579.85) I 29 30 0.4408E04 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 = ZEROPOINT 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+PVTS : 0.006745352377 0.183550369687 17.70992017 EL+E0+ET+PVTS: 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 
Oneelectron 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.
Oneelectron 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:
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 oneelectron properties from the wave function:
A complete list of properties is presented in http://www/crystal.unito.it/features.html
To run a property 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
A Mulliken population analysis can be performed by inserting the keyword PPAN. For each nonequivalent 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 
and your output is
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 XH.....Y distance.
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 BA 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 Mg^{2+} O^{2})
See Crgra2006 user manual to draw difference maps.
The total density and difference density maps look like:
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:
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 
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 XHY.
The solutions to this exercise are the following:
ECHG 0 65 ATOMS 2 0 0 1 2 0 0 0 2 1 1 0 RECTANGU MARGINS 2. 2. 2. 2. PPAN END 
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).
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.
The solutions to this exercise are the following:
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 
The calculation of density of states requires the Fock/KS eigenvectors in the k points defined by PackMonkhorst 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:
In case of a projected density onto a subset of AOs, some other data must be specified in the input file:
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
PackMonkhorst 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
Your input is
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 thinfilm (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:
All one needs to run a surface calculation with CRYSTAL is:
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 reorientation 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 reordered 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 ******************************************************************************* DIRECTLATTICE 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 DIRECTLATTICE 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 XY PLANE NEW FUNDAMENTAL DIRECTLATTICE 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 XY 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.
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:
The surface formation energy per unit area can be computed as
E_{ns} = [E(n)  nE_{bulk}] / 2A
(1)
where E(n) is the energy of an nlayer slab per formula unit, E_{bulk}
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 3layer model of the (100) surface of MgO, E_{3s} = 3.331 mHa/ Å^{2}
= 1.443 J/m^{2}
Note: 1 Ha/ Å^{2} = 435.9748 J/m^{2}.
As more layers are used in the calculation, E_{ns} will converge to the surface energy
if numerical errors in the calculation of E(n) and E_{bulk} 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 E_{bulk} is due to kspace sampling.
Layers 
SCF energy (hartree) 
MgO is a wide band gap insulator and the convergence of E_{bulk} and E(n) with respect
to kspace sampling is excellent.
E_{ns} 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:
E_{ns} = {E(n)  n [E(n)  E(n1)]} / 2A
(2)
In this expression E_{bulk} is replaced with E(n)E(n1). If you consider each additional
layer in the slab to be added as the central layer it is clear that E(n)E(n1) 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.
Repeat the calculation of the surface formation energy for MgO(100) with eq. (2).
Calculate the surface formation energy of a 3layer slab model for the (110) face of MgO.
Compare the result with E_{3s} 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.
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 2D 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) 47934804.
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 Mg^{2+} and Ca^{2+} cations (70 pm vs 100 pm). The experimental MgO 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:
E_{Ca} = E_{MgOCa}  N E_{MgO} + E_{Mg2+}  E_{Ca2+}
where:
E_{Ca }
defect formation energy
E_{MgOCa
}
total energy per supercell of the defective structure_{
N } number
of MgO units in the corresponding perfect lattice supercell_{
}
E_{MgO } energy
of bulk MgO per MgO unit_{}
E_{Mg2+}_{ } total
energy of isolated Mg^{2+} ion_{
} E_{Ca2}_{+ } total
energy of isolated Ca^{2+} 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:
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.
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 S_{nv}, being nv the ratio between the volume of the supercell and the volume of the primitive one.
Expansion matrix: from primitive S_{1} to crystallographic S_{4}: SUPERCEL 1 1 1 1 1 1 1 1 1 
Expansion matrix: from primitive S_{1} to S_{32}: SUPERCEL 2 2 2 2 2 2 2 2 2 
Expansion matrix: from primitive S_{1} to S_{8}: SUPERCEL 2 0 0 0 2 0 0 0 2 
Expansion matrix: from primitive S_{1} to S_{16}: SUPERCEL 3 1 1 1 3 1 1 1 3 
Expansion matrix: from primitive S_{1} to S_{64}: SUPERCEL 4 0 0 0 4 0 0 0 4 
Expansion matrix: from primitive S_{1} to S_{2}: 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 PackMonkhorst 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.
S_{1}  2 atoms  18 AO  space group 225  Fm3m 48 symmops  
Shrink  n. k points  n. SCF cycles  cpu time (seconds)  energy (hartree) 
2  
4  
8  
12 
S_{4}  8 atoms  72 AO  space group 221  Pm3m  48 symmops  
Shrink  n. k points  n. SCF cycles  cpu time (seconds)  energy (hartree) 
2  
4  
8  
12 
S_{8}  16 atoms  144 AO  space group 225  Fm3m  48 symmops  
Shrink  n. k points  n. SCF cycles  cpu time (seconds)  energy (hartree) 
2  
4  
8  
12 
S_{16}  32 atoms  288 AO  space group 229  Im3m  48 symmops  
Shrink  n. k points  n. SCF cycles  cpu time (seconds)  energy (hartree) 
2  
4  
8  
12 
S_{32}  64 atoms  576 AO  space group 221  Pm3m  48 symmops  
Shrink  n. k points  n. SCF cycles  cpu time (seconds)  energy (hartree) 
2  
4  
8  
12 
S_{64}  128 atoms  1152 AO  space group 225  Fm3m  48 symmops  
Shrink  n. k points  n. SCF cycles  cpu time (seconds)  energy (hartree) 
2  
4  
8  
12 
The first essential requirement of supercell method is the sizeextensivity. Each extensive property for an ncells system must be equal to n times the corresponding value obtained for 1 cell system.
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 S_{8}, S_{16 }, S_{32}, S_{64}
Type of cell  Total energy (hartree)  Energy per formula unit (hartree) 
S_{1} (primitive cell)  
S_{8}  
S_{16 }  
S_{32}  
S_{64} 
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 S_{16 }
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 01 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 01 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 01 20 O 111 20 O 0 0 0 26 O 0 0 0 26 O 1 0 0 28 O 0 0 0 28 O 01 0 1 CA 6 4.1910 7.9198 3 MG 0 0 0 3 MG 11 0 9 MG 0 0 0 9 MG 1 01 11 MG 0 0 0 11 MG 011 1 CA 24 4.6857 8.8547 22 O 0 0 0 22 O 01 0 22 O 0 0 1 ............................................ 30 O 0 0 0 30 O 011 30 O 01 0 1 CA 24 5.1329 9.6998 2 MG 111 2 MG 0 01 4 MG 0 0 0 ............................................ 1 CA 12 5.9270 11.2003 3 MG 111 3 MG 1 0 0 3 MG 01 0 3 MG 0 0 1 9 MG 111 9 MG 0 01 9 MG 1 0 0 9 MG 0 1 0 11 MG 111 11 MG 01 0 11 MG 0 01 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 01 2 MG 011 4 MG 1 0 0 .......................................... 1 CA 24 6.9500 13.1336 18 O 111 18 O 1 0 0 18 O 0 1 0 ........................................... 1 CA 8 7.2590 13.7176 1 CA 111 1 CA 1 1 1 1 CA 1 0 0 1 CA 1 0 0 1 CA 01 0 1 CA 0 1 0 1 CA 0 01 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 
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) 
Ca^{2+}  
Mg^{2+}  
O^{2} 
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 