scispace - formally typeset
Search or ask a question
Journal ArticleDOI

CRYSTAL14: A program for the ab initio investigation of crystalline solids

TL;DR: Crystal14 as discussed by the authors is an ab initio code that uses a Gaussian-type basis set: both pseudopotential and all-electron strategies are permitted; the latter is not much more expensive than the former up to the first second transition metal rows of the periodic table.
Abstract: The capabilities of the Crystal14 program are presented, and the improvements made with respect to the previous Crystal09 version discussed. Crystal14 is an ab initio code that uses a Gaussian-type basis set: both pseudopotential and all-electron strategies are permitted; the latter is not much more expensive than the former up to the first-second transition metal rows of the periodic table. A variety of density functionals is available, including as an extreme case Hartree–Fock; hybrids of various nature (global, range-separated, double) can be used. In particular, a very efficient implementation of global hybrids, such as popular B3LYP and PBE0 prescriptions, allows for such calculations to be performed at relatively low computational cost. The program can treat on the same grounds zero-dimensional (molecules), one-dimensional (polymers), two-dimensional (slabs), as well as three-dimensional (3D; crystals) systems. No spurious 3D periodicity is required for low-dimensional systems as happens when plane-waves are used as a basis set. Symmetry is fully exploited at all steps of the calculation; this permits, for example, to investigate nanotubes of increasing radius at a nearly constant cost (better than linear scaling!) or to perform self-consistent-field (SCF) calculations on fullerenes as large as (10,10), with 6000 atoms, 84,000 atomic orbitals, and 20 SCF cycles, on a single core in one day. Three versions of the code exist, serial, parallel, and massive-parallel. In the second one, the most relevant matrices are duplicated, whereas in the third one the matrices in reciprocal space are distributed for diagonalization. All the relevant vectors are now dynamically allocated and deallocated after use, making Crystal14 much more agile than the previous version, in which they were statically allocated. The program now fits more easily in low-memory machines (as many supercomputers nowadays are). Crystal14 can be used on parallel machines up to a high number of cores (benchmarks up to 10,240 cores are documented) with good scalability, the main limitation remaining the diagonalization step. Many tensorial properties can be evaluated in a fully automated way by using a single input keyword: elastic, piezoelectric, photoelastic, dielectric, as well as first and second hyperpolarizabilies, electric field gradients, Born tensors and so forth. Many tools permit a complete analysis of the vibrational properties of crystalline compounds. The infrared and Raman intensities are now computed analytically and related spectra can be generated. Isotopic shifts are easily evaluated, frequencies of only a fragment of a large system computed and nuclear contribution to the dielectric tensor determined. New algorithms have been devised for the investigation of solid solutions and disordered systems. The topological analysis of the electron charge density, according to the Quantum Theory of Atoms in Molecules, is now incorporated in the code via the integrated merge of the Topond package. Electron correlation can be evaluated at the Moller–Plesset second-order level (namely MP2) and a set of double-hybrids are presently available via the integrated merge with the Cryscor program. © 2014 Wiley Periodicals, Inc.

Summary (3 min read)

1 Introduction

  • Derivative-free optimization has experienced a renewed interest over the past decade that has encouraged a new wave of theory and algorithms.
  • The authors explore benchmarking procedures for derivative-free optimization algorithms when there is a limited computational budget.
  • ∗Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 60439.
  • Users with expensive function evaluations are often interested in a convergence test that measures the decrease in function value.
  • These performance profiles are useful to users who need to choose a solver that provides a given reduction in function value within a limit of µf function evaluations.

2 Benchmarking Derivative-Free Optimization Solvers

  • Performance profiles, introduced by Dolan and Moré [5], have proved to be an important tool for benchmarking optimization solvers.
  • Dolan and Moré define a benchmark in terms of a set P of benchmark problems, a set S of optimization solvers, and a convergence test T .
  • Once these components of a benchmark are defined, performance profiles can be used to compare the performance of the solvers.
  • In this section the authors first propose a convergence test for derivative-free optimization solvers and then examine the relevance of performance profiles for optimization problems with expensive function evaluations.

2.2 Data Profiles

  • The authors can use performance profiles with the convergence test (2.2) to benchmark optimization solvers for problems with expensive function evaluations.
  • Performance profiles compare different solvers, while data profiles display the raw data.
  • The authors illustrate the differences between performance and data profiles with a synthetic case involving two solvers.
  • Assume that solver S1 requires k1 simplex gradients to solve each of the first n1 problems, but fails to solve the remaining n2 problems.
  • The limiting value of ρs(α) as α→∞ is the percentage of problems that can be solved with µf function evaluations.

3 Derivative-Free Optimization Solvers

  • The selection of solvers S that the authors use to illustrate the benchmarking process was guided by a desire to examine the performance of a representative subset of derivative-free solvers, and thus they included both direct search and model-based algorithms.
  • The authors note that some solvers were not tested because they require additional parameters outside the scope of this investigation, such as the requirement of bounds by imfil [8, 15].
  • Other implementations of the Nelder-Mead method exist, but this code performs well and has a reasonable default for the size of the initial simplex.
  • The NEWUOA code requires an initial starting point x0, a limit on the number of function evaluations, and the initial trust region radius.
  • The authors effectively set all termination parameters to zero so that all codes terminate only when the limit on the number of function evaluations is exceeded.

4 Benchmark Problems

  • The benchmark problems the authors have selected highlight some of the properties of derivativefree solvers as they face different classes of optimization problems.
  • The problems in the benchmark set P are defined by a vector (kp, np,mp, sp) of integers.
  • These functions are twice continuously differentiable on the level set associated with x0.
  • The authors believe the noise in this type of simulation is better modeled by a function with both high-frequency and low-frequency oscillations.
  • An advantage of the benchmark problems P is that a set of piecewise-smooth problems PPS can be easily derived by setting f(x) = m∑ k=1 |fk(x)|. (4.6) These problems are continuous, but the gradient does not exist when fk(x) = 0 and ∇fk(x) 6= 0 for some index k.

5 Computational Experiments

  • The authors now present the results of computational experiments with the performance measures introduced in Section 2.
  • This data is then processed to obtain a history vector hs ∈.
  • For each problem, p ∈ P, fL was taken to be the best function value achieved by any solver within µf function evaluations, fL = mins∈S hs(xµf ).
  • The authors present the data profiles for τ = 10−k with k ∈ {1, 3, 5, 7} because they are interested in the short-term behavior of the algorithms as the accuracy level changes.
  • This accuracy level is mild compared to classical convergence tests based on the gradient.

