CRYSTAL home page CRYSTAL home page

MSSC2007 Logo

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.

Background requirements:

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?
    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.

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

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:

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:

 

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:

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:

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

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

MgO electron density: total and difference maps (linear scale).
\begin{figure}{\psfig {figure=tabelle/mgodiff.eps,height=16cm}}

\protect\end{figure}

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
    

    Si electron density: total and difference maps (linear scale).
    \begin{figure}{\psfig {figure=tabelle/Sibulk.maps.pos,height=16cm}}

\protect\end{figure}

    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
Urea electron density: total and difference maps (logarithmic scale).
\begin{figure}{\psfig {figure=tabelle/atomorde.maps.pos,height=16cm}}

\protect\end{figure}

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) 
The Brillouin zone for the face centered cubic lattice.

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: $\Gamma \rightarrow$ X
6 0 6    6 3 9   X $\rightarrow$ W
6 3 9    6 6 6   W $\rightarrow$ L
6 6 6    0 0 0   L $\rightarrow \Gamma$
0 0 0    6 3 9   $\Gamma \rightarrow$

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