CRYSTAL2003
1.0
User's Manual
August 24, 2004
Main authors of CRYSTAL 2003:
V.R. Saunders1, R. Dovesi2, C. Roetti2, R. Orlando3,2,
C. M. Zicovich-Wilson 2,4, N.M. Harrison1,5, K. Doll1,6,
B. Civalleri2, I.J. Bush1, Ph. D'Arco7,2, M. Llunell2,8
1
Computational Science & Engineering Department - CLRC Daresbury
Daresbury, Warrington, Cheshire, UK WA4 4AD
2
Theoretical Chemistry Group - University of Turin
Department of Chemistry IFM
Via Giuria 5 - I 10125 Torino - Italy
3
University of Eastern Piedmont
Department of Science and Advanced Technologies
Corso Borsalino 54 - I 15100 Alessandria - Italy
4
Departamento de Fisica, Universidad Autonoma del Estado de Morelos,
Av. Universidad 1001, Col. Chamilpa, 62210 Cuernavaca (Morelos) Mexico
5
Chemistry, Imperial College
South Kensington Campus
London, U.K.
6
Institut f¨
ur Mathematische Physik
TU Braunschweig - Mendelssohnstrasse 3
Braunschweig D-38106 Germany
7
Laboratoire de P´
etrologie, Mod´
elisation des Materiaux et Processus
Universit´
e Pierre et Marie Curie,
4 Place Jussieu, 75232 Paris CEDEX 05, France
8
Departament de Quimica Fisica, Universitat de Barcelona
Diagonal 647, Barcelona, Spain
1
2
Contents
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
Functionality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
Typographical Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
Getting Started - Installation and testing . . . . . . . . . . . . . . . . . . . . . . . .
10
Program errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
1
Wave function calculation
Basic input route
11
1.1
Geometry and symmetry information . . . . . . . . . . . . . . . . . . . . . . . .
11
Geometry input for crystalline compounds . . . . . . . . . . . . . . . . . . . . .
12
Geometry input for molecules, polymers and slabs
. . . . . . . . . . . . . . . .
12
Geometry input from external geometry editor . . . . . . . . . . . . . . . . . .
13
Comments on geometry input . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
1.2
Basis set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
1.3
General information, computational parameters,
hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
1.4
SCF input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
2
Wave function calculation - Advanced input route
22
2.1
Geometry editing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
Finite periodic electric field perturbation . . . . . . . . . . . . . . . . . . . . . .
41
Geometry optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
Notes - From periodic structure to clusters and molecules . . . . . . . . . . . .
47
Notes - The slab model
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
Notes - BSSE correction in periodic systems . . . . . . . . . . . . . . . . . . . .
47
2.2
Basis set input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
Effective core pseudo-potentials . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
Pseudopotential libraries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
Valence Basis set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
Conversion of Stuttgart-Dresden pseudopotentials . . . . . . . . . . . . . . . . .
53
Conversion of Stevens et al. pseudopotentials . . . . . . . . . . . . . . . . . . .
54
2.3
General information, computational parameters,
hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
2.4
SCF input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
3
Properties
82
3.1
Preliminary calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
3.2
Properties keywords . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
3.3
Spontaneous polarization and piezoelectricity . . . . . . . . . . . . . . . . . . . 112
4
Input examples
115
4.1
Standard geometry input
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
CRYSTAL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
SLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
POLYMER . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
3
MOLECULE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.2
Basis set input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Valence only basis set input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.3
SCF options - SPINEDIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
4.4
Geometry optimization - OPTCOORD . . . . . . . . . . . . . . . . . . . . . . . 127
5
Basis set
133
5.1
Molecular BSs performance in periodic systems . . . . . . . . . . . . . . . . . . 133
5.2
Core functions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
5.3
Valence functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Molecular crystals
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Covalent crystals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Ionic crystals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
From covalent to ionics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
Metals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
5.4
Hints on crystalline basis set optimization . . . . . . . . . . . . . . . . . . . . . 136
5.5
Check on basis-set quasi-linear-dependence
. . . . . . . . . . . . . . . . . . . . 137
6
Theoretical framework
139
6.1
Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
6.2
Remarks on the evaluation of the integrals . . . . . . . . . . . . . . . . . . . . . 140
6.3
Treatment of the Coulomb series . . . . . . . . . . . . . . . . . . . . . . . . . . 141
6.4
The exchange series
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
6.5
Bipolar expansion approximation of Coulomb and exchange integrals . . . . . . 143
6.6
Exploitation of symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
Symmetry-adapted Crystalline Orbitals
. . . . . . . . . . . . . . . . . . . . . . 144
6.7
Reciprocal space integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
6.8
Electron momentum density and related quantities . . . . . . . . . . . . . . . . 146
6.9
Elastic Moduli of Periodic Systems . . . . . . . . . . . . . . . . . . . . . . . . . 147
Examples of
matrices for cubic systems . . . . . . . . . . . . . . . . . . . . . . 149
Bulk modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
6.10 Spontaneous polarization through the Berry phase approach . . . . . . . . . . . 152
Spontaneous polarization through the localized crystalline orbitals approach . . 152
6.11 Piezoelectricity through the Berry phase approach . . . . . . . . . . . . . . . . 153
Piezoelectricity through the localized crystalline orbitals approach
. . . . . . . 153
A Symmetry groups
155
A.1 Labels and symbols of the space groups . . . . . . . . . . . . . . . . . . . . . . 155
A.2 Labels of the layer groups (slabs) . . . . . . . . . . . . . . . . . . . . . . . . . . 158
A.3 Labels of the rod groups (polymers) . . . . . . . . . . . . . . . . . . . . . . . . 159
A.4 Labels of the point groups (molecules) . . . . . . . . . . . . . . . . . . . . . . . 162
A.5 From conventional to primitive cells: transforming matrices . . . . . . . . . . . 163
B Summary of input keywords
164
C DFT integration through an auxiliary basis set fitting
172
C.1 DFT input example - fitting method . . . . . . . . . . . . . . . . . . . . . . . . 173
D Reciprocal lattice sampling
176
E Printing options
177
F External format
181
G Utility programs
182
H CRYSTAL2003 versus CRYSTAL98
183
4
I
Relevant strings
186
J
Acronyms
187
Bibliography
188
Subject index
194
5
Introduction
The CRYSTAL package performs ab initio calculations of the ground state energy, energy
gradient, electronic wave function and properties of periodic systems. Hartree-Fock or Kohn-
Sham Hamiltonians (that adopt an Exchange- Correlation potential following the postulates of
Density-Functional theory) can be used. Systems periodic in 0 (molecules, 0D), 1 (polymers,
1D), 2 (slabs, 2D), and 3 dimensions (crystals, 3D) are treated on an equal footing. In each
case the fundamental approximation made is the expansion of the single particle wave functions
('Crystalline Orbital', CO) as a linear combination of Bloch functions (BF) defined in terms
of local functions (hereafter indicated as 'Atomic Orbitals', AOs). See Chapter 6.
The local functions are, in turn, linear combinations of Gaussian type functions (GTF) whose
exponents and coefficients are defined by input (section 1.2). Functions of symmetry s, p, d
and f can be used (see page 18. Also available are sp shells (s and p shells, sharing the same
set of exponents). The use of sp shells can give rise to considerable savings in CPU time.
The program can automatically handle space symmetry: 230 space groups, 80 layer groups,
99 rod groups, 45 point groups are available (Appendix A). In the case of polymers it cannot
treat helical structures (translation followed by a rotation around the periodic axis). However,
when commensurate rotations are involved, a suitably large unit cell can be adopted.
Point symmetries compatible with translation symmetry are provided for molecules.
Input tools allow the generation of slabs (2D system) or clusters (0D system) from a 3D
crystalline structure, the elastic distortion of the lattice, the creation of a supercell with a
defect and a large variety of structure editing. See Section 2.1
Previous releases of the software in 1988 (CRYSTAL88, [1]), 1992 (CRYSTAL92, [2]), 1996
(CRYSTAL95, [3]) and 1998 (CRYSTAL98, [4]), have been used in a wide variety of research
with notable applications in studies of stability of minerals, oxide surface chemistry, and defects
in ionic materials.
The CRYSTAL package has been developed over a number of years. Many of the algorithms
developed and reviews of the applications of the code have been previously published (for
instance, [5, 6, 7, 8, 9, 10, 11, 12, 13]).
The required citation for this work is:
V.R. Saunders, R. Dovesi, C. Roetti, R. Orlando, C. M. Zicovich-Wilson,
N.M. Harrison, K. Doll, B. Civalleri, I. Bush, Ph. D'Arco, M. Llunell
CRYSTAL2003 User's Manual, University of Torino, Torino, 2003
CRYSTAL2003 output will display the references relevant to the property computed, when
necessary.
The following sites will supply updated information on the CRYSTAL code:
http://www.cse.dl.ac.uk/Activity/CRYSTAL
http://www.crystal.unito.it
6
Functionality
The basic functionality of the code is outlined below.
The single particle potential
Restricted Hartree-Fock Theory
Unrestricted Open Shell Hartree-Fock Theory
Density Functional Theory for Exchange and Correlation
Spin Density Functional Theory
Hybryds HF-DFT (B3LYP-B3PW)
Effective Core Pseudopotentials
Finite field perturbation added to the Hamiltonian
*
Algorithms
Parallel processing (replicated data)
Massive Parallel Processing (distributed data)
*
Traditional SCF
Full Direct SCF
*
Structural Editing
Use of space, layer, rod and point group symmetry
Deformation of the crystallographic cell
Removal, insertion, deletion and substitution of atoms
Displacement of atoms
Rotation of groups of atoms
Extraction of surface models from a 3D crystal structure
Cluster generation from a 3D crystal
Cluster of molecules from molecular crystal
Properties
Band structure
Density of states
Electronic charge density maps
Electronic charge density on a 3D grid
Mulliken population analysis
Spherical harmonic atom and shell multipoles
X-ray structure factors
Electron momentum distributions
Compton profiles
First order density matrix
*
Reciprocal form factors
*
Electrostatic potential, field and field gradients
Spin polarized generalization of properties
Hyperfine electron-nuclear spin tensor
A posteriori Density Functional correlation energy
Localization of Crystalline Orbitals
*
Spontaneous polarization through Berry phase approach
Spontaneous polarization through localized orbitals approach *
Piezoelectricity through Berry phase approach
*
Piezoelectricity through localized orbitals approach
*
Optical dielectric constant
*
Analytic nuclear coordinates gradient of the energy
*
Geometry optimizer
*
7
Conventions
In the description of the input data which follows, the following notation is adopted:
-
·
new record
-
free format record
-
An
alphanumeric datum (first n characters meaningful)
-
atom label
sequence number of a given atom in the primitive cell, as
printed in the output file after reading of the geometry input
-
symmops
symmetry operators
-
, [ ]
default values.
-
italic
optional input
-
optional input records follow
II
-
additional input records follow
II
Part of the code is written in fortran 77. The name of the variables is associated with the type
of data, following the fortran 77 convention: if the first letter of the name is I, J, K, L, M or
N, the type is integer. Otherwise the type is real.
Arrays are read in with a simplified implied DO loop instruction of Fortran 77:
(dslist, i=m1,m2)
where: dslist is an input list; i is the name of an integer variable, whose value ranges from m1
to m2.
Example (page 25): LB(L),L=1,NL
NL integer data are read in and stored in the first NL position of the array LB.
All the keywords are entered with an A format (case insensitive); the keywords must be typed
left-justified, with no leading blanks.
conventional atomic number (usually callet NAT) is used to associate a given basis set
with an atom. The real atomic number is the remainder of the division NAT/100.
8
Acknowledgements
Embodied in the present code are elements of programs distributed by other groups.
In particular: the atomic SCF package of Roos et al. [14], the GAUSS70 gaussian integral
package and STO-nG basis set due to Hehre et al. [15], the code of Burzlaff and Hountas for
space group analysis [16] and Saunders' ATMOL gaussian integral package [17].
We take this opportunity to thank these authors. Our modifications of their programs have
sometimes been considerable. Responsibility for any erroneous use of these programs therefore
remains with the present authors.
It is our pleasure to thank Piero Ugliengo for continous help, useful suggestions, rigorous test-
ing.
We are also in debt with Cesare Pisani, for his constant support of and interest in the devel-
opment of the new version of the CRYSTAL program.
Financial support for this research has been provided by the italian MURST (Ministero della
Universit`
a e della Ricerca Scientifica e Tecnologica), and the United Kingdom CCLRC (Council
for the Central Laboratories of the Research Council).
9
Getting Started
See http://www.crystal.unito.it documentation README.
Program errors
A very large number of tests have been performed by researchers of a few laboratories, that
had access to a test copy of CRYSTAL2003. We tried to check as many options as possible,
but not all the possible combinations of options have been checked. We have no doubts that
errors remain.
The authors would greatly appreciate comments, suggestions and criticisms by the users of
CRYSTAL; in case of errors the user is kindly requested to contact the authors, sending a
copy of both input and output by E-mail to the Torino group (crystal@unito.it) or to the
Daresbury team (crystal@dl.ac.uk).
10
Chapter 1
Wave function calculation
Basic input route
1.1
Geometry and symmetry information
The first record of the geometry definition must contain one of the keywords:
CRYSTAL
3D system
SLAB
2D system
POLYMER
1D system
MOLECULE
0D system
EXTERNAL
geometry from external file
DLVINPUT
geometry from DLV [18] Graphical User Interface.
Three input schemes are used. The first is for crystalline systems, and is specified by the
keyword CRYSTAL. The second is for slabs, polymers and molecules as specified by the key-
words SLAB, POLYMER or MOLECULE respectively. In the third scheme, with keyword
EXTERNAL or DLVINPUT, the unit cell, atomic positions and symmetry operators may
be provided directly (see Appendix F, page 181). Such an input file can be prepared by the
keyword EXTPRT (input block 1, page 31; properties). Sample input decks for a number
of structures are provided in section 4.1, page 115.
11
Geometry input for crystalline compounds
rec
variable
value
meaning
·
IFLAG
convention for space group identification (Appendix A.1, page 155):
0
space group sequential number(1-230)
1
Hermann-Mauguin alphanumeric code
IFHR
type of cell for rhombohedral groups (meaningless for non-
rhombohedral crystals):
0
hexagonal cell
1
rhombohedral cell
IFSO
setting for the origin of the crystal reference frame:
0
origin derived from the symbol of the space group: where there
are two settings, the second setting of the International Tables is
chosen.
1
standard shift of the origin: when two settings are allowed, the first
setting is chosen
>1
non-standard shift of the origin given as input
·
space group identification code (following IFLAG value):
IGR
space group sequence number (IFLAG=0)
or
A
AGR
space group alphanumeric symbol (IFLAG=1)
if IFSO > 1 insert
II
·
IX,IY,IZ
non-standard shift of the origin coordinates (x,y,z) in fractions of
the crystallographic cell lattice vectors times 24 (to obtain integer
values).
·
a,[b],[c],
minimal set of crystallographic cell parameters:
[],[]
translation vector[s] length [°
Angstrom],
[]
crystallographic angle[s] (degrees)
·
NATR
number of atoms in the asymmetric unit.
insert NATR records
II
·
NAT
conventional atomic number 1
X,Y,Z
atom coordinates in fractional units of crystallographic lattice vec-
tors
optional keywords terminated by END/ENDGEOM or STOP
II
Geometry input for molecules, polymers and slabs
When the geometrical structure of 2D, 1D and 0D systems has to be defined, attention should
be paid in the input of the atom coordinates, that are expressed in different units, fractionary
(direction with translational symmetry) or °
Angstrom (non periodic direction).
translational
unit of measure of coordinates
symmetry
X
Y
Z
3D
fraction
fraction
fraction
2D
fraction
fraction
°
Angstrom
1D
fraction
°
Angstrom °
Angstrom
0D
°
Angstrom °
Angstrom °
Angstrom
12
rec variable
meaning
· IGR
point, rod or layer group of the system:
0D - molecules (Appendix A.4, page 162)
1D - polymers (Appendix A.3, page 159)
2D - slabs (Appendix A.2, page 158)
if polymer or slab, insert
II
· a,[b],
minimal set of lattice vector(s)- length in °
Angstrom
(b for rectangular lattices only)
[]
AB angle (degrees) - triclinic lattices only
· NATR
number of non-equivalent atoms in the asymmetric unit
insert NATR records
II
· NAT
conventional atomic number
X,Y,Z
atoms coordinates. Unit of measure:
0D - molecules: x,y,z in °
Angstrom
1D - polymers : y,z in °
Angstrom, x in fractional units of crystallographic
cell translation vector
2D - slabs : z in °
Angstrom, x, y in fractional units of crystallographic cell
translation vectors
optional keywords terminated by END or STOP
II
Geometry input from external geometry editor
The keywords EXTERNAL and DLVINPUT select the third input scheme. They work
for molecules, polymers, slabs and crystals. The input data are read from Fortran unit 34.
The unit cell, atomic positions and symmetry operators are provided directly according to the
format described in Appendix F, page 181. Coordinates in °
Angstrom. Such an input file is
written when OPTBERNY route for geometry optimization is chosen, and can be prepared by
the keyword EXTPRT (input block 1, page 31; properties).
The geometry so defined can be modified by inserting any geometry editing keyword (page 22)
after EXTERNAL.
Comments on geometry input
1. All coordinates in °
Angstrom. In geometry editing the unit of measure of coordinates may
be modified by entering the keywords FRACTION (page 31) or BOHR (page 27).
2. The geometry of a system is defined by the crystal structure ([19], Chapter 1 of ref. [20]).
Reference is made to the International Tables for Crystallography [21] for all definitions.
The crystal structure is determined by the space group, by the shape and size of the unit
cell and by the relative positions of the atoms in the asymmetric unit.
3. The lattice parameters represent the length of the edges of the cell (a,b,c) and the angles
between the edges ( = b c; = a c; = a b). They determine the cell volume and
shape.
4. Minimal set of lattice parameters to be defined in input:
cubic
a
hexagonal
a,c
rhombohedral
hexagonal cell
a,c
rhombohedral cell
a,
tetragonal
a,c
orthorhombic
a,b,c
monoclinic
a,b,c, (b unique)
a,b,c, (c unique)
a,b,c, (a unique - non standard)
triclinic
a,b,c, , ,
13
5. The asymmetric unit is the largest subset of atoms contained in the unit-cell, where
no atom pair related by a symmetry operator can be found. Usually several equivalent
subsets of this kind may be chosen so that the asymmetric unit needs not be unique.
The asymmetric unit of a space group is a part of space from which, by application of
all symmetry operations of the space group, the whole of space is filled exactly.
6. The crystallographic, or conventional cell, is used as the standard option in input. It
may be non-primitive, which means it may not coincide with the cell of minimum volume
(primitive cell) which contains just one lattice point. The matrices which transform the
conventional (as given in input) to the primitive cell (used by CRYSTAL) are given in
Appendix A.5, page 163, and are taken from table 5.1 of the International Tables [21].
Examples. A cell belonging to the face-centred cubic Bravais lattice has a volume four
times larger than that of the corresponding primitive cell, and contains four lattice points
(see page 38, keyword SUPERCEL). A unit cell belonging to the hexagonal Bravais
lattice has a volume three times larger than that of the rhombohedral primitive cell (R
Bravais lattice), and contains three lattice points.
7. The use of the International Tables to identify the symmetry groups requires some prac-
tice. The examples given below may serve as a guide. The printout of geometry informa-
tion (equivalent atoms, fractional and Cartesian atomic coordinates etc.) allows a check
on the correctness of the group selected. To obtain a complete neighborhood analysis
for all the non-equivalent atoms, a complete input deck must be read in (blocks 1-4),
and the keyword TESTPDIM inserted in block 3, to stop execution after the symmetry
analysis.
8. Different settings of the origin may correspond to a different number of symmetry oper-
ators with translational components.
Example: bulk silicon - Space group 227 - 1 irreducible atom per cell.
setting of the origin
Si coordinates
symmops with
translational component
2nd (default)
1/8
1/8
1/8
36
1st
0.
0.
0.
24
NB With 2nd setting, the position 0., 0., 0. has multiplicity 4.
The choice is important when generating a supercell, as the first step is the removal of the
symmops with translational component. The keyword ORIGIN (input block 1, page
34) translates the origin in order to minimize the number of symmops with translational
component.
9. When coordinates are obtained from experimental data or from geometry optimization
with semi-classical methods, atoms in special positions, or related by symmetry are not
correctly identified, as the number of significative digits is lower that the one used by
the program crystal to recognize equivalence or special positions.
In that case the
coordinates must be edited by hand (see FAQ at www.crystal.unito.it).
10. The symbol of the space group for crystals (IFLAG=1) is given precisely as it appears
in the International Tables, with the first letter in column one and a blank separating
operators referring to different symmetry directions. The symbols to be used for the
groups 221-230 correspond to the convention adopted in editions of the International
Tables prior to 1983: the 3 axis is used instead of 3. See Appendix A.1, page 155.
Examples:
Group number
input symbol
137 (tetragonal)
P 42/N M C
10 (monoclinic)
P 1 2/M 1
(unique axis b, standard setting)
P 1 1 2/M
(unique axis c)
14
P 2/M 1 1
(unique axis a)
25 (orthorhombic)
P M M 2
(standard setting)
P 2 M M
P M 2 M
11. In the monoclinic and orthorhombic cases, if the group is identified by its number (3-74),
the conventional setting for the unique axis is adopted. The explicit symbol must be
used in order to define an alternative setting.
12. For the centred lattices (F, I, C, A, B and R) the input cell parameters refer to the
centred conventional cell; the fractional coordinates of the input list of atoms are in a
vector basis relative to the centred conventional cell.
13. It is sufficient to supply the coordinates of only one of a group of atoms equivalent under
centring translations (eg: for space group Fm3m only the parameters of the face-centred
cubic cell, and the coordinates of one of the four atoms at (0,0,0), (0, 1 , 1 ), ( 1 ,0, 1 ) and
2 2
2
2
( 1 , 1 ,0) are required).
2 2
The coordinates of only one atom among the set of atoms linked by centring translations
are printed. The vector basis is relative to the centred conventional cell. However when
Cartesian components of the direct lattice vectors are printed, they are those of the
primitive cell.
14. The conventional atomic number NAT is used to associate a given basis set with an
atom (see Basis Set input, Section 1.2, page 16). The real atomic number is given by the
remainder of the division of the conventional atomic number by 100 (Example: NAT=237,
Z=37; NAT=128, Z=28). Atoms with the same atomic number, but in non-equivalent
positions, can be associated with different basis sets, by using different conventional
atomic numbers: e.g. 6, 106 (all electron basis set for carbon atom); 206, 306 (core
pseudo-potential for carbon atom, Section 2.2, page 51).
If the remainder of the division is 0, a "ghost" atom is identified, to which no nuclear
charge corresponds (it may have electronic charge). This option may be used for enriching
the basis set by adding bond basis function [22], or to allow build up of charge density on
a vacancy. A given atom may be transformed into a ghost after the basis set definition
(input block 2, keyword GHOSTS, page 50).
15. The keyword SLAB (Geometry editing input, page 37) allows the creation of a slab
(2D) of given thickness from the 3D perfect lattice. See for comparison test4-test24;
test5-test25; test6-test26; test7- test27.
16. For slabs (2D), when two settings of the origin are indicated in the International Tables
for Crystallography, setting number 2 is chosen. The setting can not be modified.
17. Conventional orientation of slabs and polymers: Polymers are oriented along the x axis.
Slabs are parallel to the xy plane.
18. The keywords MOLECULE (for molecular crystals only; page 33) and CLUSTER
(for any n-D structure; page 28) allow the creation of a non-periodic system (molecule(s)
or cluster) from a periodic one.
15
1.2
Basis set
rec variable value
meaning
· NAT
n
conventional atomic number
<200
all-electron basis set
>200
valence electron basis set. ECP (Effective Core Pseudopotential)
must be defined (page 51)
=99
end of basis set input section
NSHL
n
number of shells
0
end of basis set input (when NAT=99)
if NAT > 200 insert ECP input (page 51)
II
NSHL sets of records - for each shell
· ITYB
type of basis set to be used for the specified shell:
0
general BS, given as input
1
Pople standard STO-nG (Z=1-54)
2
Pople standard 3(6)-21G (Z=1-54(18)) Standard polarization func-
tions are included.
LAT
shell type:
0
1 s AO (S shell)
1
1 s + 3 p AOs (SP shell)
2
3 p AOs (P shell)
3
5 d AOs (D shell)
NG
Number of primitive Gaussian Type Functions (GTF) in the con-
traction for the basis functions (AO) in the shell
1NG10 for ITYB=0 and LAT 2
1NG6
for ITYB=0 and LAT = 3
2NG6
for ITYB=1
6
6-21G core shell
3
3-21G core shell
2
n-21G inner valence shell
1
n-21G outer valence shell
CHE
formal electron charge attributed to the shell
SCAL
scale factor (if ITYB=1 and SCAL=0., the standard Pople scale
factor is used for a STO-nG basis set.
if ITYB=0 (general basis set insert NG records
II
· EXP
exponent of the normalized primitive GTF
COE1
contraction coefficient of the normalized primitive GTF:
LAT=0,1 s function coefficient
LAT=2 p function coefficient
LAT=3 d function coefficient
COE2
LAT=1 p function coefficient
optional keywords terminated by END/ENDB or STOP
II
The choice of basis set is the most critical step in performing ab initio calculations of periodic
systems, with Hartree-Fock or Kohn-Sham Hamiltonians. Optimization criteria are discussed in
Chapter 3.2. When an effective core pseudo-potential is used, the basis set must be optimized
with reference to that potential (Section 2.2, page 51).
1. A basis set (BS) must be given for each atom with different conventional atomic number
defined in the crystal structure input. If atoms are removed (geometry input, keyword
ATOMREMO, page 25), the corresponding basis set input can remain in the input
stream.
2. The basis set for each atom has NSHL shells whose constituent AO basis functions
are built from a linear combination ('contraction') of individually normalized primitive
Gaussian-type functions (GTF) (Chapter 6, page 139).
3. A conventional atomic number NAT links the basis set with the atoms defined in the
16
crystal structure. The atomic number Z is given by the remainder of the division of the
conventional atomic number by 100 (Example: NAT=108, Z=8, all electron; NAT=228,
Z=28, ECP). See point 5 below.
4. A conventional atomic number 0 defines ghost atoms, that is points in space with an
associated basis set, but lacking a nuclear charge (vacancy). See test 28.
5. Atoms with equal conventional atomic number are associated with the same basis set.
NAT < 200: all electron basis set. A maximum of two different basis sets
may be given for the same chemical species in different positions:
NAT=Z, NAT=Z+100.
NAT > 200: valence electron basis set.
A maximum of two different BS
may be given for the same chemical species in positions not
symmetry-related: NAT=Z+200, NAT=Z+300. A core pseudo-
potential must be defined. See Section 2.2, page 51, for informa-
tion on core pseudo-potentials.
Suppose we have four non-equivalent carbon atoms in the unit cell. Conventional atomic
numbers 6 106 206 306 mean that carbon atoms (real atomic number 6) unrelated by
symmetry are to be associated with different basis sets: the first two (6, 106) all-electron,
the second two (206, 306) valence only, with pseudo-potential.
6. The basis set input ends with the card:
99
0
conventional atomic number 99, 0 shell.
The optional keywords may follow.
In summary:
1. CRYSTAL can use the following all electrons basis sets:
a)
general basis sets, including s,p,d functions (given in input);
b)
standard Pople basis sets [23] (internally stored)
STOnG,
Z=1 to 54
6-21G,
Z=1 to 18
3-21G,
Z=1 to 54
The standard basis sets b) are stored as internal data in the CRYSTAL code. They are
all electron basis sets, and can not be combined with ECP.
2. Warning The standard scale factor is used for STO-nG basis set when the input datum
SCAL is 0.0 in basis set input. All the atoms of the same row are attributed the same
Pople STO-nG basis set when the input scale factor SCAL is 1.
3. Standard polarization functions can be added to 6(3)-21G basis sets of atoms up to Z=18,
by inserting a record describing the polarization shell (ITYB=2, LAT=2, p functions on
hydrogen, or LAT=3, d functions on 2-nd row atoms; see test 12).
H
Polarization functions exponents
He
1.1
1.1
__________
______________________________
Li
Be
B
C
N
O
F
Ne
0.8
0.8
0.8
0.8
0.8
0.8
0.8
--
___________
______________________________
Na
Mg
Al
Si
P
S
Cl
Ar
0.175
0.175
0.325 0.45 0.55 0.65 0.75 0.85
_____________________________________________________________________
The formal electron charge attributed to a polarization function must be zero.
4. The shell types available are :
17
shell shell
n.
order of internal storage
code type
AO
0
S
1
s
1
SP
4
s, x, y, z
2
P
3
x, y, z
3
D
5
2z2 - x2 - y2, xz, yz, x2 - y2, xy
4
F
7
(2z2 - 3x2 - 3y2)z, (4z2 - x2 - y2)x, (4z2 - x2 - y2)y,
(x2 - y2)z, xyz, (x2 - 3y2)x, (3x2 - y2)y
The order of internal storage of the AO basis functions is an information necessary to
read certain quantities calculated by the program properties. See Chapter 3: Mul-
liken population analysis (PPAN, page 78), electrostatic multipoles (POLI, page 105),
projected density of states (DOSS,page 91) and to provide an input for some options
(EIGSHIFT, input block 4, page 74).
5. Spherical harmonics d-shells consisting of 5 AOs are used.
6. The formal shell charges CHE, the number of electrons attributed to each shell, are
assigned to the AO following the rules:
shell
shell
max
rule to assign the shell charges
code
type
CHE
0
S
2.
CHE to S functions
1
SP
8.
if CHE>2, 2 to S and (CHE-2) to P functions,
if CHE2, CHE to S function
2
P
6.
CHE to P functions
3
D
10.
CHE to D functions
4
F
14.
CHE to F functions
5
G
18.
CHE to G functions
7. A maximum of one open shell for each of the s, p and or d atomic symmetries is allowed
in the electronic configuration defined in the input. The atomic energy expression is not
correct for all possible double open shell couplings of the form pmdn. Either m must
equal 3 or n must equal 5 for a correct energy expression in such cases. A warning
will be printed if this is the case. However, the resultant wave function (which is a
superposition of atomic densities) will usually provide a reasonable starting point for the
periodic density matrix.
8. When extended basis sets are used, all the functions corresponding to symmetries (an-
gular quantum numbers) occupied in the isolated atom are added to the atomic basis
set for atomic wave function calculations, even if the formal charge attributed to that
shell is zero. Polarization functions are not included in the atomic basis set; their input
occupation number should be zero.
9. The formal shell charges are used only to define the electronic configuration of the atoms
to compute the atomic wave function. The initial density matrix in the scf step may be
built as a superposition of atomic (or ionic) density matrices (default option, GUESS-
PAT 2.4). When a different guess is required (GUESSF or GUESSP), the shell charges
are not used, but checked for electron neutrality when the basis set is entered.
10. Each atom in the cell may have an ionic configuration, when the sum of formal shell
charges (CHE) is different from the nuclear charge. When the number of electrons in
the cell, that is the sum of the shell charges CHE of all the atoms, is different from the
sum of nuclear charges, the reference cell is non-neutral. This is not allowed for periodic
systems, and in that case the program stops. In order to remove this constraint, it is
necessary to introduce a uniform charged background of opposite sign to neutralize the
system [24]. This is obtained by entering the keyword CHARGED (page 48) after the
standard basis set input.
18
11. It may be useful to allow atoms with the same basis set to have different electronic
configurations (e.g, for an oxygen vacancy in MgO one could use the same basis set for
all the oxygens, but begin with different electronic configuration for those around the
vacancy). The formal shell charges attributed in the basis set input may be modified for
selected atoms by inserting the keyword CHEMOD (input block 2, page 48).
12. The energies given by an atomic wave function calculation with a crystalline basis set
should not be used as a reference to calculate the formation energies of crystals. The
external shells should first be re-optimized in the isolated atom by adding a low-exponent
Gaussian function, in order to provide and adequate description of the tails of the isolated
atom charge density [25] (keyword ATOMHF, input block 3, page 57).
Optimized basis sets for periodic systems used in published papers are available on WWW:
http://www.crystal.unito.it
1.3
General information, computational parameters,
hamiltonian
No input is required if the default values are used. Note however that END/ENDM or
STOP, to close the section, are always needed. If no Hamiltonian is specified, RHF (Restricted
Hartree-Fock Hamiltonian) is assumed.
1.4
SCF input
rec variable value
meaning
if the system is periodic insert
II
· IS
Shrinking factor in reciprocal space (Section 6.7, page 145)
IDUM
not used.
ISP
Shrinking factor for a denser k point net (Gilat net) in the For peri-
evaluation of the Fermi energy and density matrix.
if IS = 0 insert
II
· IS1,IS2,IS3
Shrinking factors along B1,B2,B3 (reciprocal lattice vectors);
to be used when the unit cell is highly anisotropic
optional keywords terminated by END or STOP
II
odic systems, 1D, 2D, 3D, the mandatory input information is the shrinking factor, IS, to
generate a commensurate grid of k points in reciprocal space, according to Pack-Monkhorst
method. The Hamiltonian matrix computed in direct space, Hg, is Fourier transformed for
each k value, and diagonalized, to obtain eigenvectors and eigenvalues:
Hk =
Hgeigk
g
HkAk = SkAkEk
A second shrinking factor, ISP, defines the sampling of k points, i "Gilat net", used for the cal-
culation of the density matrix and the determination of Fermi energy in the case of conductors,
when bands are not fully occupied.
In 3D crystals, the sampling points belong to a lattice i (called the Pack-Monkhorst net), with
basis vectors:
b1/is1, b2/is2, b3/is3
is1=is2=is3=IS, unless otherwise stated
where b1, b2, b3 are the reciprocal lattice vectors, and is1, is2, is3 are integers "shrinking
factors".
19
In 2D crystals, IS3 is set equal to 1; in 1D crystals both IS2 and IS3 are set equal to 1. Only
points ki of the Pack-Monkhorst net belonging to the irreducible part of the Brillouin Zone
(IBZ) are considered, with associated a geometrical weight, wi. The choice of the reciprocal
space integration parameters to compute the Fermi energy is a delicate step for metals. See
Section 6.7, page 145. Two parameters control the accuracy of reciprocal space integration for
Fermi energy calculation and density matrix reconstruction:
IS shrinking factor of reciprocal lattice vectors.
The value of IS determines the number
of k points at which the Fock/KS matrix is diagonalized. Multiples of 2 or 3 should
be used, according to the point symmetry of the system (order of principal axes). The
k-points net is automatically made anisotropic for 1D and 2D systems.
The figure presents the reciprocal lattice cell of 2D graphite (rhombus), the first
Brillouin zone (hexagon), the irreducible part of Brillouin zone (in gray), and the
coordinates of the ki points according to a Monkhorst-Pack sampling, with shrinking
factor 3 and 6.
ISP shrinking factor of reciprocal lattice vectors in the Gilat net (see [7], Chapter II.6). ISP
is used in the calculation of the Fermi energy and density matrix. Its value can be equal
to IS for insulating systems and equal to 2*IS for conducting systems.
Note. The value used in the calculation is ISP=IS*NINT(MAX(ISP,IS)/IS)
1. When an anisotropic net is user defined (IS=0), the ISP input value is taken as ISP1
(shrinking factor of Gilat net along first reciprocal lattice) and ISP2 and ISP3 are set to:
ISP2=(ISP*IS2)/IS1,
ISP3=(ISP*IS3)/IS1.
2. User defined anisotropic net is not compatible with SABF (Symmetry Adapted Bloch
Functions). See NOSYMADA, page 78.
Some tools for accelerating convergence are given through the keywords LEVSHIFT (page
77 and tests 29, 30, 31, 32, 38), FMIXING (page 75), SMEAR (page 79), BROYDEN
(page 73) and ANDERSON (page 72).
At each SCF cycle the total atomic charges, following a Mulliken population analysis scheme,
and the total energy are printed.
The defaulty value of the parameters to control the exit from the SCF cycle (E < 10-6
hartree, maximum number of SCF cycles: 50) may be modified entering the keywords:
TOLSCF (tolerance on change in eigenvalues and total energy),
TOLENE (tolerance on change in total energy),
TOLDEP (tolerance on SQM in density matrix elements),
MAXCYCLE (maximum number of cycles).
20
Spin-polarized system
By default the orbital occupancies are controlled according to the 'Aufbau' principle.
To obtain a spin polarized solution an open shell Hamiltonian must be defined (block3, UHF
or DFT/SPIN). A spin-polarized solution may then be computed after definition of the (-)
electron occupancy. This can be performed by the keywords SPINLOCK and BETALOCK.
21
Chapter 2
Wave function calculation -
Advanced input route
2.1
Geometry editing
The following keywords allow editing of the crystal structure, printing of extended information,
generation of input data for visualization programs. Processing of the input block 1 only
(geometry input) is allowed by the keyword TESTGEOM.
Each keyword operates on the geometry active when the keyword is entered. For instance, when
a 2D structure is generated from a 3D one through the keyword SLABCUT, all subsequent
geometry editing operates on the 2D structure. When a dimer is extracted from a molecular
crystal through the keyword MOLECULE, all subsequent editing refers to a system without
translational symmetry.
The keywords can be entered in any order: particular attention should be paid to the action
of the keywords KEEPSYMM and BREAKSYM, that allow maintaining or breaking the
symmetry while editing the structure. These keywords behave as a switch, and require no
further data. Under control of the BREAKSYM keyword (the default), subsequent mod-
ifications of the geometry are allowed to alter (reduce: the number of symmetry operators
cannot be increased) the point-group symmetry. The new group is a subgroup of the original
group and is automatically obtained by CRYSTAL. However if a KEEPSYMM keyword
is presented, the program will endeavor to maintain the number of symmetry operators, by
requiring that atoms which are symmetry related remain so after a geometry editing (key-
words:ATOMSUBS, ATOMINSE, ATOMDISP, ATOMREMO).
The space group of the system may be modified after editing.
For 3D systems,
the file FINDSYM.DAT is written.
This file is input to the program findsym
(http://physics.byu.edu/ stokesh/isotropy.html), that finds the space-group symmetry of a
crystal given the coordinates of the atoms.
Geometry keywords
Symmetry information
ATOMSYMM
printing of point symmetry at the atomic positions
27
MAKESAED
printing of symmetry allowed elastic distortions (SAED)
32
PRSYMDIR
printing of displacement directions allowed by symmetry.
35
SYMMDIR
printing of symmetry allowed geom opt directions
40
SYMMOPS
printing of point symmetry operators
40
TENSOR
tensor of physical properties
40
I
Symmetry information and control
22
BREAKSYM
allow symmetry reduction following geometry modifications
28
KEEPSYMM
maintain symmetry following geometry modifications
32
MODISYMM
removal of selected symmetry operators
32
I
PURIFY
cleans atomic positions so that they are fully consistent with the 35
group
SYMMREMO
removal of all symmetry operators
40
TRASREMO
removal of symmetry operators with translational components
40
Modifications without reduction of symmetry
ATOMORDE
reordering of atoms in molecular crystals
25
NOSHIFT
no shift of the origin to minimize the number of symmops with 34
translational components before generating supercell
ORIGIN
shift of the origin to minimize the number of symmetry operators 34
with translational components
PRIMITIV
crystallographic cell forced to be the primitive cell
35
REDEFINE
definition of a new cell, with xy
to a given plane
36
I
Atoms and cell manipulation (possible symmetry reduction (BREAKSYMM)
ATOMDISP
displacement of atoms
25
I
ATOMINSE
addition of atoms
25
I
ATOMREMO
removal of atoms
25
I
ATOMROT
rotation of groups of atoms
26
I
ATOMSUBS
substitution of atoms
27
I
ELASTIC
distortion of the lattice
30
I
USESAED
given symmetry allowed elastic distortions, reads
41
I
SUPERCEL
generation of supercell - input refers to primitive cell
38
I
SUPERCON
generation of supercell - input refers to conventional cell
38
I
From crystals to slabs
SLABCUT
generation of a slab parallel to a given plane (3D2D)
37
I
From periodic structure to clusters
CLUSTER
cutting of a cluster from a periodic structure (3D0D)
28
I
HYDROSUB
border atoms substituted with hydrogens (0D0D)
31
I
Molecular crystals
MOLECULE
extraction of a set of molecules from a molecular crystal 33
I
(3D0D)
MOLEXP
variation of lattice parameters at constant symmetry and molec- 33
I
ular geometry (3D3D)
MOLSPLIT
periodic structure of non interacting molecules (3D3D)
34
RAYCOV
modification of atomic covalent radii
35
I
BSSE correction
MOLEBSSE
counterpoise method for molecules (molecular crystals only) 32
I
(3D0D)
ATOMBSSE
counterpoise method for atoms (3D0D)
25
I
Auxiliary and control keywords
23
ANGSTROM
sets inputs unit to °
Angstrom
24
BOHR
sets input units to bohr
27
BOHRANGS
input bohr to °
A conversion factor (0.5291772083 default value) 27
I
BOHRCR98
bohr to °
A conversion factor is set to 0.529177 (CRYSTAL98
value)
END/ENDG
terminate processing of geometry
FRACTION
sets input unit to fractional
31
NEIGHBOR
number of neighbours in geometry analysis
34
I
PARAMPRT
printing of parameters controlling dimensions of static allocation 35
arrays
PRINTOUT
setting of printing options by keywords
35
I
SETINF
setting of inf array options
36
I
SETPRINT
setting of printing options
37
I
STOP
execution stops immediately
38
TESTGEOM
stop after checking the geometry
40
Output of data on external units
COORPRT
coordinates of all the atoms in the cell
29
EXTPRT
generation of file as CRYSTAL input
31
MOLDRAW
generation of file for the program MOLDRAW
32
STRUCPRT
cell parameters and coordinates of all the atoms in the cell
38
External electric field - modified Hamiltonian
FIELD
external field applied (2D-3D systems only)
41
I
Geometry optimization
OPTCOORD
Atom coordinates optimization
43
I
Initial Hessian
HESGUESS
initial guess for the Hessian
I
Convergence criteria
TOLDEG
RMS of the gradient [0.0003]
I
TOLDEX
RMS of the displacement [0.0012]
I
TOLDEE
energy difference between two steps [10-7]
I
MAXOPTC
max number of optimization steps
I
Optimization control
ATOMFREE
partial geometry optimization
I
RESTART
data from previous run
FINALRUN
Wf single point with optimized geometry
I
Gradient calculation control
NUMGRAD
numerical first deivatives
Printing options
PRINTFORCES atomic gradients
PRINTHESS
Hessian
PRINTOPT
optimization procedure
PRINT
verbose printing
ANGSTROM - unit of measure
The unit of length in geometry editing is set to °
Angstrom, (default value).
24
ATOMBSSE - counterpoise for closed shell atoms and ions
rec variable meaning
· IAT
label of the atom in the reference cell
NSTAR maximum number of stars of neighbors included in the calculation.
RMAX
maximum distance explored searching the neighbors of the atom.
A cluster is defined including the selected atom and the basis functions belonging to the NSTAR
sets of neighbours, when their distance R from the central atom is smaller than RMAX. The
atomic wave function is not computed by the atomic package, but by the standard CRYSTAL
route for 0D, 1 atom system. UHF and SPINLOCK must be used to define a reasonable
orbital occupancy. It is suggested to compute the atomic wave function using a program
properly handling the electronic configuration of open shell atoms.
Warning. The system is 0D. No reciprocal lattice information is required in the scf input
(Section 1.4, page 19).
ATOMDISP
rec variable
meaning
· NDISP
number of atoms to be displaced
insert NDISP records
II
· LB
label of the atom to be moved
DX,DY,DZ
increments of the coordinates in the primitive cell [°
A].
Selected atoms are displaced in the primitive cell. The point symmetry of the system may be
altered (default value BREAKSYM, page 28). Increments are in °
Angstrom, unless otherwise
requested (keyword BOHR, FRACTION, page 24). See tests 17, 20, 37.
ATOMINSE
rec variable
meaning
· NINS
number of atoms to be added
insert NINS records
II
· NA
conventional atomic number
X,Y,Z
coordinates [°
A] of the inserted atom. Coordinates refer to the primitive cell.
New atoms are added to the primitive cell. Coordinates are in °
Angstrom, unless otherwise
requested (keyword BOHR, FRACTION, page 24). Remember that the original symmetry
of the system is maintained, applying the symmetry operators to the added atoms if the
keyword KEEPSYMM (page 28) was previously entered. The default is BREAKSYM
(page 28). Attention should be paid to the neutrality of the cell (see CHARGED, page 48).
See tests 16, 35, 36.
ATOMORDE
After processing the standard geometry input, the symmetry equivalent atoms in the reference
cell are grouped. They may be reordered, following a chemical bond criterion. This simplifies
the interpretation of the output when the results of bulk molecular crystals are compared with
those of the isolated molecule. See option MOLECULE (page 33) and MOLSPLIT (page
34). No input data are required.
ATOMREMO
rec variable
meaning
· NL
number of atoms to remove
· LB(L),L=1,NL label of the atoms to remove
25
Selected atoms, and related basis set, are removed from the primitive cell (see test 16). A
vacancy is created in the lattice. The symmetry can be maintained (KEEPSYMM), by
removing all the atoms symmetry-related to the selected one, or reduced (BREAKSYM).
Attention should be paid to the neutrality of the cell (see CHARGED, page 48).
NB. The keyword GHOSTS (basis set input, page 50) allows removal of selected atoms,
leaving the related basis set.
ATOMROT
rec variable value
meaning
· NA
0
all the atoms of the cell are rotated and/or translated
>0
only NA selected atoms are rotated and/or translated.
<0
the atom with label |NA| belongs to the molecule to be rotated. The
program selects all the atoms of the molecule on the base of the sum of
their atomic radii (Table on page 35).
if NA > 0, insert NA data
II
·
LB(I),I=1,NA label of the atoms to be rotated and/or translated.
· ITR
>0
translation performed. The selected NA atoms are translated by -r,
where r is the position of the ITR-th atom. ITR is at the origin after
the translation.
0
a general translation is performed. See below.
=999
no translation.
IRO
> 0
a rotation around a given axis is performed. See below.
< 0
a general rotation is performed. See below.
=999
no rotation.
if ITR<0 insert
II
· X,Y,Z
Cartesian components of the translation vector [°
A]
if ITR=0 insert
II
· N1,N2
label of the atoms defining the axis.
DR
translation along the axis defined by the atoms N1 and N2, in the di-
rection N1 N2 [°
A].
if IRO<0 insert
II
· A,B,G
Euler rotation angles (degree).
IPAR
defines the origin of the Cartesian system for the rotation
0
the origin is the barycentre of the NAT atoms
>0
the origin is the atom of label IPAR
if IRO>0 insert
II
· N1,N2
label of the atoms that define the axis for the rotation
ALPHA = 0.
rotation angle around the N1N2 axis (degrees)
0.
the selected atoms are rotated anti-clockwise in order to orientate the
N1N2 axis parallel to the z axis.
This option allows to rotate and/or translate the specified atoms. When the rotation of a
molecule is required (NA < 0), the value of the atomic radii must be checked, in order to
obtain a correct definition of the molecule. It is useful to study the conformation of a molecule
in a zeolite cavity, or the interaction of a molecule (methane) with a surface (MgO).
The translation of the selected group of atoms can be defined in three different ways:
1. Cartesian components of the translation vector (ITR < 0);
2. modulus of the translation vector along an axis defined by two atoms (ITR = 0);
3. sequence number of the atom to be translated to the origin. All the selected atoms are
subjected to the same translation (ITR > 0).
The rotation can be performed in three different ways:
26
1. by defining the Euler rotation angles , , and the origin of the rotating system (IRO
< 0). The axes of the rotating system are parallel to the axes of the Cartesian reference
system. (The rotation is given by: RzRxRz, where R are the rotation matrices).
2. by defining the rotation angle around an axis defined by two atoms A and B. The
origin is at A, the positive direction AB.
3. by defining a z' axis (identified by two atoms A and B). The selected atoms are rotated,
in such a way that the AB z' axis becomes parallel to the z Cartesian axis. The origin
is at A and the positive rotation anti clockwise (IRO>0, =0).
The selected atoms are rotated according to the defined rules, the cell orientation and the
cartesian reference frame are not modified. The symmetry of the system is checked after the
rotation, as the new geometry may have a different symmetry.
See tests 15, rotation of the N H3 molecule in a zeolite cavity, and 16, rotation of the H2O
molecule in the zeolite cavity.
ATOMSUBS
rec variable
meaning
· NSOST
number of atoms to be substituted
insert NSOST records
II
· LB
label of the atom to substitute
NA(LB)
conventional atomic number of the new atom
Selected atoms are substituted in the primitive cell (see test 17, 34, 37). The symmetry can be
maintained (KEEPSYMM), by substituting all the atoms symmetry-related to the selected
one, or reduced (BREAKSYM). Attention should be paid to the neutrality of the cell: a non-
neutral cell will cause an error message, unless allowed by entering the keyword CHARGED,
page 48.
ATOMSYMM
The point group associated with each atomic position and the set of symmetry related atoms
are printed. No input data are required. This option is useful to find the internal coordinates
to be relaxed when the unit cell is deformed (see ELASTIC, page 30).
BOHR
The keyword BOHR sets the unit of distance to bohr. When the unit of measure is modified,
the new convention is active for all subsequent geometry editing.
The conversion factor °
Angstrom/bohr is 0.5291772083 (CODATA 1998). This value can be
modified by entering the keyword BOHRANGS and the desired value in the record following.
The keyword BOHRCR98 sets the conversion factor to 0.529177, as in the program CRYS-
TAL98.
CRYSTAL88 default value was 0.529167).
BOHRANGS
rec variable
meaning
· BOHR
conversion factor °
Angstrom/bohr
The conversion factor °
Angstrom/bohr can be user-defined.
In CRYSTAL88 the default value was 0.529167.
In CRYSTAL98 the default value was 0.529177.
27
BOHRCR98
The conversion factor °
Angstrom/bohr is set to 0.529177, as in CRYSTAL98. No input data
required.
BREAKSYM
Under control of the BREAKSYM keyword (the default), subsequent modifications of the
geometry are allowed to alter (reduce: the number of symmetry operators cannot be increased)
the point-group symmetry. The new group is a subgroup of the original group and is automat-
ically obtained by CRYSTAL.
The symmetry may be broken by attributing different spin (ATOMSPI, block4, page 73) to
atoms symmetry related by geometry.
Example: When a CO molecule is vertically adsorbed on a (001) 3-layer MgO slab, (D4h
symmetry), the symmetry is reduced to C4v, if the BREAKSYM keyword is active. The
symmetry operators related to the h plane are removed. However, if KEEPSYMM is
active, then additional atoms will be added to the underside of the slab so as to maintain the
h plane (see page 25, keyword ATOMINSE).
CLUSTER - a cluster (0D) from a periodic system
The CLUSTER option allows one to cut a finite molecular cluster of atoms from a periodic
lattice. The size of the cluster (which is centred on a specified 'seed point' A) can be controlled
either by including all atoms within a sphere of a given radius centred on A, or by specifying
a maximum number of symmetry-related stars of atoms to be included.
The cluster includes the atoms B (belonging to different cells of the direct lattice) satisfying
the following criteria:
1. those which belong to one of the first N (input data) stars of neighbours of the seed point
of the cluster.
and
2. those at a distance RAB from the seed point which is smaller then RMAX (input datum).
The resulting cluster may not reproduce exactly the desired arrangement of atoms, particularly
in crystals with complex structures such as zeolites, and so it is possible to specify border
modifications to be made after definition of the core cluster.
Specification of the core cluster:
rec variable value
meaning
· X, Y, Z
coordinates of the centre of the cluster [°
A] (the seed point)
NST
maximum number of stars of neighbours explored in defining the core
cluster
RMAX
radius of a sphere centred at X,Y,Z containing the atoms of the core
cluster
· NNA
= 0
print nearest neighbour analysis of cluster atoms (according to a radius
criterion)
NCN
0
testing of coordination number during hydrogen saturation carried out
only for Si (coordination number 4), Al (4) and O(2)
N
N user-defined coordination numbers are to be defined
if NNA = 0 insert 1 record
II
· RNNA
radius of sphere in which to search for neighbours of a given atom in
order to print the nearest neighbour analysis
if NCN = 0 insert NCN records
II
· L
conventional atomic number of atom
MCONN(L)
coordination number of the atom with conventional atomic number L.
MCONN=0, coordination not tested
28
Border modification:
rec variable value
meaning
· NMO
number of border atoms to be modified
if NMO > 0 insert NMO records
II
· IPAD
label of the atom to be modified (cluster sequence)
NVIC
number of stars of neighbours of atom IPAD to be added to the cluster
IPAR
= 0
no hydrogen saturation
= 0
cluster border saturated with hydrogen atoms
BOND
bond length Hydrogen-IPAD atom (direction unchanged).
if NMO < 0 insert
II
· IMIN
label of the first atom to be saturated (cluster sequence)
IMAX
label of the last atom to be saturated (cluster sequence)
NVIC
number of stars of neighbours of each atom to be added to the cluster
IPAR
= 0
no hydrogen saturation
= 0
cluster border saturated with hydrogen atoms
BOND
H-cluster atom bond length (direction unchanged).
The two kinds of possible modification of the core cluster are (a) addition of further stars of
neighbours to specified border atoms, and (b) saturation of the border atoms with hydrogen.
This latter option can be essential in minimizing border electric field effects in calculations for
covalently-bonded systems.
(Substitution of atoms with hydrogen is obtained by HYDROSUB).
The hydrogen saturation procedure is carried out in the following way. First, a coordination
number for each atom is assumed (by default 4 for Si, 4 for Al and 2 for O, but these may
be modified in the input deck for any atomic number). The actual number of neighbours of
each specified border atom is then determined (according to a covalent radius criterion) and
compared with the assumed connectivity. If these two numbers differ, additional neighbours are
added. If these atoms are not neighbours of any other existing cluster atoms, they are converted
to hydrogen, otherwise further atoms are added until the connectivity allows complete hydrogen
saturation whilst maintaining correct coordination numbers.
The label of the IPAD atoms refers to the generated cluster, not to the original unit cell. The
preparation of the input thus requires two runs:
1. run using the CLUSTER option with NMO=0, in order to generate the sequence number
of the atoms in the core cluster. The keyword TESTGEOM should be inserted in the
geometry input block. Setting NNA = 0 in the input will print a coordination analysis of
all core cluster atoms, including all neighbours within a distance RNNA (which should
be set slightly greater than the maximum nearest neighbour bond length). This can be
useful in deciding what border modifications are necessary.
2. run using the CLUSTER option with NMO = 0, to perform desired border modifica-
tions.
Note that the standard CRYSTAL geometry editing options may also be used to modify the
cluster (for example by adding or deleting atoms) placing these keywords after the specification
of the CLUSTER input.
Warning. The system is 0D. No reciprocal lattice information is required in the scf input
(Section 1.4, page 19). See test 16.
COORPRT
Geometry information is printed: cell parameters, fractionary coordinates of all atoms in the
reference cell, symmetry operators.
A formatted file is written (in append mode) in fortran unit 33. See Appendix F, page 181.
No input data are required.
Fortran unit 33 has the right format for the program Xmol [26].
Download from http://biotech.icmb.utexas.edu/mime/xmol.html
29
or by the program MOLDEN [27] which can be downloaded from:
www.cmbi.kun.nl/ schaft/molden/molden.html
ELASTIC
An elastic deformation of the lattice may be defined in terms of the Z or
strain tensors defined
in section 6.9, page 147.
rec variable
value
meaning
· IDEF
±1
deformation through equation 6.36, Z matrix.
±2
deformation through equation 6.35:
matrix.
> 0
volume conserving deformation (equation 6.37).
< 0
not volume conserving (equation 6.36 or 6.35).
· D11 D12 D13
first row of the matrix.
· D21 D22 D23
second row of the matrix.
· D31 D32 D33
third row of the matrix.
The elastic constant is V -1 2E |
2
i =0 , where V
is the volume of the primitive unit cell.
i
The symmetry of the system is defined by the symmetry operators in the new crystallographic
cell. The keyword MAKESAED gives information on symmetry allowed elastic distortions.
The calculation of the elastic constants with CRYSTAL requires the following sequence of
steps:
1. select the ij matrix elements to be changed ( for example, 4 23 + 32), and set the
others j to zero;
2. perform calculations with different values of the selected matrix element(s) i: 0.02, 0.01,
0.001, -0.001, -0.01, -0.02, for example, and for each value compute the total energy E;
3. perform a polynomial fit of E as a function of i.
is adimensional, Z in °
A(default) or in bohr (page 24). The suggested value for IDEF is
-2 (deformation through equation 6.35, not volume conserving). The examples refer to this
setting.
Example
Geometry input deck to compute one of the energy points used for the evaluation of the C44
(page 150) elastic constants of Li2O [28].
CRYSTAL
0 0 0
3D code
225
3D space group number
4.5733
lattice parameter (°
A)
2
2 non equivalent atoms in the primitive cell
8 0.0 0.0 0.0
Z=8, Oxygen; x, y, z
3 .25 .25 .25
Z=3, Lithium; x, y, z
ATOMSYMM
printing of the point group at the atomic positions
ELASTIC
-2
deformation not volume conserving through equation 6.35
0. 0.03 0.03
matrix input by rows
0.03 0. 0.03
0.03 0.03 0.
ATOMSYMM
printing of the point group at the atomic positions after the defor-
mation
. . . . . . .
A rhombohedral deformation is obtained, through the
matrix. The printout gives information
on the crystallographic and the primitive cell, before and after the deformation:
LATTICE PARAMETERS (ANGSTROMS AND DEGREES) OF
(1) ORIGINAL PRIMITIVE CELL
30
(2) ORIGINAL CRYSTALLOGRAPHIC CELL
(3) DEFORMED PRIMITIVE CELL
(4) DEFORMED CRYSTALLOGRAPHIC CELL
A
B
C
ALPHA
BETA
GAMMA
VOLUME
(1)
3.233811
3.233811
3.233811
60.000000
60.000000
60.000000
23.912726
(2)
4.573300
4.573300
4.573300
90.000000
90.000000
90.000000
95.650903
(3)
3.333650
3.333650
3.333650
56.130247
56.130247
56.130247
23.849453
(4)
4.577414
4.577414
4.577414
86.514808
86.514808
86.514808
95.397811
After the deformation of the lattice, the point symmetry of the Li atoms is C3v, where the C3
axis is along the (x,x,x) direction. The Li atoms can be shifted along the principal diagonal,
direction (x,x,x) of the primitive cell without altering the point symmetry, as shown by the
printing of the point group symmetry obtained by the keyword ATOMSYMM (page 27).
See test 20 for complete input deck, including shift of the Li atoms.
END
Terminate processing of block 1, geometry definition, input. Execution continues. Subsequent
input records are processed, if required.
EXTPRT
A formatted input deck with explicit structural/symmetry information is written in fortran
unit 34. If the keyword is entered many times, the data are overwritten. The last geometry is
recorded. The deck may be used as input of the crystal geometry to CRYSTAL through the
EXTERNAL keyword (final optimized geometry, geometry obtained by editing who modified
the original space group). See Appendix F, page 181. No input data are required.
FRACTION
The keyword FRACTION means input coordinates given as fraction of the lattice parameter
in subsequent input, along the direction of translational symmetry:
x,y,z
crystals (3D)
x,y
slabs (2D; z in °
Angstrom or bohr)
x
polymers (1D; y,z in °
Angstrom or bohr)
no action for 0D.
When the unit of measure is modified, the new convention is active for all subsequent geometry
editing.
HYDROSUB - substitution with hydrogen atoms
rec variable
meaning
· NSOST
number of atoms to be substituted with hydrogen
insert NSOST records
II
· LA
label of the atom to substitute
LB
label of the atom linked to LA
BH
bond length B-Hydrogen
Selected atoms are substituted with hydrogens, and the bond length is modified. To be used
after CLUSTER.
31
KEEPSYMM
n any subsequent editing of the geometry, the program will endeavour to maintain the number
of symmetry operators, by requiring that atoms which are symmetry related remain so after
geometry editing (keywords:ATOMSUBS, ATOMINSE, ATOMDISP, ATOMREMO)
or the basis set (keywords CHEMOD, GHOSTS).
Example: When a CO molecule is vertically adsorbed on a (001) 3-layer MgO slab, (D4h
symmetry) (see page 25, keyword ATOMINSE), the symmetry is reduced to C4v, if the
BREAKSYM keyword is active. The symmetry operators related to the h plane are re-
moved. However, if KEEPSYMM is active, then additional atoms will be added to the
underside of the slab so as to maintain the h plane.
MAKESAED
This generates symmetry allowed elastic distortions. No input data are required.
MODISYMM
rec variable
meaning
· N
number of atoms to be attached a flag
· LA,LF(LA),L=1,N atom labels and flags (n couples of integers in 1 record).
The point symmetry of the lattice is lowered by attributing a different "flag" to atoms related
by geometrical symmetry. The symmetry operators linking the two atoms are removed and the
new symmetry of the system is analyzed. For instance, when studying spin-polarized systems, it
may be necessary to apply different spins to atoms which are related by geometrical symmetry.
MOLDRAW
A formatted input deck for the visualization program MOLDRAW [29] is written on fortran
unit 93. If the keyword is entered many times, the data are overwritten. The last geometry
can be visualized.
No input data are required. See:
http://www.moldraw.unito.it .
MOLEBSSE - counterpoise for molecular crystals
rec variable meaning
· NMOL
number of molecules to be isolated
insert NMOL records
· ISEED
label of one atom in the n-th molecule
J,K,L
integer coordinates (direct lattice) of the primitive cell containing the ISEED
atom
· NSTAR maximum number of stars of neighbours included in the calculation
RMAX
maximum distance explored searching the neighbours of the atoms belonging
to the molecule(s)
The counterpoise method is applied to correct the Basis Set Superposition Error in molecular
crystals. A molecular calculation is performed, with a basis set including the basis functions
of the selected molecules and the neighbouring atoms. The program automatically finds all
the atoms of the molecule(s) containing atom(s) ISEED (keyword MOLECULE, page 33).
The molecule is reconstructed on the basis of the covalent radii reported in Table on page 35.
32
They can be modified by running the option RAYCOV, if the reconstruction of the molecule
fails. The radius of the hydrogen atom is very critical when intermolecular hydrogen bonds
are present.
All the functions of the neighbouring atoms in the crystal are added to the basis set of the
selected molecule(s) such that both the following criteria are obeyed:
1. the atom is within a distance R lower than RMAX from at least one atom in the molecule
and
2. the atom is within the NSTAR-th nearest neighbours of at least one atom in the molecule.
Warning. The system obtained is 0D. No reciprocal lattice information is required in the scf
input (Section 1.4, page 19). See test 19.
MOLECULE - Extraction of n molecules from a molecular crystal
rec variable meaning
· NMOL
number of molecules to be isolated
insert NMOL records
II
· ISEED
label of one atom in the nth molecule
J,K,L
integer coordinates (direct lattice) of the primitive cell containing the
ISEED atom
The option MOLECULE isolates one (or more) molecules from a molecular crystal on the
basis of chemical connectivity, defined by the sum of the covalent radii (Table on page 35).
The program stops after printing full neighbouring analysis of the non-equivalent atoms, up to
n neighbours (default value 3; keyword NEIGHBOR, page 34 to modify it).
The input order of the atoms (atoms symmetry related are grouped) is modified, according
to the chemical connectivity. The same order of the atoms in the bulk crystal is obtained by
entering the keyword ATOMORDE (see Section 2.1, page 25). The total number of electrons
attributed to the molecule is the sum of the shell charges attributed in the basis set input (input
block 2, Section 1.2, page 16).
The keyword GAUSS98, entered in input block 2 (basis set input), writes an input deck to
run Gaussian 98 (see page 49)
Warning. The system is 0D. No reciprocal lattice information is required in the scf input
(Section 1.4, page 19).
Test 18 - Oxalic acid. In the 3D unit cell there are four water and two oxalic acid molecules.
The input of test 18 refers to a cluster containing a central oxalic acid molecule surrounded by
four water molecules.
MOLEXP - Variation of lattice parameters at constant symmetry
and molecular geometry
rec
variable
meaning
·
a,[b],[c], increments of the minimal set of crystallographic cell parameters:
[],[]
translation vectors length [°
Angstrom],
[]
crystallographic angles (degrees)
The cell parameters (the minimum set, see page 13) are modified, according to the increments
given in input. The volume of the cell is then modified. The symmetry of the lattice and the
geometry (bond lengths and bond angles) of the molecules within the cell is kept. The fractional
33
coordinates of the barycentre of the molecules are kept constant, the cartesian coordinates
redefined according to the modification of the lattice parameters. Optimization of the geometry
with reference to the compactness of the lattice is allowed, keeping constant the geometry of
the molecules. When there are very short hydrogen bonds linking the molecules in the lattice,
it may be necessary a modification of the atomic radii to allow proper identification of the
molecules (see option RAYCOV, page 35)
MOLSPLIT - Periodic lattice of non-interacting molecules
In order to compare bulk and molecular properties, it can be useful to build a density ma-
trix as a superposition of the density matrices of the isolated molecules, arranged in the same
geometry as in the crystal. The keyword MOLSPLIT (no additional input required) per-
forms an expansion of the lattice, in such a way that the molecules of the crystal are at an
"infinite" distance from each other. The crystal coordinates are scaled so that the distances
inside the molecule are fixed, and the distances among the molecules are expanded by a factor
100, to avoid molecule-molecule interactions. The 3D translational symmetry is not changed.
Reciprocal lattice information is required in the scf input (Section 1.4, page 19).
A standard wave function calculation of the expanded crystal is performed. The density matrix
refers to the non-interacting subsystems. Before running properties, the lattice is automatically
contracted to the bulk situation given in input. If a charge density or electrostatic potential
map is computed (ECHG, POTM options), it corresponds to the superposition of the charge
densities of the isolated molecules in the bulk geometry.
This option must be used only for molecular crystals only (no charged fragments).
See test 21.
NEIGHBOR/NEIGHPRT
rec
variable meaning
·
INEIGH number of neighbours of each non-equivalent atom to be printed
The option is active when analyzing the crystal structure (bond lengths and bond angles) and
when printing the bond populations following Mulliken analysis. Full input deck must be given
(block 1-2-3-4),in order to obtain neighbors analysis of all the non-equivalent atoms
For each non-equivalent atom information on the first INEIGH neighbours is printed: number,
type, distance, position (indices of the direct lattice cell).
Warning: the neighbors analysis is performed after the symmetry analysis and the screening
of the integrals. If very soft tolerances for the integrals screening are given in input, it may
happen that the information is not given for all the neighbors requested, as their are not taken
into account when truncation criteria are applied.
NOSHIFT
To be used before SUPERCEL keyword. It avoids shift of the origin in order to minimize the
number of symmetry operators with finite translation component. No input data are required.
ORIGIN
The origin is moved to minimize the number of symmetry operators with finite translation
components. Suggested before cutting a slab from a 3D structure (option SLABCUT, page
37) No input data are required.
34
PARAMPRT - printing of parametrized dimensions
The parameters controlling the dimensions of the static allocation arrays of the program are
printed.
No input data are required.
PRIMITIV
Some properties (XFAC, EMDL, EMDP, PROF) input the oblique coordinates of the k
points in the reciprocal lattice with reference to the conventional cell, though the computation
refers to the primitive one. This option allows entering directly the data with reference to the
primitive cell. The transformation matrix from primitive to crystallographic (Appendix A.5,
page 163) is set to the identity. No effect on the CPU time: CRYSTAL always refers to the
primitive cell. No input data are required.
PRINTOUT - Setting of printing environment
Extended printout can be obtained by entering selected keywords in a printing environment
beginning with the keyword PRINTOUT and ending with the keyword END. The possible
keywords are found in the fifth column of the table on page 179.
Extended printing request can be entered in any input block. Printing requests are not trans-
ferred from wave function to properties calculation.
See Appendix E, page 177.
PRSYMDIR
Printing of displacement directions allowed by symmetry. The printing is done after the neigh-
bor analysis, before computing the wave function. Full input must be supplied (4 blocks). Test
run allowed with the keyword TESTPDIM.
No input data required.
PURIFY
This cleans up the atomic positions so that they are fully consistent with the group (to within
machine rounding error). No input data are required.
RAYCOV - covalent radii modification
rec variable meaning
· NCOV
number of atoms for which the covalent radius is redefined
insert NCOV records
II
· NAT
atomic number (0 NAT 92)
RAY
covalent radius of the atom with atomic number NAT ([°
A], default,
or bohr, if the keyword BOHR precedes in the deck)
The option RAYCOV allows modification of the covalent radius default value for a given
atom.
35
Table of covalent radii (Angstrom)
H
He
0.68
1.47
---------
-----------------------------
Li
Be
B
C
N
O
F
Ne
1.65 1.18
0.93 0.81 0.78 0.78 0.76 1.68
---------
-----------------------------
Na
Mg
Al
Si
P
S
Cl
Ar
2.01 1.57
1.50 1.23 1.15 1.09 1.05 0.97
-----------------------------------------------------------------------------------------
K
Ca
Sc
Ti
V
Cr
Mn
Fe
Co
Ni
Cu
Zn
Ga
Ge
As
Se
Br
Kr
2.31 2.07 1.68 1.47 1.41 1.47 1.47 1.47 1.41 1.41 1.41 1.41 1.36 1.31 1.21 1.21 1.21 2.10
-----------------------------------------------------------------------------------------
Rb
Sr
Y
Zr
Ni
Mo
Tc
Ru
Rh
Pd
Ag
Cd
In
Sn
Sb
Te
I
Xe
2.31 2.10 1.94 1.60 1.52 1.52 1.42 1.36 1.42 1.47 1.68 1.62 1.62 1.52 1.52 1.47 1.47 2.66
-----------------------------------------------------------------------------------------
Cs
Ba
La
Hf
Ta
W
Re
Os
Ir
Pt
Au
Hg
Tl
Pb
Bi
Po
At
Rn
2.73 2.10 1.94 1.63 1.63 1.63 1.63 1.63 1.63 1.63 1.63 1.63 1.99 1.89 1.68 1.42 1.42 1.62
-----------------------------------------------------------------------------------------
The choice of the covalent radius of hydrogen may be very critical when extracting a molecule
from a hydrogen bonded molecular crystal. See test 15.
REDEFINE - 3D unit cell redefinition
rec variable meaning
· h,k,l
Crystallographic (Miller) indices of the basal layer of the new 3D unit cell
1. A new unit cell is defined, with two lattice vectors perpendicular to the [hkl] direction.
The indices refer to the Bravais lattice of the crystal; the hexagonal lattice is used for
the rhombohedral systems, the cubic lattice for cubic systems (non primitive).
2. A new Cartesian reference system is defined, with the xy plane parallel to the (hkl) plane.
3. The atoms in the reference cell are re-ordered according to their z coordinate, in order
to recognize the layered structure, parallel to the (hkl) plane.
4. The layers of atoms are numbered. This information is necessary for generating the input
data for the SLAB option.
5. After neighboring analysis, the program stops.
6. the keyword ORIGIN can be used to shift the origin after the rotation of the cell, and
minimize the number of symmetry operators with translational component. Useful to
maximize the point group of the 2D system that can be generated from 3D using the
keyword SLAB (page 37).
SETINF - Setting of INF values
rec variable
meaning
· NUM
number of INF vector positions to set
· J,INF(J),I=1,NUM
position in the vector and corresponding value
The keyword SETINF allows setting of a value in the INF array. It can be entered in any
input section.
36
SETPRINT - Setting of printing options
rec variable
meaning
· NPR
number of LPRINT vector positions to set
· J,LPRINT(J),I=1,NPR
prtrec; position in the vector and corresponding value
The keyword SETPRINT allows setting of a value in the LPRINT array, according to the
information given in Appendix E, page 179. It can be entered in any input section.
SLABCUT (SLAB)
rec variable meaning
· h, k, l
crystallographic (Miller) indices of the plane parallel to the surface
· ISUP
label of the surface layer
NL
number of atomic layers in the slab
The SLABCUT option is used to create a slab of given thickness, parallel to the given plane
of the 3D lattice.
A new Cartesian frame, with the z axis orthogonal to the (hkl) plane, is defined. A layer is
defined by a set of atoms with same z coordinate, with reference to the new Cartesian frame.
The thickness of the slab, the 2D system, is defined by the number of layers. No reference is
made to the chemical units in the slab. The neutrality of the slab is checked by the program.
1. The crystallographic (Miller) indices of the plane refer to the conventional cell (cubic and
hexagonal systems).
2. A two-sided layer group is derived from the 3D symmetry group of the original crystal
structure: the origin may be shifted to maximize the order of the layer group (keyword
ORIGIN, page 34).
3. The unit cell is selected with upper and lower surface parallel to the (hkl) plane.
4. The 2D translation vectors a1 and a2 are chosen according to the following criteria:
(a) minimal cell area;
(b) shortest translation vectors;
(c) minimum |cos()|, where is the angle between a1 and a2.
5. The surface layer ISUP may be found from an analysis of the information printed by
the REDEFINE (page 36) option. This information can be obtained by a test run,
inserting in the geometry input block the keyword TESTGEOM (page 40). Only the
geometry input block is processed, then the program stops.
Two separate runs are required in order to get the information to prepare the input for a full
SLAB option run:
1. keyword REDEFINE: Rotation of the 3D cell, to have the z axis perpendicular to the
(hkl) place, with numbering of the atomic layers in the rotated reference cell, according
to their z coordinate of the atoms (insert STOP after REDEFINE to avoid further
processing).
2. keyword SLAB: Definition of the 2D system, a slab of given thickness (NL, number of
atomic layers) parallel to the (hkl) crystallographic plane, with the ISUP-th atom on the
surface layer
37
The SLABCUT option, combined with ATOMINSE (page 25), ATOMDISP (page 25),
etc. can be used to create a slab of given thickness, with an atom (or group of atoms) adsorbed
at given position. This is achieved by adding new atoms to the 2D structure, obtained after
executing the SLAB option.
Test cases 5-6-7 refer to a 2D system; test cases 25-26-27 refer to the same system, but generated
from the related 3D one. See also tests 35, 36, 37.
STOP
Execution stops immediately. Subsequent input records are not processed.
STRUCPRT
A formatted deck with cell parameters and atoms coordinates in cartesian reference is written
in the file STRUC.INCOOR. See appendix F.
SUPERCEL
rec variable meaning
· E
expansion matrix E (IDIMxIDIM elements, input by rows: 9 reals (3D); 4 reals
(2D); 1 real (1D)
A supercell is obtained by defining the new unit cell vectors as linear combinations of the
primitive cell unit vectors (use SUPERCON for conventional cell vectors reference). The
point symmetry is defined by the number of symmetry operators in the new cell. It may be
reduced, not increased.
The new translation vectors b , b , b are defined in terms of the old vectors b
1
2
3
1, b2, b3 and of
the matrix E, read in input by rows, as follows:
b =
e
1
11 · b1
+
e12 · b2
+
e13 · b3
b =
e
2
21 · b1
+
e22 · b2
+
e23 · b3
b =
e
3
31 · b1
+
e32 · b2
+
e33 · b3
The symmetry is automatically reduced to the point symmetry operators without translational
components and a further reduction of the symmetry is also possible.
Before building the supercell, the origin is shifted in order to minimize the number of sym-
metry operators with translational components (see page 14). To avoid this operation, insert
NOSHIFT before SUPERCEL
Atoms that are related by translational symmetry in the unit cell are considered inequivalent
in a supercell.
The supercell option is a useful starting point for the study of defective systems, of chemisorp-
tion and anti ferromagnetism, by combining the SUPERCELoption with the options de-
scribed in this chapter: ATOMREMO (page 25), ATOMSUBS (page 27), ATOMINSE
(page 25), ATOMDISP (page 25), SLAB (page 37).
To study anti ferromagnetic (AFM) states, it may be necessary to generate a supercell, and
then attribute different spin to atoms related by translational symmetry (ATOMSPIN, input
block 4, page 73). See tests 17, 30, 31, 34, 37.
Example.
Construction of supercells of face-centred cubic 3D system (a = 5.42 °
A).
The crystallographic cell is non-primitive, the expansion matrix refers to primitive
cell vectors. The E matrix has 9 elements:
38
PRIMITIVE CELL
DIRECT LATTICE VECTORS COMPONENTS
X
Y
Z
B1
.000
2.710
2.710
B2
2.710
.000
2.710
B3
2.710
2.710
.000
2 UNITS SUPERCELL (a)
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
.000
1.000
1.000
B1
5.420
2.710
2.710
E2
1.000
.000
1.000
B2
2.710
5.420
2.710
E3
1.000
1.000
.000
B3
2.710
2.710
5.420
2 UNITS SUPERCELL (b)
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
1.000
1.000
-1.000
B1
.000
.000
5.420
E2
.000
.000
1.000
B2
2.710
2.710
.000
E3
1.000
-1.000
.000
B3
-2.710
2.710
.000
4 UNITS SUPERCELL (c) crystallographic cell
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
-1.000
1.000
1.000
B1
5.420
.000
.000
E2
1.000
-1.000
1.000
B2
.000
5.420
.000
E3
1.000
1.000
-1.000
B3
.000
.000
5.420
8 UNITS SUPERCELL
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
2.000
.000
.000
B1
.000
5.420
5.420
E2
.000
2.000
.000
B2
5.420
.000
5.420
E3
.000
.000
2.000
B3
5.420
5.420
.000
16 UNITS SUPERCELL
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
3.000
-1.000
-1.000
B1
-5.420
5.420
5.420
E2
-1.000
3.000
-1.000
B2
5.420
-5.420
5.420
E3
-1.000
-1.000
3.000
B3
5.420
5.420
-5.420
27 UNITS SUPERCELL
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
3.000
.000
.000
B1
.000
8.130
8.130
E2
.000
3.000
.000
B2
8.130
.000
8.130
E3
.000
.000
3.000
B3
8.130
8.130
.000
32 UNITS SUPERCELL
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
-2.000
2.000
2.000
B1
10.840
.000
.000
E2
2.000
-2.000
2.000
B2
.000
10.840
.000
E3
2.000
2.000
-2.000
B3
.000
.000
10.840
a), b) Different double cells
c) quadruple cell. It corresponds to the crystallographic, non-primitive cell, whose parameters
are given in input (page 14).
Example.
Construction of supercells of hexagonal R¯
3 (corundum lattice) cubic 3D system.
The crystallographic cell is non-primitive: CRYSTAL refer to the primitive cell, with volume
1/3 of the conventional one. The E matrix has 9 elements:
GEOMETRY INPUT DATA:
LATTICE PARAMETERS
(ANGSTROMS AND DEGREES) - CONVENTIONAL CELL
A
B
C
ALPHA
BETA
GAMMA
4.76020
4.76020
12.99330
90.00000
90.00000
120.00000
TRANSFORMATION WITHIN CRYSTAL CODE FROM CONVENTIONAL TO PRIMITIVE CELL:
LATTICE PARAMETERS
(ANGSTROMS AND DEGREES) - PRIMITIVE CELL
A
B
C
ALPHA
BETA
GAMMA
VOLUME
5.12948
5.12948
5.12948
55.29155
55.29155
55.29155
84.99223
3 UNITS SUPERCELL crystallographic cell
EXPANSION MATRIX
DIRECT LATTICE VECTORS
E1
1.000
-1.000
.000
B1
4.122
-2.380
.000
E2
.000
1.000
-1.000
B2
.000
4.760
.000
E3
1.000
1.000
1.000
B3
.000
.000
12.993
LATTICE PARAMETERS (ANGSTROM AND DEGREES)
39
A
B
C
ALPHA
BETA
GAMMA
VOLUME
4.76020
4.76020 12.99330
90.000
90.000
120.000
254.97670
SUPERCON
A supercell is obtained by defining the new unit cell vectors as linear combinations of the
conventional cell vectors. The point symmetry is defined by the number of symmetry operators
in the new cell. It may be reduced, not increased. See SUPERCEL, page 38 for input
instructions.
SYMMOPS
Point symmetry operator matrices are printed in the Cartesian representation. No input data
are required.
SYMMREMO
All the point group symmetry operators are removed. Only the identity operator is left. The
wave function can be computed. No input data are required.
Warning: the CPU time may increase by a factor MVF (order of point-group), both in the
integral calculation and in the scf step. The size of the bielectronic integral file may increase
by a factor MVF2.
SYMMDIR
The symmetry allowed directions, corresponding to internal degrees of freedom are printed.
No input data are required.
TENSOR
rec
variable meaning
·
IORD
order of the tensor ( 4)
This option evaluates and prints the non zero elements of the tensor of physical properties up
to order 4.
TESTGEOM
Execution stops after reading the geometry input block and printing the coordinates of the
atoms in the conventional cell. Neighbours analysis, as requested by the keyword NEIGH-
BOR, is not executed. The geometry input block must end with the keyword END or ENDG.
No other input blocks (basis set etc) are required.
TRASREMO
Point symmetry operators with fractional translation components are removed. It is suggested
to previously add the keyword ORIGIN (page 34), in order to minimize the number of sym-
metry operators with finite translation component. No input data are required.
40
USESAED
rec variable
meaning
· (i),i=1,nsaed for each distortion
Given the symmetry allowed elastic distortion (SAED), (printed by the keyword MAKE-
SAED, page 32) for the allowed distortion are given in input.
FIELD - Electric field in the hamiltonian
rec variable value
meaning
· E0MAX
electric field intensity E0 (in atomic units)
if CRYSTAL 3D calculation insert
II
· MUL
number of term in Fourier expansion for triangular electric potential
ISYM
+1
triangular potential is symmetric with respect to the z = 0 plane
-1
triangular potential is anti-symmetric with respect to the z = 0 plane
A perturbation due to a finite periodic electric field (E) is added to the Hamiltonian (Fock or
Kohn-Sham):
^
H = ^
H0 + ^
H1(E)
(2.1)
where ^
H0 is the unperturbed Hamiltonian and ^
H1(E) the electric-dipole interaction.
The crystalline orbitals are relaxed under the effect of the field, leading to a perturbed wave
function and charge density.
In the three-dimensional case, the applied electric field is periodic, in order to maintain the
translational symmetry. The corresponding electric potential has a triangular form ("saw-
tooth").
When the field is along the z direction, the perturbation ^
H1(E) takes the form:
^
H±(E
1
z ) = V (z) = -qE0 × f ±(z)
(2.2)
where the f + or the f - function is developed in a Fourier series and is chosen according to
the symmetry of the supercell in the direction of the applied field:
+
2C
1
2(2k + 1)z
f +(z) =
cos
(2.3)
2
(2k + 1)2
C
k=0
+
2C
(-1)k
2(2k + 1)z
f -(z) =
sin
(2.4)
2
(2k + 1)2
C
k=0
C is the cell size in the field direction.
Not implemented fo 0-D (molecules) and 1-D (polymers) systems.
1. A supercell approach (see keywords SUPERCEL/SUPERCON, page 38 allows to
minimize edge effect and obtain a better convergence of results.
2. The direction of electric field is along z axis, In 2D systems (slab) it is always perpendic-
ular to the slab. The symmetry of the system can be reduced due to the anti-symmetric
nature of the field, even if intensity E0 = 0.
41
3. In 3D systems, as the direction of electric field is always along the z axis, the keyword
ROTATE (page 36 allows definition of a new crystallographic cell, where in the new
Cartesian reference system the z axis is perpendicular to a selected (hkl) plane.
4. In 3D-crystal, the electric potential takes a triangular form to maintain translational
symmetry and electric neutrality of cell. The symmetry of triangular potential has two
options:
a) ISYM=+1, triangular potential is symmetric with respect to the center of supercell,
along z axis. Use this option if there is a symmetry plane orthogonal to z axis.
b) ISYM=-1, triangular potential is anti-symmetric. This option can be used when
the supercell does not have a symmetry plane orthogonal to z axis.
5. MUL, the number of term in Fourier expansion, can take values between 1 and 60
(LIM060).
MUL=40 is sufficient to perfectly reproduce the triangular shape of the potential.
Note - CPU time is not MUL-dependent.
6. When |E0| takes high value, the system can reach a non converged conducting state during
the SCF process. The threshold |Emax
0
| value depends on the dielectric susceptibility
of the system and on the gap width. For very narrow gap cases, the eigenvalue level
shifting technique (keyword LEVSHIFT, page 77) can be useful to avoid the SCF fall
in a spurious conducting state
7. When an external field is applied, the system can become conducting during the SCF
procedure. In order to avoid convergence problems, it is advisable to set the shrinking
factor of the Gilat net ISP equal to 2× IS, where IS is the Monkhorst net shrinking factor
(see SCF input, page 19.
Conversion factors for electric field
1 AU = 1.71527E+07 ESU·CM-2 = 5.72152E+01 C ·M-2 = 5.14226E+11 V·M-1
42
OPTCOORD - Atom coordinates optimization
A modified conjugate gradient algorithm (H.B. Schlegel, J. Comp. Chem. 3, 214, 1982) is
used [30] to optimize the atomic coordinates and locate minima on the potential energy surface
(PES). A cell fixed unconstrained optimization of the atomic coordinates can be carried out.
Lattice parameters must be optimized numerically point by point.
HF and DFT (pure and hybrid functionals) analytical gradients are used for insulators, both
for all-electron and ECP calculations. For conducting systems numerical first derivatives must
be computed (keyword NUMGRAD).
Geometry optimization can be performed on systems with any periodicity: molecules, poly-
mers, slabs, and crystals.
Gradients are evaluated every time the total energy is computed; the second derivative matrix
is built from the gradients. At each step, a one dimensional minimization using a quadratic
polynomial is carried out, followed by an n-dimensional search using the Hessian matrix.
The OPTCOORD keyword must be the last keyword in the geometry input section. The
optimization input block must be closed by the keyword END (or ENDOPT).
No input data are required, as default values for all parameters controlling the optimization
are supplied. The default value for SCF convergence criteria on total energy is set to 10-7.
Geometry optimization is performed in symmetrized cartesian coordinates, in order to exploit
the point group symmetry of the lattice. The keyword PRSYMDIR (input block 1, page
35), may be used to print the symmetrized directions adopted in the geometry optimization.
If there are no symmetry allowed directions, the program prints a warning message and stops.
The output from the wave function and gradient calculation is printed at the first step only.
The output from the other steps is routed to fortran unit 65.
The fortran unit 33 contains the cartesian coordinates of the atoms in the unit cell for each
geometry optimization step in a simple xyz format (see keyword COORPRT, input block1 1,
page 29). This file is suitable to be read by molecular graphics programs (e.g. Molden, XMol,
...) to display the animation of the geometry optimization run.
Fortran unit 34 contains complete geometry input deck, with the last defined geometry. To
run crystal with that geometry, insert EXTERNAL keyword to define the geometry in the
input stream (see keyword EXTPRT, input block 1, page 31, for explanation of the format).
Fortran unit 68 contains information to restart optimization. (see keyword RESTART in
OPTCOORD input block).
ATOMFREE
Partial geometry optimization (default: global optimization)
· NL
number of atoms "free"
· LB(L),L=1,NL
label of the atoms to move
This keyword allows partial geometry optimization, limited to an atomic fragment rather than
the whole system. Symmetrized cartesian coordinates are generated according to the list of
atoms allowed to move. Note that no advantage is taken in the gradient calculation to reduce
the number of atoms, i.e. gradients are calculated on the whole system. The symmetrized
forces are then computed by using the new set of symmetrized coordinates.
43
FINALRUN
action after geometry optimization - integrals classification is based on the
last geometry
· ICODE
Action type code:
0
the program stops (default)
1
single-point energy calculation
2
single-point energy and gradient calculation
3
single-point energy and gradient calculation - if convergence criteria on
gradients are not satisfied, optimization restarts
Truncation of infinite series (Coulomb, exchange, etc.), based on the overlap between two
atomic functions (see chapter 6.11), depends on the geometry of a crystal.
With default
setting of thresholds slightly different truncation levels correspond to different geometries that
often introduce small discontinuities in the PES, producing artificial noise in the optimization
process. To avoid noise in interpolation of PES, FIXINDEX option is always active during
optimization. The adopted truncation pattern refers to the starting geometry.
If equilibrium geometry is significantly different from starting point, reference truncation pat-
tern may be inappropriate and the use of proper truncation becomes mandatory.
Since both total energy and gradients are affected by the integrals classification, a single-point
energy calculation ought to be run always with the final structure, and integrals classified
according to the new geometry, to calculate correct total energy and gradients.
If the convergence test on the forces is not satisfied, optimization has to be restarted, keeping
the integrals classification based on the new geometry.
The three different options of FINALRUN allow the following actions, after classification of
integrals:
1. single-point energy calculation (correct total energy),
2. single-point energy and gradient calculation (correct total energy and gradients),
3. single-point energy and gradient computation, followed by a new optimization process,
starting from the final geometry of the previous one, (used to classify the integrals), if
the convergence test is not satisfied.
If the starting and final geometry are not far away, the energy and gradient calculated from
the final geometry, but the integrals classification based on the initial geometry, and the values
computed with integrals classification based on the final geometry are not very different. In
some cases (e.g. optimization of the geometry of a surface, with reconstruction) the two ge-
ometries are very different, and a second optimization cycle is almost mandatory (ICODE=3).
HESGUESS
defines the initial guess for the Hessian
· ICODE
Initial guess code:
-1
identity matrix (default)
0
numerical estimate (two-points formula)
1
external guess (read from fort.66
By default a unit matrix is adopted as initial Hessian. The Hessian matrix is stored in fortran
unit 66 during the optimization at each step. This may be useful to restart the optimization
from a previous run performed at a lower level of theory (e.g. a smaller basis set).
An initial Hessian can also be obtained as numerical first-derivative, but this process is very
expensive.
MAXOPTC
· MAX
maximum number of optimization steps (default 100)
44
NOGUESS scf guess at each geometry point: superposition of atomic densities at each
scf calculation
At each geometry point the default guess for scf is the density matrix calculated at the end
of the previous run. If the solution does not correpond at real convergence, but to an en-
ergy stabilization due to the techiques applied to help convergence (LEVSHIFT, FMIXING,
BROYDEN..) when the final density matrix, obtained from the last cycle eigenvectors, is
used to build the first hamiltonian matrix, the hamiltonian eigenvalues can be unphysical,
no chances to recover the scf process. In those cases it may be better an atomic guess.
NUMGRAD numerical first-derivatives are used during the geometry optimization
The nuclear coordinate gradients of the energy can also be computed numerically. A three-point
numerical derivative formula is adopted. A finite positive (and then negative) displacement
is applied to the desired coordinate and a full SCF calculation is performed. The gradient is
then computed as
Ex - E
g
i
-xi
i =
2 xi
where xi is the finite displacement along the i-coordinate.
Such a computation is very expensive compared to analytical gradients, since the cost is 2·N ·t,
where N is the number of coordinates to be optimized and t the cost of the SCF calculation. Nu-
merical first-derivatives should be avoided whenever possible, but sometimes they are the only
way to obtain gradients (i.e. for metals) and therefore to optimize the atom coordinates.
PRINTFORCES printing atomic gradients
PRINTHESS
printing Hessian inormation
PRINTOPT
prints information on optimization process
PRINT
verbose printing
RESTART restart geometry optimization from a previous run
Full restart of geometry optimization is possible for a job which is abruptly terminated (e.g.
number of steps exceeded, available cpu time exceeded,...). The optimization restarts from the
last complete step of the previous run. Needed information is read from fortran unit 68, and
saved in the same file at each step .
The same input deck as for the initial geometry optimization must be used when the RESTART
keyword is added.
Convergence criteria
A stationary point on the potential energy surface is found when the forces acting on atoms
are zero. So, a geometry optimization is usually completed when the gradients are below a
given threshold.
In Crystal, the optimization convergence is checked on the root-mean-square and the absolute
value of the largest component of both the gradients and the estimated displacements. When
these four conditions are all satisfied at a time, optimization is considered complete.
In the current implementation the following defaults are used (in a.u.):
45
TOLDEG
RMS GRADIENT
0.000300 default value - MAX GRADIENT 1.5*TOLDEG
TOLDEX
RMS DISPLAC.
0.001200 default value - MAX DISPLAC. 1.5*TOLDEX
The maximum gradient and the maximum displacement threshold are defined as 1.5 times
TOLDEG and TOLDEX, respectively.
In case of flat surfaces (e.g. a system formed by two weakly-interacting moieties, a molecule
adsorbed on a surface at large distance) it may happen that the convergence criteria on dis-
placements are not satisfied, though the energy does not change. Additional combined test on
gradient and energy has been adopted:
1. If the gradient criteria are satisfied (but not the displacement criteria) and the energy
difference between two steps is below a given threshold (see TOLDEE), the optimization
stops with a warning message;
2. If both the gradient and displacements criteria are not satisfied, but the energy does not
change (TOLDEE parameter) for four subsequent steps, the optimization stops with a
warning message.
TOLDEE
threshold on the energy change between optimization steps
· IG
|E| < 10-IG (default: 7)
The value of IG must be larger or equal to the threshold adopted for the SCF convergence.
The value is checked when input block 4, defining the SCF convergence criteria, is processed.
TOLDEG
convergence criterion on the RMS of the gradient
· TG
max RMS of the gradient (default: 0.0003)
TOLDEX
convergence criterion on the RMS of the displacement
· TX
max RMS of the displacement (default: 0.0012)
Full optimization: cell parameters and atomic coordinates
A full optimization of a crystalline structure includes also the optimization of the lattice pa-
rameters. However, analytical gradients for the cell constants are not yet implemented; thus,
a two-step iterative process must be adopted:
1. the cell parameters are optimized at fixed atomic positions. This must be accomplished
by hand or by means of external drivers (e.g. the LoptCG script).
2. atomic positions are relaxed by keeping the lattice parameters fixed at the previously
optimum values.
The process is then iterated until convergence is attained, that is cell constants no longer
change and convergence criteria on nuclear forces are satisfied.
46
Notes - From periodic structure to molecules and clusters
The geometry editing described in this section allow the generation of finite (non-periodic)
systems, derived from periodic structures: clusters, molecules, groups of molecules (for exam-
ple, dimers), functional groups. The atoms in a molecular crystal are assigned to the same
molecule when their distance is shorter than the sum of the covalent radii, according to table
at page 35.
The possibility of computing the wave function of a finite cluster allows the comparison of
data with those computed for the original periodic structure. In this way the exact nature of
the 'border effect' introduced by the cluster approximation in solid-state calculations may be
analyzed [31, 32].
The default value of covalent radii can be modified for selected atoms by the option RAYCOV
(page 35). This modification may be necessary to force the correct definition of a fragment in
the MOLECULE and CLUSTER options, as well as in the MOLEBSSE option (page 47).
Warning
The result of the geometry editing:
CLUSTER, MOLECULE, ATOMBSSE and
MOLEBSSE is a 0D system. No reciprocal lattice information is required in the scf in-
put (Section 1.4, page 19).
ATOMORDE can be applied to molecular crystals only. The molecules ordered in the lattice
are defined as groups of atoms with a separation less than the sum of their covalent radii (Table
on page 35). When there are very short hydrogen bonds linking the molecules in the lattice,
it may be necessary a modification of the atomic radii, to allow proper identification of the
molecules (see option RAYCOV, page 35)
Notes - From 3D periodic structure to 2D
Slab model for surfaces
to be written
at least:
test on thickness of the slab
non-polar surfaces,
calculation of surface energy
Notes - BSSE correction - counterpoise scheme
The counterpoise method [33] is applied to correct the Basis Set Superposition Error for atoms,
ions or molecules. An atom or cluster of atoms surrounded by the basis functions of the
neighbours are extracted from the periodic structure.
A BSSE estimate, maintaining the periodic symmetry of the system, can be obtained trans-
forming into ghosts selected atoms in the lattice (keyword GHOSTS, page 50, test 36).
example: urea bulk
47
2.2
Basis set input
Symmetry control
ATOMSYMM
printing of point symmetry at the atomic positions
27
Basis set modification
CHEMOD
modification of the electronic configuration
48
I
GHOSTS
eliminates nuclei and electrons, leaving BS
50
I
Auxiliary and control keywords
CHARGED
allows non-neutral cell
48
PARAMPRT
printing of parameters controlling code dimensions
35
PRINTOUT
setting of printing options
35
I
SETINF
setting of inf array options
36
I
SETPRINT
setting of printing options
37
I
STOP
execution stops immediately
38
SYMMOPS
printing of point symmetry operators
40
END/ENDB
terminate processing of basis set definition keywords
Output of data on external units
GAUSS98
printing of an input file for the GAUSS94/98 package
49
ATOMSYMM
See input block 1, page 27
CHARGED - charged reference cell
The unit cell of a periodic system must be neutral. This option forces the overall system to
be neutral even when the number of electrons in the reference cell is different from the sum
of nuclear charges, by adding a uniform background charge density to neutralize the charge in
the reference cell.
CHEMOD - modification of electronic configuration
rec variable
meaning
· NC
number of configurations to modify
· LA
label of the atom with new configuration
CH(L),L=1,NS shell charges of the LA-th atom. The number NS of shells must coincide
with that defined in the basis set input.
The CHEMOD keyword allows modifications of the shell charges given in the basis set input,
which are used in the atomic wave function routines. The original geometric symmetry is
checked, taking the new electronic configuration of the atoms into account. If the number of
symmetry operators should be reduced, information on the new symmetry is printed, and the
program stops. No automatic reduction of the symmetry is allowed. Using the information
printed, the symmetry must be reduced by the keyword MODISYMM (input block 1, page
32).
See test 37. MgO supercell, with a Li defect. The electronic configuration of the oxygen nearest
to Li corresponds to O-, while the electronic configuration of those in bulk MgO is O2-. The
48
basis set of oxygen is unique, while the contribution of the two types of oxygen to the initial
density matrix is different.
END
Terminate processing of block 2, basis set, input. Execution continues. Subsequent input
records are processed, if required.
GAUSS98 - Printing of input file for GAUSS98 package
The keyword GAUSS98 writes on fortran unit 92 an input deck to run Gaussian 94 (or
Gaussian 98) [34, 35]. The deck can be prepared without the calculation of the wave function
by entering the keyword TESTPDIM in input block 3 (page 69). For periodic systems,
coordinates and basis set for all the atoms in the reference cell are written.
If the keyword is entered many times, the data are overwritten. The fortran file fort.92 contains
the data corresponding to the last call.
No input data required.
1. The route card specifies:
method
HF
basis set
GEN
type of job
SP
geometry
UNITS=AU
2. The title card is the same as in CRYSTAL input.
3. The molecule specification defines the molecular charge as the net charge in the reference
cell. If the system is not closed shell, the spin multiplicity is indicated with a string "??",
and must be defined by the user.
4. Input for effective core pseudopotentials is not written. In the route card PSEUDO =
CARDS is specified; the pseudopotential parameters used for the crystal calculation are
printed in the crystal output.
5. The scale factors of the exponents are all set to 1., as the exponents are already scaled.
6. the input must be edited when different basis sets are used for atoms with the same
atomic number (e.g., CO on MgO, when the Oxygen basis set is different in CO and in
MgO)
To compare molecular energies obtained with GAUSSIAN and CRYSTAL, default computa-
tional parameters of CRYSTAL must be modife by inserting the keywords:
NOBIPOLA remove bipolar exoansio to compute coulomb interaction of non overlapping
pseudo-charge distributions;
BOHRANGS follwed by the value adopteb by GAUSSIAN for the conversion factor bohr/
°
Angstrom
Warning: Only for 0D systems! The programs does not stop when the keyword GAUSS94
is entered for 1-2-3D systems. Coordinates and basis set of all the atoms in the primitive cell
are written, formatted, on fortran unit 92, following Gaussian 94 scheme.
49
GHOSTS
rec variable
meaning
· NA
number of atoms to be transformed into ghosts
· LA(L),L=1,NA
label of the atoms to be transformed.
Selected atoms may be transformed into ghosts, by deleting the nuclear charge and the shell
electron charges, but leaving the basis set centred at the atomic position. The conventional
atomic number is set to zero.
If the system is forced to maintain the original symmetry (KEEPSYMM), all the atoms
symmetry related to the given one are transformed into ghosts.
Useful to create a vacancy (Test 37), leaving the variational freedom to the defective region
and to evaluate the basis set superposition error (BSSE), in a periodic system. The periodic
structure is maintained, and the energy of the isolated components computed, leaving the basis
set of the other one(s) unaltered. For instance, the energy of a mono-layer of CO molecules on
top of a MgO surface can be evaluated including the basis functions of the first layer of MgO,
or, vice-versa, the energy of the MgO slab including the CO ad-layer basis functions. See test
36.
Warning Do not use with ECP
Warning The keyword ATOMREMO (input block 1, page 25) creates a vacancy, removing
nuclear charge, electron charge, and basis functions. The keyword GHOSTS creates a vacancy,
but leaves the basis functions at the site, so allowing better description of the electron density
in the vacancy.
PARAMPRT - Printing of parametrized dimensions
See input block 1, page 35.
PRINTOUT - Setting of printing environment
See input block 1, page 35.
SETINF - Setting of INF values
See input block 1, page 36.
SETPRINT - Setting of printing options
See input block 1, page 37.
STOP
Execution stops immediately. Subsequent input records are not processed.
SYMMOPS
See input block 1, page 40
50
Effective core pseudo-potentials.
rec
variable value
meaning
· A
PSN
pseudo-potential keyword:
HAYWLC
Hay and Wadt large core ECP.
HAYWSC
Hay and Wadt small core ECP.
BARTHE
Durand and Barthelat ECP.
DURAND
Durand and Barthelat ECP.
INPUT
free ECP - input follows.
if PSN = INPUT insert
II
·
ZNUC
effective core charge (ZN in eq. 2.6).
M
Number of terms in eq. 2.7
M0
Number of terms in eq. 2.8 for
=0.
M1
Number of terms in eq. 2.8 for
=1.
M2
Number of terms in eq. 2.8 for
=2.
M3
Number of terms in eq. 2.8 for
=3.
insert M+M0+M1+M2+M3 records
II
·
ALFKL
Exponents of the Gaussians: k .
CGKL
Coefficient of the Gaussians: Ck .
NKL
Exponent of the r factors: nk .
Valence-electron only calculations can be performed with the aid of effective core pseudo-
potentials (ECP). The ECP input must be inserted into the basis set input of the atoms with
conventional atomic number > 200.
The form of pseudo-potential Wps implemented in CRYSTAL is a sum of three terms: a
Coulomb term (C), a local term (W0) and a semi-local term (SL):
Wps = C + W 0 + SL
(2.5)
where:
C = -ZN /r
(2.6)
M
W 0 =
rnk Cke-kr2
(2.7)
k=1
3
M
SL =
[
rnk Ck e-k r2 ]P
(2.8)
=0 k=1
ZN is the effective nuclear charge, equal to total nuclear charge minus the number of electrons
represented by the ECP, P is the projection operator related to the
angular quantum number,
and M, nk, k, M , nk , Ck , k are atomic pseudo-potential parameters.
1. Hay and Wadt (HW) ECP ([36, 37]) are of the general form 2.5. In this case, the NKL
value given in the tables of ref. [36, 37] must be decreased by 2 (2 0, 1 -1, 0 -2).
2. Durand and Barthelat (DB) ([38] - [39], [40], [41]), and Stuttgart-Dresden [42] ECPs
contain only the Coulomb term C and the semi-local SL term.
3. In Durand and Barthelat ECP the exponential coefficient in SL depends only on
(i.e.
it is the same for all the Mk terms).
3
M
SL =
e- r2 [
rnk Ck ]P
(2.9)
=0
k=1
51
The core orbitals replaced by Hay and Wadt large core and Durand-Barthelat ECPs are as
follows:
Li-Ne
= [He]
Na-Ar
= [Ne]
first series
= [Ar]
second series
= [Kr]
third series
= [Xe]4f 14.
The core orbitals replaced by Hay and Wadt small core ECPs are as follows:
K-Cu
= [Ne]
Rb-Ag
= [Ar] 3d10
Cs-Au
= [Kr] 4d10 .
The program evaluates only those integrals for which the overlap between the charge distri-
bution 0 g (page 139) and the most diffuse Gaussian defining the pseudopotential is larger
µ
than a given threshold Tps (the default value is 10-5). See also TOLPSEUD (Section 1.3).
Pseudopotential libraries
The following periodic tables show the effective core pseudo-potentials included as internal
data in the CRYSTAL code. The f terms have not yet been implemented. The program stops
when an ECP including f terms ( =3, M=3) is entered.
HAY AND WADT LARGE CORE ECP. CRYSTAL92 DATA
-------
------------------
Na Mg
Al Si P
S
Cl Ar
------------------------------------------------------
K
Ca Sc Ti V
Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
------------------------------------------------------
Rb Sr Y
Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I
Xe
------------------------------------------------------
Cs Ba
Hf Ta W
Re Os Ir Pt Au Hg Tl Pb Bi
------------------------------------------------------
HAY AND WADT SMALL CORE ECP. CRYSTAL92 DATA
-------------------------------------------------------
K
Ca Sc Ti V
Cr Mn Fe Co Ni Cu
-------------------------------------------------------
Rb Sr Y
Zr Nb Mo Tc Ru Rh Pd Ag
-------------------------------------------------------
Cs Ba
Hf Ta W
Re Os Ir Pt Au
-------------------------------------------------------
DURAND AND BARTHELAT'S LARGE CORE ECP - CRYSTAL92 DATA
------
------------------
Li Be
B
C
N
O
F
Ne
------
------------------
Na Mg
Al Si P
S
Cl Ar
-------------------------------------------------------
K
Ca Sc Ti V
Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
-------------------------------------------------------
Rb
Y
Ag
In Sn Sb
I
-------------------------------------------------------
Tl Pb Bi
-------------------------------------------------------
52
BARTHE, HAYWSC and HAYWLC pseudopotential coefficients and exponents are in-
serted as data in the CRYSTAL code. The data defining the pseudo-potentials where in-
cluded in CRYSTAL92, and never modified. The keyword INPUT allows entering updated
pseudo-potentials, when available. An a posteriori check has been possible for HAYWLC
and HAYWSC only, as the total energy of the atoms for the suggested configuration and
basis set has been published [36, 43]. Agreement with published atomic energies data is satis-
factory (checked from Na to Ba) for Hay and Wadt small core and large core pseudo-potentials,
when using the suggested basis sets. The largest difference is of the order of 10-3 hartree.
For Durand and Barthelat the atomic energies are not published, therefore no check has been
performed. The printed data should be carefully compared with those in the original papers.
The authors of the ECP should be contacted in doubtful cases.
Valence Basis set and pseudopotentials
Hay and Wadt ([36, 43]) have published basis sets suitable for use with their small and large core
pseudopotentials. and in those basis set the s and p gaussian functions with the same quantum
number have different exponent. It is common in CRYSTAL to use sp shells, where basis
functions of s and p symmetry share the same set of Gaussian exponents, with a consequent
considerable decrease in CPU time. The computational advantage of pseudopotentials over
all-electron sets may thus be considerably reduced.
Basis set equivalent to those suggested by Hay and Wadt can be optimized by using CRYSTAL
as an atomic package (page 57), or any