5.1 Smooth Problems

  • The data profiles in Figure 5.1 show that NEWUOA solves the largest percentage of problems for all sizes of the computational budget and levels of accuracy τ .
  • With a budget of 10 simplex gradients and τ = 10−5, NEWUOA solves almost 35% of the problems, while both NMSMAX and APPSPACK solve roughly 10% of the problems.
  • A benefit of the data profiles is that they can be useful for allocating a computational budget.
  • Performance differences are also of interest in this case.
  • Both plots in Figure 5.2 show that the performance difference between solvers decreases as the performance ratio increases.

5.2 Noisy Problems

  • The authors now present the computational results for the noisy problems PN as defined in Section 4.
  • An interesting difference between the data profiles for the smooth and noisy problems is that solver performances for large computational budgets tend to be closer than in the smooth case.
  • For τ = 10−5, NEWUOA is the fastest solver on about 60% of the noisy problems, while it was the fastest solver on about 70% of the smooth problems.
  • The performance differences between the solvers are about the same.

5.3 Piecewise-Smooth Problems

  • The computational experiments for the piecewise-smooth problems PPS measure how the solvers perform in the presence of non-differentiable kinks.
  • The data profiles for the piecewise-smooth problems, shown in Figure 5.5, show that these problems are more difficult to solve than the noisy problems PN and the smooth problems PS .
  • Another interesting observation on the data profiles is that APPSPACK solves more problems than NMSMAX with τ = 10−5 for all sizes of the computational budget.
  • The same behavior can be seen in the performance profile with τ = 10−1, but now the initial difference in performance is larger, more than 40%.

6 Concluding Remarks

  • The authors interest in derivative-free methods is motivated in large part by the computationally expensive optimization problems that arise in DOE’s SciDAC initiative.
  • This convergence test relies only on the function values obtained by the solver and caters to users with an interest in the short-term behavior of the solver.
  • Performance of derivative-free solvers for larger problems is of interest, but this would require a different set of benchmark problems.
  • The authors computational experiments used default input and algorithmic parameters, but the authors are aware that performance can change for other choices.

Acknowledgments

  • The authors are grateful to Christine Shoemaker for useful conversations throughout the course of this work and to Josh Griffin and Lúıs Vicente for valuable comments on an earlier draft.
  • Lastly, the authors are grateful for the developers of the freely available solvers they tested for providing us with source code and installation assistance: Genetha Gray and Tammy Kolda (APPSPACK[10]), Nick Higham (MDSMAX and NMSMAX [13]), Tim Kelley (nelder and imfil [15]), Michael Powell (NEWUOA [20] and UOBYQA [19]), and Ana Custódio and Lúıs Vicente (SID-PSM [4]).

Did you find this useful? Give us your feedback

Figures (29)

Content maybe subject to copyright    Report

CRYSTAL14: A Program for the Ab initio Investigation of Crystalline Solids
Roberto Dovesi,
1,
Roberto Orlando,
1
Alessandro Erba,
1
Claudio M. Zicovich-Wilson,
2
Bartolomeo Civalleri,
1
Silvia Casassa,
1
Lorenzo Maschio,
1
Matteo Ferrab one,
1
Marco De la
Pierre,
1
Philippe D’Arco,
3, 4
Yves No¨el,
3, 4
Mauro Caus´a,
5
Michel erat,
6
and B ernard Kirtman
7
1
Dipartimento di Chimica, and Centre of Excellence NIS (Nanostructured I nterfaces and Surfaces),
Universit`a di Torino, via Giuria 5, I-10125 Torino (Italy)
2
Facultad de Ciencias, Universidad Aut´onoma del Estado de Morelos,
Av. Universidad, 1001, Col. Chamilpa, 62209 Cuernavaca, Morelos (Mexico)
3
Sorbonne Universit´es, UPMC Uni v Paris 06, UMR 7193,
Institut des Sciences de la Terre Paris (iSTeP), F-75005 Paris, (France)
4
CNRS, UMR 7193, Institut des Sciences de la Terre Paris (iSTeP), F-75005 Paris, (France)
5
Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale,
Universita’ di Napoli Federico II , Napoli (Italy)
6
Equipe de Chimie Physique, IPREM UMR5254,
Universit´e de Pau et des Pays de l’Adour, 64000 Pau (France)
7
Department of Chemistry and Biochemistry, University of California, Santa Barbara, California 93106 (USA)
(Dated: February 12, 2014)
The capabilities of the Crystal14 program are presented, and the improvements made with respect
to the previous Crystal09 version discussed. Crystal14 is an ab initio code that uses a Gaussian-
type basis set: both pseudopotential and all-electron strategies are permitted; the latter is not much
more expen sive than the former up to the first-second transition metal rows of th e period ic table.
A variety of density functionals is available, including as an extreme case Hartree-Fock; hybrids
of various nature (global, range-separated, double) can be used. In particular, a very efficient
implementation of global hybrids, such as popular B3LYP and PBE0 prescriptions, allows for such
calculations to be performed at relatively low computational cost. The program can treat on th e
same grounds 0D (molecules), 1D (polymers), 2D (slabs), as well as 3D (crystals) systems. No
spurious three-dimensional periodicity is required for low-dimensional systems as happens when
plane-waves are used as a basis set. Symmetry is fully exploited at all steps of the calculation;
this permits, for example, to investigate nanotubes of incr easing radius at a nearly constant cost
(better than linear scaling!) or to perform Self-Consistent-Field (SCF) calculations on fullerenes
as large as (10,10), with 6000 atoms, 84000 atomic orbitals and 20 SCF cycles, on a single core
in one day. Three versions of the code exist, serial, parallel and massive-parallel. In the second
one the most relevant matrices are duplicated whereas in th e third one the matrices in reciprocal
space are distributed for diagonalization. All the relevant vectors are now dynamically allocated
and deallocated after use, making Crystal14 much more agile than the previous version, in which
they were statically allocated. The program now fits more easily in low-memory machines (as many
supercomputers nowadays are). Crystal14 can be used on parallel machines up to a high number
of cores (benchmarks up to 10240 cores are documented) with good scalability, the main limitation
remaining the diagonalization step. Many tensorial p roperties can be evaluated in a fully automated
way by using a single input keyword: elastic, piezoelectric, photoelastic, dielectric, as well as first
and second hyperpolarizabilies, electric field gradients, Born tensors, etc. Many tools permit a
complete analysis of the vibrational properties of crystalline compoun ds. The infrared and Raman
intensities are now computed analytically and related spectra can be generated. I sotopic shifts are
easily evaluated, frequencies of only a fragment of a large system computed and nuclear contribution
to the dielectric tensor determined. New algorithms have been devised for the investigation of solid
solutions and disordered systems. The topological analysis of the electron charge d ensity, according
to the Q uantum Theory of Atoms in Molecules, is now incorporated in the code via the integrated
merge of the Topond package. Electron correlation can be evaluated at the oller-Plesset second
order level (namely MP2) and a set of double-hybrids are presently available via the integrated
merge with the Cryscor p rogram.
I. INTRODUCTION
Quantum mechanical ab initio simulation is rapidly
gaining a more important ole in many scientific c ommu-
nities due to decr e asing computational c ost, as well as the
availability of computer programs of increasing capability
and ease of use. In particular, the number of computer
codes devoted to periodic systems has been rapidly grow-
ing. Crystal88 was the first such code to be distributed
publicly 25 years ago through the Quantum Chemistry
Program Exchange.
1,2
Since then six other releases have
followed in 1992, 1995, 1998, 2003, 2006, 2009 a nd now
there is the curre nt one, Crystal14. This computer
program can be used to study the properties of many
types of compounds characterized by pe riodicity in one
dimension (quasilinear and helical polymers, nanotubes),
Dovesi, R. and Orlando, R. and Erba, A. and Zicovich-Wilson, C. and Civalleri, B. and Casassa, S. and Maschio, L. et al. 2014.
CRYSTAL14: A program for the ab initio investigation of crystalline solids. International Journal of Quantum Chemistry.
114 (19): pp. 1287-1317.

2
two dimensions (monolayers, slabs), or three dimensions
(crystals, solid solutions, substitutionally disordered sys-
tems). As a limiting case, molecules can also be studied.
Despite the many improvements and generalizations
that have been introduced s ince its first release, the
basic aspe cts of Crystal have remained the same.
Thus, this program computes the electronic structure
of periodic s ystems within the Hartree-Fock (originally)
and DFT single particle models using Bloch functions.
A special feature of the code is that the crystal or-
bitals are expanded as linear combinations of atom-
centered Gaussian-type functions. Powerful scr e e ning
techniques are employed to exploit real space local-
ity, which is another distinguishing characteristic of
Crystal. Restricted (closed shell) and unres tricted
(spin-polarized) calculations can be performed with all-
electron or with valence -only basis sets using effective
core pseudo- potentials.
Another unique feature is the extensive exploitation of
symmetry to achieve computational efficiency. Besides
the 230 space groups, 80 two-sided plane groups, 99 rod
groups, and 32 crystallographic point gr oups, there is
provision for molecular point group symmetry (e.g. icosa-
hedral) as well as helical symmetry. Automatic tools al-
low users to obtain lower dimensionality systems from
3D structures by specification of a few geometrical pa-
rameters. Slabs (2D periodic), nanorods (1D periodic)
and nanoc rystals (0D) ar e easily generated from 3D crys -
talline structures; nanotubes (1D) and fullerenes (0D)
can be co ns tructed from 2D sheets o r multi-layered slabs
(Section II A). Full use of symmetry involves all steps of a
Crystal calculation, leading to dra stically reduced com-
putation time and required memory, as well as improved
task farming in parallel ca lculations (Section II).
3–7
Sym-
metry is also used to select the independent elements of
tensor properties for computation.
Several algorithms of Crystal14 now rely on the nu-
merical accur acy of the geometry optimizer: search for
equilibrium structures and transition states, volume- or
pressure-constrained minimizations for the determina-
tion of the equation of state of bulk crystals, nuclea r
relaxation of strained lattices for the computation of elas-
tic, piezoelectric, photoelastic tensors, etc. Recent devel-
opments of the algorithms of the geometry optimizer are
illustrated in Section III.
A wide variety of crystal properties can now be com-
puted automatically. They include third and fourth rank,
as well as first and second rank tensors. Amongst the
former are the fourth rank elastic tensor along with the
related seismic wave velocities;
8–11
the third rank di-
rect and co nverse piezoelectric tensors,
12
and the fourth
rank photoelastic Pockels’ tensor (Section IV).
13
In ad-
dition to the dielectric (or polarizability) tensor (static
and frequency-dependent), the second- a nd third-order
electric sus c eptibilities (or third- and fourth-order hyper-
polarizabilities) are computed analytically via the Cou-
pled Perturbed Hartree-Fock/Kohn-Sham (CPHF/KS)
method (Section V).
14–19
TABLE I: Physical properties that can be compu ted with Crystal14. For each property, its formula and tensor rank are given
along with the general method of computation, which may be either analytical (A) or involve a mixed analytical/numerical
(A/N) scheme. Here, t, u, v, w = x, y, z represent Cartesian indices.
Property Tensor rank Formu la Definitions Method
atomic gradient 1 g
a
t
=
E
r
a,t
Energy E and atomic position vector r
a
, of atom a; A
cell gradient 2 G
tu
=
E
A
tu
Energy E and lattice metric matrix A = (a
1
, a
2
, a
3
); A
polarizability 2 α
tu
=
2
E
ǫ
t
ǫ
u
Energy E and electric field ǫ; A
Born charge 2 Z
a
tu
=
2
E
ǫ
t
r
a,u
Energy E, electric field ǫ and atomic position vector r
a
; A
electric field gradient 2 q
tu
=
ǫ
t
r
u
Electric field ǫ and position vector r; A
Hessian 2 H
ab
tu
=
2
E
r
a,t
r
b,u
Energy E and atomic position vector r
a
; A/N
direct piezoelectricity 3 e
tuv
=
P
t
η
uv
Polarization P and rank-2 strain tensor η; A/N
converse piezoelectricity 3 d
tuv
=
η
uv
ǫ
t
Electric field ǫ and rank-2 strain ten sor η; A/N
first hyperpolarizability 3 β
tuv
=
3
E
ǫ
t
ǫ
u
ǫ
v
Energy E and electric field ǫ; A
Raman intensity 3 I
a
tuv
=
3
E
ǫ
t
ǫ
u
r
a,v
Energy E, electric field ǫ and atomic position vector r
a
; A
elasticity 4 C
tuvw
=
2
E
η
tu
η
vw
Energy E and rank-2 strain tensor η; A/N
photoelasticity 4 p
tuvw
=
ǫ
1
tu
η
vw
Rank-2 dielectric tensor ǫ and rank-2 strain tensor η; A/N
second hyperpolarizability 4 γ
tuvw
=
4
E
ǫ
t
ǫ
u
ǫ
v
ǫ
w
Energy E and electric field ǫ; A

3
Previous versions of the program, starting with
Crystal03,
20,21
have included efficient computation
of v ibration frequencies (Hessian matrix) and related
properties. Now we have added analytical calculation
of infrared (IR) and Raman intensities through the
CPHF/KS scheme
22–25
and automated computation of
IR and Raman spectra.
26–28
These vibrationa l properties
are listed in Table I, together with those discussed in the
previous paragraph, as the main tensor properties avail-
able in Crystal14. In addition, the program now con-
tains improved algorithms for calculating phonon disper-
sion and anisotropic displacement parameters (ADP)
29,30
(see Section VI).
New algorithms have rec e ntly been developed for the
study of solid solutions and, mor e generally, disordered
systems.
31,32
As far as solid state solutions are concerned,
for any substitution fraction x within a given series , the
program finds the total number of atomic configur ations
and determines the symmetry-irreducible configurations
among them. Symmetry irreducible configurations can,
then, be explored (i.e. optimized) automatically. Section
VII is devoted to the presentation of these new tools.
We have also developed new tools for the analysis of
electron charge and momentum densities. Section VIII
deals with i) topological analysis of the electron charge
density, performed a c c ording to Bader’s pr e scriptions us -
ing Gatti’s Topond package,
33
that is now integrated
with the Crystal14 program; ii) computation of ADPs
and Debye-Waller thermal damping for X-ray structure
factors,
29
and iii) new algorithms for analyzing the elec-
tron momentum density.
34,35
For DFT calculations, the loca l-density and
generalized-gradient approximations (LDA and GGA)
already available in previous versions of the program
have been augmented w ith exchange-correlation func-
tionals corresponding to the third, fourth and fifth rung
of “Jacob’s ladder” . The semilocal and hybrid meta-
GGA functionals, range-separa ted hybrids and double
hybrids that have been implemented are discussed in
Section IX.
Crystal may be run either sequentially on a single
processor or in parallel. Parallel process ing, in turn,
can be done either through a r e plicated data pro cedure
(Pcrystal), wherein all the most relevant quantities are
copied on each node, or according to a MPP (massively
parallel) strategy, in which the la rge matrices are pa r-
titioned and distributed amongst the cor e s. The MP-
Pcrystal version of the program, first r e le ased in 2010,
is advantageous for systems with a large unit c e ll and low
symmetry. Since then, we have improved performance so
that c alculations now scale linearly up to thousands of
cores. These recent advance s are described in Section X.
The paper is organized as follows: each section corre-
sp onds to a particular capability of the Crystal14 pr o-
gram; the newly developed features are illustrated from
a general point of view and a few examples are given of
their application to systems of interest. Technical details
about how to run specific calculations and extract the
corresponding information (input/output structure) can
be found in the CRYSTAL14 User’s Manual
36
and in the
many tutorials available a t the Crystal website.
37
II. USE OF SYMMETRY IN CRYSTAL
Symmetry plays a crucial rˆole in the study of crys-
talline compounds and its use in simulation can be enor-
mously advantageous because: i) it greatly simplifies the
set of data to be given as input; ii) it improves perfor-
mance dramatically; iii) the amo unt of memory storage
is reduced significantly; and iv) it aids in comparing with
exp eriment.
The data needed to define a crystalline structure is
greatly reduced if its space group is known because the
unit cell can be generated automatically from knowl-
edge of the asymmetric unit only. The s ame rule ap-
plies to slabs, rods, nanotubes and molecules (including
fullerenes). Two s ignificant examples are nanotubes and
fullerenes, as obtained by geometrical construction from
graphene (see Section II A for further details and exam-
ples). The same kind of simplification applies to many
structure manipulations that modify either translational
or point symmetry, or bo th. For instance, a 2D graphene
monolayer can be automatically cut out of a 3D graphite
crystal simply by using the keyword SLABCUT (see below).
A. Low Dimensionality Systems
Twenty or thirty years ago, theoretical studies could
be classified into three main categories: standard molec-
ular calculations, fully periodic (bulk) calculations, and
cluster calculations. The first and last kind of calcula-
tions were performed with molecular” codes, the second
kind with “solid state” codes, the two using completely
different tools and technologies: plane waves vs atomic
orbitals, DFT vs HF, all-electro n vs pseudopotential ba-
sis set. Crystal, on the contrary, maintains full com-
patibility along the 3D 2D 1D 0D periodicity
series. At the same time, various tools have been imple-
mented for an automated, consistent, easy construction
of low-dimensional systems, such as:
Slabs (3D 2D) - This option allows for cutting a 2D
periodic sla b out of a 3D bulk structure, like a graphene
monolayer out of the 3D graphite crystal simply by spec-
ifying: i) the crystallographic Miller’s indices of the cut-
ting plane (0 0 1 in this case); ii) the label of the layer
corresponding to the bottom cut (any layer in graphite);
iii) the number of layers forming the s lab (1 in the case
of graphene). The origin of the 2D unit cell is s e t auto-
matically so as to maximize symmetry. In principle, this
option allows for any termination of the surface, any stoi-
chiometry and sla b polarity beca use intrinsically polar”
surfaces exist, where the electrical dipole perpendicular
to the surface can only be re moved by altering the surface
structure, for instance, by br e aking stoichiometry.

4
FIG. 1: Structure of (a) a single-walled carbon nanotube and
(b) the (10,10) icosahedral carbon fullerene, made up of 6000
atoms.
Nanotubes (2D 1D) - Nanotubes (see Figure 1 a)
of any s iz e and order of symmetry can be built automat-
ically by rolling up a 2D slab.
38
That is done by specify-
ing a pair of integers to define a roll-up vector in terms
of the slab unit vector comp onents. The rolling vector
is perpendicular to the nanotube axis and its modulus is
the nanotube circ umference . In order to build the (120,0)
carbon nanotube from a graphene sheet, for instance, the
value of those two integers is: 120 and 0.
Fullerenes (2 D 0D) - A fullerene cage can be
obtained starting from any hexagonal sheet;
39
available
shapes are: icosahedron, octahedron and tetrahedron (all
consisting of e quilateral triangular fac e s). The required
information includes a pair of integers representing the
components of 2D unit cell vectors defining the side of the
triangular face, the fullerene point group and polyhedron
type. For example, a (10,10) icosahedral fullerene (see
Figure 1 b) is obtained from a graphene sheet by speci-
fying the following data: 10, 10, IH, ICOSA (standing for
icosahedro n).
Nanorods (3D 1D) - To build crystals in 1 D”
starting from the corr e sponding 3D bulk s tructure.
Nanorods can be used as models for crystal edges (rele-
vant to catalysis), and real nano-objects. The required
input steps to build a nanorod are similar to those used to
create slabs: information about the crystal plane, termi-
nation a nd thickness of the two crystal planes defining it.
At this stage, a nanorod can exhibit sharp and unphys-
ical edges that can be smoothed out by further cutting
the rod along additional crystal planes parallel to the pe-
riodic direction. The origin of the nanorod cell is shifted
automatically to maximize symmetry of the rod group.
See Figure 2 (a) for an example.
Nanocrystals (3D 0D) - Non-periodic nanosys-
tems are defined as nanocr ystals, provided they pr eserve
a crystalline structure. Very small nanosystems loos-
ing the parent crystal structure and s toichiometry are
commonly referred to as nanoclusters. A new feature in
Crystal14 allows users to create nanocrystals from bulk
in a controlled way by generating a super c e ll with faces
parallel (2 by 2) to a set of three crystal planes (see Fig-
ure 2 (b) for an example). The same set of information
FIG. 2: (a) MgCl
2
nanorod built on the (104) and (001)
surfaces;
40
(b) the Mg
216
O
216
nanocrystal defined by the
(100), (010) and (001) surfaces of bulk MgO, as further edited
by the (111) family of planes; (c) Wulff polyhedron of CeO
2
:
only (111) and (331) surfaces are represented, with contribu-
tions to the total area of 86% and 14%, respectively. Surface
energies are reported in Ref. 41.
used to create slabs (and nanorods ) needs to be replicated
three times in this case. As for nanorods, a nanocr ystal
can be smoo thed by cutting e dges along other crystal
planes, so breaking the original chemical stoichiometry.
For this re ason, chemical composition is calculated and
printed at each step of the procedure. The or igin of the
nanocrystal is shifted automatically to maximize point
symmetry.
Wulff polyhedron construction (2D 0D) - After
a sy stematic study of a se lec ted set of c rystal surfaces,
where a set of well converged sur face formation energies is
obtained, the thermodynamic equilibrium crystal sha pe
can be calculated using the Gibbs-Wulff theory to build
a Wulff polyhedron.
B. Symmetry and efficiency
Use of trans lation invariance is manda tory to trans-
form a cluster calculation for a piece of a bulk sys tem into
a periodic bulk calculation. Bloch Functions (BFs) are
employed, in that regard, to transform an infinite Hamil-
tonian matrix into block-diagonal fo rm. Each block cor-
responds to a k point in the first Brillouin Z one (BZ) (see
the first two panels in Figure 3). All periodic codes are
based on the use of BFs.
On the other hand, use of point symmetry is not e ssen-
tial to make a periodic calculation affordable and, indeed,
most periodic co des neglect it. However, point symmetry
can be very advantageous in per iodic calculations (much
more than in mo lecular codes). It is exhibited by most
crystalline systems and can be invoked at several steps
of an ab initio calculation, thereby drastically reducing
computation time and resources. The key steps where
savings occurs are given in the following for the case
where the BFs a re co nstructed from a local basis:
Diagonalization of the Fock matrix is restricted to
the s ubset of k points in the irreducible BZ. The
eigenvalues in a star of k points (a set of points

5
FIG. 3: Block-factorization of the Fock matrix for periodic
systems. F
g
: in the basis of AOs, non-packed form (borders
are blurry to indicate that such matrix is infinite in princi-
ple); F
k
: in the basis of BFs orbitals; F
k
: in the basis of
SACOs. Also in the last two cases the matrix is infite, but
block diagonal
that are symmetry related) are the same and the
eigenvectors can be generated by symmetry opera-
tions. This kind of symmetry is used in Crystal
since its first re le ase (Crystal88).
Time required for the calculation of one- and two-
electron integrals is reduced by a factor of up to
the number of point symmetry operators in the
group. Aga in, such use of symmetry dates back
to Crystal88.
3
Diagonalization of the Fock matrix can be speeded
up dramatically in the basis of the Symmetry
Adapted Crystalline Orbitals (SACO)/Symmetry
Adapted Molecular O rbitals (SAMO), shown in the
last panel of Figure 3). SACO/SAMO are gener-
ated automatically in Crystal from the selected
basis set of Atomic Orbitals (AO) in the unit cell,
with no need for a dditional information about irre-
ducible representations (irrep) or characters. This
part of the code was implemented about 15 years
ago.
4,5
The savings factor in computation time is
roughly proportional to the third power of the ratio
between the number of AOs in the basis set (N
AO
)
and the size of the largest block in the Fock matrix ,
when repre sented in the SACO/SAMO basis.
Construction of the density matrix scales with the
third power of the basis set size (N
3
AO
), as each of
the N
2
AO
matrix elements is obtained by summing
over all occupied crystalline orbitals (very roughly
N
AO
). The advantage is obtained by building the
matrix in the SACO/SAMO basis first.
Reduction of computing time thanks to the use of sym-
metry is not the only issue when handling very large
unit cell cases. Memory requirements can also become
a bottleneck, if not properly managed at every step of
a calculation. Storage of the Fock, overlap and density
matrices as full square matrices in the AO-Bloch func-
tion basis represent the main bottleneck and need to be
avoided. Since both one- and two-electron integrals are
evaluated in the AO basis, a set of back and forth trans-
formations (SACO to AO and AO to SACO)is required.
However, it is possible to switch directly from the AO
to the SACO/SAMO basis and use the latter through-
out. T his has been implemented in Crystal14.
6,7
with
a drastic r eduction of b oth running time and memory al-
location. The larger the point group size, the bigger the
reduction of computational resources (see Section X for
details).
In the following we document the savings from full use
of symmetry in the case of Carbon Nanotubes (CNT), in
particular for the (n,0) family. The number of symmetry
operators in CNTs increases with the tube size, i.e. with
n. Thus, in principle, there is no limit to the number
of non-purely translational e lements. Neglecting mirror
planes parallel and perpendicular to the tube axis, (n,0)
nanotubes pos sess 2n roto-trans lation symmetry opera-
tors. The largest nanotub e considered here, with n =
320, has 640 such operators and 1280 atoms in its unit
cell. This is a huge number compared to a maximum of 48
point operators for standard cubic crystalline s ystems or
120 for fullerenes, in the molecular context. Thus, large
nanotubes are ex pected to show the maximum symmetry
savings factor in terms of computation time and memory
allocation. They also represent a severe test for the flexi-
bility of our code and its algorithms, as will be discuss e d
below. A more detailed description of the effect of sym-
metry c an be found in Refs. 6 and 7. Large nanotubes
are, clearly, an extreme case. Evidence of the efficiency
of the code for systems with lower point symmetry (from
1 to 48 operators) is given in one of the next sections.
Tables II and I II show the effect of full use of symmetry
on peak memory usa ge and on running time. The various
steps of the calculation can be divided into three units
(corresponding to the three blocks in Tables II and III):
a) Preliminary calculations for symmetry analysis (in-
cluding co ns truction of the character table and
SACOs), index mapping and screening of integrals
(init). This unit includes also the generation of a
set of orthogonal functions through Cholesky de-
composition.
b) All steps involved in the SCF cycle, which are it-
erated as many times as nece ssary to achieve con-
vergence (15 iterations in the present case). In the
“direct SCF” strategy this unit includes: calcula-
tion of 1- and 2-electron integrals, a s well as in-
tegrals approximated by multipola r expa ns ion of
the electron charge density (pole); trans fo rmation
of the Fock matrix from the AO to the SACO basis
and its diagonalization (Fock) + (diag ); construc-
tion of the density matrix in the AO basis (dens);
and integration of the exchange-correlatio n density
functional (dft).
c) Calculation of the total energy gradients with re-
sp e c t to the nuclear coordinates (gr ad) at the end
of the SCF proc e ss.

Citations
More filters
Journal ArticleDOI
TL;DR: The Crystal program as discussed by the authors adopts atom-centered Gaussian-type functions as a basis set, which makes it possible to perform all-electron as well as pseudopotential calculations.
Abstract: The latest release of the Crystal program for solid-state quantum-mechanical ab initio simulations is presented. The program adopts atom-centered Gaussian-type functions as a basis set, which makes it possible to perform all-electron as well as pseudopotential calculations. Systems of any periodicity can be treated at the same level of accuracy (from 0D molecules, clusters and nanocrystals, to 1D polymers, helices, nanorods, and nanotubes, to 2D monolayers and slab models for surfaces, to actual 3D bulk crystals), without any artificial repetition along nonperiodic directions for 0–2D systems. Density functional theory calculations can be performed with a variety of functionals belonging to several classes: local-density (LDA), generalized-gradient (GGA), meta-GGA, global hybrid, range-separated hybrid, and self-consistent system-specific hybrid. In particular, hybrid functionals can be used at a modest computational cost, comparable to that of pure LDA and GGA formulations, because of the efficient implementation of exact nonlocal Fock exchange. Both translational and point-symmetry features are fully exploited at all steps of the calculation, thus drastically reducing the corresponding computational cost. The various properties computed encompass electronic structure (including magnetic spin-polarized open-shell systems, electron density analysis), geometry (including full or constrained optimization, transition-state search), vibrational properties (frequencies, infrared and Raman intensities, phonon density of states), thermal properties (quasi-harmonic approximation), linear and nonlinear optical properties (static and dynamic [hyper]polarizabilities), strain properties (elasticity, piezoelectricity, photoelasticity), electron transport properties (Boltzmann, transport across nanojunctions), as well as X-ray and inelastic neutron spectra. The program is distributed in serial, parallel, and massively parallel versions. In this paper, the original developments that have been devised and implemented in the last 4 years (since the distribution of the previous public version, Crystal14, occurred in December 2013) are described.

1,108 citations


Cites background or methods from "CRYSTAL14: A program for the ab ini..."

  • ...In Figure 2 we report the effect of the DIIS accelerator as compared with the Fock mixing + level-shifting scheme that was default in CRYSTAL14, benchmarked on our test set....

    [...]

  • ...In particular, the new developments and improvements with respect to the previous major release of the program, CRYSTAL14 (Dovesi et al., 2014), are discussed into detail....

    [...]

  • ...In particular, the new features of the code with respect to the previous major release, namely CRYSTAL14, include: (a) implementation of a DIIS scheme for SCF and CPHF/ KS convergence acceleration (now active by default); (b) fully automated implementation of the QHA for volume-dependent thermal properties; (c) input parameter-free implementation of Grimme’s DFT-D3 correction for weak dispersive interactions as well as a semiempirical composite method for molecular crystals; (d) calculation of elastic constants under pressure, directional elastic wave velocities, the piezo-optic tensor, and the piezoelectric tensor through an analytical CPHF/KS approach; (e) calculation of dynamic linear polarizabilities, SHG, and Pockels effect; (f ) self-consistent system-specific hybrid functionals; (g) improved scalability of the massively parallel version of the program, MPPCRYSTAL; (h) implementation of new tools for magnetic systems (evaluation of spin contamination, restricted open-shell HF, and noninteger spin locking); (i) atomic partition of the electron charge density according to the Hirshfeld-I scheme; (j) calculation of total and projected PDOS as well as its neutron-weighted counterpart (for simulation of INS); (k) calculation of X-ray diffraction spectra; and (l) evaluation of several electronic transport properties (electrical conductivity, Seebeck coefficient, electronic contribution to thermal conductivity, and transport across nanojunctions)....

    [...]

  • ...The CPHF/KS scheme was further extended to second-order in the perturbed wavefunction in CRYSTAL14, thus allowing for the calculation of static nonlinear properties (namely, second hyperpolarizabilities and n + 1 rule first hyperpolarizabilities) (Ferrero, Rérat, Orlando, Dovesi, & Bush, 2008; Orlando, Ferrero, Rérat, Kirtman, & Dovesi, 2009)....

    [...]

  • ...In the following of this section, we summarize the main developments made for these properties in the CRYSTAL17 version of the program with respect to CRYSTAL14....

    [...]

Journal ArticleDOI
TL;DR: The capabilities and design philosophy of the current version of the PySCF package are document, which is as efficient as the best existing C or Fortran‐based quantum chemistry programs.
Abstract: Python-based simulations of chemistry framework (PySCF) is a general-purpose electronic structure platform designed from the ground up to emphasize code simplicity, so as to facilitate new method development and enable flexible computational workflows. The package provides a wide range of tools to support simulations of finite-size systems, extended systems with periodic boundary conditions, low-dimensional periodic systems, and custom Hamiltonians, using mean-field and post-mean-field methods with standard Gaussian basis functions. To ensure ease of extensibility, PySCF uses the Python language to implement almost all of its features, while computationally critical paths are implemented with heavily optimized C routines. Using this combined Python/C implementation, the package is as efficient as the best existing C or Fortran-based quantum chemistry programs. In this paper, we document the capabilities and design philosophy of the current version of the PySCF package. WIREs Comput Mol Sci 2018, 8:e1340. doi: 10.1002/wcms.1340 This article is categorized under: Structure and Mechanism > Computational Materials Science Electronic Structure Theory > Ab Initio Electronic Structure Methods Software > Quantum Chemistry

1,042 citations

Journal ArticleDOI
TL;DR: This Review describes the recent developments (including some historical aspects) of dispersion corrections with an emphasis on methods that can be employed routinely with reasonable accuracy in large-scale applications.
Abstract: Mean-field electronic structure methods like Hartree–Fock, semilocal density functional approximations, or semiempirical molecular orbital (MO) theories do not account for long-range electron correlation (London dispersion interaction). Inclusion of these effects is mandatory for realistic calculations on large or condensed chemical systems and for various intramolecular phenomena (thermochemistry). This Review describes the recent developments (including some historical aspects) of dispersion corrections with an emphasis on methods that can be employed routinely with reasonable accuracy in large-scale applications. The most prominent correction schemes are classified into three groups: (i) nonlocal, density-based functionals, (ii) semiclassical C6-based, and (iii) one-electron effective potentials. The properties as well as pros and cons of these methods are critically discussed, and typical examples and benchmarks on molecular complexes and crystals are provided. Although there are some areas for furthe...

932 citations

Journal ArticleDOI
TL;DR: The potential of metal–organic frameworks (MOFs) and covalent organic frameworks (COFs) as platforms for the development of heterogeneous single-site catalysts is reviewed thoroughly.
Abstract: Heterogeneous single-site catalysts consist of isolated, well-defined, active sites that are spatially separated in a given solid and, ideally, structurally identical. In this review, the potential of metal–organic frameworks (MOFs) and covalent organic frameworks (COFs) as platforms for the development of heterogeneous single-site catalysts is reviewed thoroughly. In the first part of this article, synthetic strategies and progress in the implementation of such sites in these two classes of materials are discussed. Because these solids are excellent playgrounds to allow a better understanding of catalytic functions, we highlight the most important recent advances in the modelling and spectroscopic characterization of single-site catalysts based on these materials. Finally, we discuss the potential of MOFs as materials in which several single-site catalytic functions can be combined within one framework along with their potential as powerful enzyme-mimicking materials. The review is wrapped up with our personal vision on future research directions.

785 citations

Journal ArticleDOI
13 Oct 2015-ACS Nano
TL;DR: This work presents an important progress for selective and reversible NO2 sensing by demonstrating an economical sensing platform based on the charge transfer between physisorbed NO2 gas molecules and two-dimensional (2D) tin disulfide (SnS2) flakes at low operating temperatures.
Abstract: Nitrogen dioxide (NO2) is a gas species that plays an important role in certain industrial, farming, and healthcare sectors. However, there are still significant challenges for NO2 sensing at low detection limits, especially in the presence of other interfering gases. The NO2 selectivity of current gas-sensing technologies is significantly traded-off with their sensitivity and reversibility as well as fabrication and operating costs. In this work, we present an important progress for selective and reversible NO2 sensing by demonstrating an economical sensing platform based on the charge transfer between physisorbed NO2 gas molecules and two-dimensional (2D) tin disulfide (SnS2) flakes at low operating temperatures. The device shows high sensitivity and superior selectivity to NO2 at operating temperatures of less than 160 °C, which are well below those of chemisorptive and ion conductive NO2 sensors with much poorer selectivity. At the same time, excellent reversibility of the sensor is demonstrated, which has rarely been observed in other 2D material counterparts. Such impressive features originate from the planar morphology of 2D SnS2 as well as unique physical affinity and favorable electronic band positions of this material that facilitate the NO2 physisorption and charge transfer at parts per billion levels. The 2D SnS2-based sensor provides a real solution for low-cost and selective NO2 gas sensing.

579 citations

References
More filters
Journal ArticleDOI
TL;DR: An efficient scheme for calculating the Kohn-Sham ground state of metallic systems using pseudopotentials and a plane-wave basis set is presented and the application of Pulay's DIIS method to the iterative diagonalization of large matrices will be discussed.
Abstract: We present an efficient scheme for calculating the Kohn-Sham ground state of metallic systems using pseudopotentials and a plane-wave basis set. In the first part the application of Pulay's DIIS method (direct inversion in the iterative subspace) to the iterative diagonalization of large matrices will be discussed. Our approach is stable, reliable, and minimizes the number of order ${\mathit{N}}_{\mathrm{atoms}}^{3}$ operations. In the second part, we will discuss an efficient mixing scheme also based on Pulay's scheme. A special ``metric'' and a special ``preconditioning'' optimized for a plane-wave basis set will be introduced. Scaling of the method will be discussed in detail for non-self-consistent and self-consistent calculations. It will be shown that the number of iterations required to obtain a specific precision is almost independent of the system size. Altogether an order ${\mathit{N}}_{\mathrm{atoms}}^{2}$ scaling is found for systems containing up to 1000 electrons. If we take into account that the number of k points can be decreased linearly with the system size, the overall scaling can approach ${\mathit{N}}_{\mathrm{atoms}}$. We have implemented these algorithms within a powerful package called VASP (Vienna ab initio simulation package). The program and the techniques have been used successfully for a large number of different systems (liquid and amorphous semiconductors, liquid simple and transition metals, metallic and semiconducting surfaces, phonons in simple metals, transition metals, and semiconductors) and turned out to be very reliable. \textcopyright{} 1996 The American Physical Society.

81,985 citations

Journal ArticleDOI
TL;DR: A detailed description and comparison of algorithms for performing ab-initio quantum-mechanical calculations using pseudopotentials and a plane-wave basis set is presented in this article. But this is not a comparison of our algorithm with the one presented in this paper.

47,666 citations

Journal ArticleDOI
TL;DR: In this paper, the Hartree and Hartree-Fock equations are applied to a uniform electron gas, where the exchange and correlation portions of the chemical potential of the gas are used as additional effective potentials.
Abstract: From a theory of Hohenberg and Kohn, approximation methods for treating an inhomogeneous system of interacting electrons are developed. These methods are exact for systems of slowly varying or high density. For the ground state, they lead to self-consistent equations analogous to the Hartree and Hartree-Fock equations, respectively. In these equations the exchange and correlation portions of the chemical potential of a uniform electron gas appear as additional effective potentials. (The exchange portion of our effective potential differs from that due to Slater by a factor of $\frac{2}{3}$.) Electronic systems at finite temperatures and in magnetic fields are also treated by similar methods. An appendix deals with a further correction for systems with short-wavelength density oscillations.

47,477 citations

Journal ArticleDOI
TL;DR: The M06-2X meta-exchange correlation function is proposed in this paper, which is parametrized including both transition metals and nonmetals, and is a high-non-locality functional with double the amount of nonlocal exchange.
Abstract: We present two new hybrid meta exchange- correlation functionals, called M06 and M06-2X. The M06 functional is parametrized including both transition metals and nonmetals, whereas the M06-2X functional is a high-nonlocality functional with double the amount of nonlocal exchange (2X), and it is parametrized only for nonmetals.The functionals, along with the previously published M06-L local functional and the M06-HF full-Hartree–Fock functionals, constitute the M06 suite of complementary functionals. We assess these four functionals by comparing their performance to that of 12 other functionals and Hartree–Fock theory for 403 energetic data in 29 diverse databases, including ten databases for thermochemistry, four databases for kinetics, eight databases for noncovalent interactions, three databases for transition metal bonding, one database for metal atom excitation energies, and three databases for molecular excitation energies. We also illustrate the performance of these 17 methods for three databases containing 40 bond lengths and for databases containing 38 vibrational frequencies and 15 vibrational zero point energies. We recommend the M06-2X functional for applications involving main-group thermochemistry, kinetics, noncovalent interactions, and electronic excitation energies to valence and Rydberg states. We recommend the M06 functional for application in organometallic and inorganometallic chemistry and for noncovalent interactions.

22,326 citations

Journal ArticleDOI
TL;DR: The authors assess various approximate forms for the correlation energy per particle of the spin-polarized homogeneous electron gas that have frequently been used in applications of the local spin density a...
Abstract: We assess various approximate forms for the correlation energy per particle of the spin-polarized homogeneous electron gas that have frequently been used in applications of the local spin density a...

17,531 citations

Frequently Asked Questions (16)
Q1. What have the authors contributed in "Crystal14: a program for the ab initio investigation of crystalline solids" ?

In this paper, Dovesi et al. presented a paper on the NIS ( Nanostructured Interfaces and Surfaces ). 

Due to the fact that core and inner-valence electrons of atoms follow the movement of the respective nuclei, when ECD and related X-ray structure factors (XSF) are considered, it is mandatory to account for the effect of finite temperature, for instance by means of atomic harmonic Debye-Waller thermal factors. 

Convergence of the calculated optical properties to the infinite periodic polymer limit is essentially achieved for a chain length of 50 monomers. 

Geometry optimization is mostly employed in quantum chemical calculations to obtain a nuclear configuration that is either• a stable structure of a given chemical species, or• can be used to estimate the transition configuration along a reaction path leading to determination of the rate for the corresponding process. 

The three acoustic wave velocities, also referred to as seismic velocities, can be labeled as quasi-longitudinal vp, slow quasi-transverse vs1 and fast quasi-transverse vs2, depending on their polarization with respect to the propagation direction. 

Diminishing memory requirements has been a main challenging issue in extending the capabilities of Crystal to handle large unit-cell systems in the HPC context where the general trend implies reduced memory availability both in terms of GB/core and bandwidth. 

The aim of using Eq. (19) as a startingpoint for their treatment here is top avoid explicit gradients of coefficients that would require the solution of an additional set of coupled–perturbed equations. 

12,66Due to its peculiar piezoelectric properties, α-SiO2 is another material that is widely utilized in the electronics industry. 

By imposing equality10between the elastic strain energy of the monolayer and the corresponding vibrational energy of the nanotube, the slope of frequency versus 1/n for the latter can be related to the elastic constants of the monolayer. 

Victor R. Saunders is acknowledged for his invaluable contribution in the development of the code and in making it computationally efficient and numerically stable. 

The recent development of an algorithm for computing CPs from the density matrix of the system, rather than from the crystalline orbitals, made possible the investigation of the effect of the adopted computational method on momentum space properties, even beyond the one-electron approximation, with the MP2 approach implemented in the Cryscor program. 

The rate of increase relative to the corresponding static electronic property is dictated primarily by the number of static fields used to characterize the process. 

Some of the new features of Crystal14 as regards electron densities are: i) the complete topological analysis of the ECD by means of the automated integration of the Topond package into Crystal14 (see Section VIII A); ii) the parallelization, with linear speed-up, of all the algorithms related to ECD and EMD; iii) calculation of anisotropic displacement parameters and DebyeWaller thermal factors for dynamical X-ray structure factors (see Section VIII B); iv) new algorithms for the analysis of the EMD (see Section VIII C). 

As an example, the 239 MBytes for the grad step in the case of the largest nanotube is about 1/10 of the memory commonly available in a single core standard machine. 

For small bandgap polymers, the precision of coupled perturbed (hyper)polarizability calculations is mainly determined by two computational parameters, namely1. 

This is the reason why the memory-occupation curve in the Figure appears as to be tending to an asymptotic value for large values of n, which represents the maximum amount of replicated data being stored to memory during an SCF+G calculation for MCM-41/X1.