Description
Quantum technology is a new field of physics and engineering, which transitions some of the stranger features of quantum mechanics, especially quantum entanglement and most recently quantum tunneling, into practical applications such as quantum computing, quantum cryptography, quantum simulation, quantum metrology, quantum sensing, and quantum imaging.
QUANTUM ESPRESSO: a modular and open-source software project
for quantum simulations of materials
Paolo Giannozzi,
1, 2
Stefano Baroni,
1, 3
Nicola Bonini,
4
Matteo Calandra,
5
Roberto Car,
6
Carlo Cavazzoni,
7, 8
Davide Ceresoli,
4
Guido L. Chiarotti,
9
Matteo Cococcioni,
10
Ismaila Dabo,
11
Andrea Dal Corso,
1, 3
Stefano de Gironcoli,
1, 3
Ralph Gebauer,
1, 12
Uwe Gerstmann,
13
Christos Gougoussis,
5
Anton Kokalj,
14
Michele Lazzeri,
5
Layla Martin Samos Colomer,
1
Nicola
Marzari,
4
Francesco Mauri,
5
Stefano Paolini,
3, 9
Alfredo Pasquarello,
15
Lorenzo Paulatto,
1, 3
Carlo Sbraccia
?
,
1
Sandro
Scandolo,
1, 12
Gabriele Sclauzero,
1, 3
Ari P. Seitsonen,
5
Alexander Smogunov,
12
Paolo Umari,
1
and Renata M. Wentzcovitch
10
1
CNR-INFM Democritos National Simulation Center, 34100 Trieste, Italy
2
Dipartimento di Fisica, Università degli Studi di Udine, via delle Scienze 208, 33100 Udine, Italy
3
SISSA – Scuola Internazionale Superiore di Studi Avanzati, via Beirut 2-4, 34014 Trieste Grignano, Italy
4
Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139 USA
5
Institut de Minéralogie et de Physique des Milieux Condensés,
Université Pierre et Marie Curie, CNRS, IPGP, 140 rue de Lourmel, 75015 Paris, France
6
Department of Chemistry, Princeton University, Princeton NJ 08544, USA
7
CINECA National Supercomputing Center, Casalecchio di Reno, 40033 Bologna, Italy
8
CNR-INFM S3 Research Center, 41100 Modena, Italy
9
SPIN s.r.l. via del Follatoio 12, 34148 Trieste, Italy
10
Department of Chemical Engineering and Materials Science, University of Minnesota,
151 Amundson Hall, 421 Washington Avenue SE, Minneapolis MN 55455, USA
11
Université Paris-Est, CERMICS, Projet Micmac ENPC-INRIA,
6-8 avenue Blaise Pascal, 77455 Marne-la-Vallée Cedex 2, France
12
The Abdus Salam International Centre for Theoretical Physics, Strada Costiera 11, 34014 Trieste Grignano, Italy
13
Theoretische Physik, Universität Paderborn, D-33098 Paderborn, Germany
14
Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia
15
Ecole Polytechnique Fédérale de Lausanne (EPFL), Institute of Theoretical Physics,
and Institut Romand de Recherche Numérique en Physique des Matériaux (IRRMA), CH-1015 Lausanne, Switzerland
QUANTUM ESPRESSO is an integrated suite of computer codes for electronic-structure calculations and
materials modeling, based on density-functional theory, plane waves, and pseudopotentials (norm-conserving,
ultrasoft, and projector-augmented wave). QUANTUM ESPRESSO stands for opEn Source Package for Re-
search in Electronic Structure, Simulation, and Optimization. It is freely available to researchers around the
world under the terms of the GNU General Public License. QUANTUM ESPRESSO builds upon newly-
restructured electronic-structure codes that have been developed and tested by some of the original authors
of novel electronic-structure algorithms and applied in the last twenty years by some of the leading materials
modeling groups worldwide. Innovation and ef?ciency are still its main focus, with special attention paid to
massively-parallel architectures, and a great effort being devoted to user friendliness. QUANTUM ESPRESSO
is evolving towards a distribution of independent and inter-operable codes in the spirit of an open-source project,
where researchers active in the ?eld of electronic-structure calculations are encouraged to participate in the
project by contributing their own codes or by implementing their own ideas into existing codes.
I. INTRODUCTION
The combination of methodological and algorithmic in-
novations and ever-increasing computer power is deliver-
ing a simulation revolution in materials modeling, starting
from the nanoscale up to bulk materials and devices [1].
Electronic-structure simulations based on density-functional
theory (DFT) [2–4] have been instrumental to this revolu-
tion, and their application has now spread outside a restricted
core of researchers in condensed-matter theory and quantum
chemistry, involving a vast community of end users with very
diverse scienti?c backgrounds and research interests. Sus-
taining this revolution and extending its bene?cial effects to
the many ?elds of science and technology that can capital-
?
Present address: Constellation Energy Commodities Group 7th Floor, 61
Aldwich, London, WC2B 4AE, United Kingdom
ize on it represents a multifold challenge. In our view it
is also a most urgent, fascinating and fruitful endeavor, able
to deliver new forms for scienti?c exploration and discovery,
where a very complex infrastructure—made of software rather
than hardware—can be made available to any researcher, and
whose capabilities continue to increase thanks to the method-
ological innovations and computing power scalability alluded
to above.
Over the past few decades, innovation in materials simula-
tion and modeling has resulted from the concerted efforts of
many individuals and groups worldwide, often of small size.
Their success has been made possible by a combination of
competences, ranging from the ability to address meaningful
and challenging problems, to a rigorous insight into theoret-
ical methods, ending with a marked sensibility to matters of
numerical accuracy and algorithmic ef?ciency. The readiness
to implement new algorithms that utilize novel ideas requires
total control over the software being used–for this reason,
the physics community has long relied on in-house computer
2
codes to develop and implement new ideas and algorithms.
Transitioning these development codes to production tools is
nevertheless essential, both to extensively validate new meth-
ods and to speed up their acceptance by the scienti?c com-
munity. At the same time, the dissemination of codes has to
be substantial, to justify the learning efforts of PhD students
and young postdocs who would soon be confronted with the
necessity of deploying their competences in different research
groups. In order to sustain innovation in numerical simulation,
we believe there should be little, if any, distinction between
development and production codes; computer codes should
be easy to maintain, to understand by different generations of
young researchers, to modify, and extend; they should be easy
to use by the layman, as well as general and ?exible enough
to be enticing for a vast and diverse community of end users.
One easily understands that such con?icting requirements can
only be tempered, if anything, within organized and modular
software projects.
Software modularity also comes as a necessity when com-
plex problems in complex materials need to be tackled with
an array of different methods and techniques. Multiscale ap-
proaches, in particular, strive to combine methods with differ-
ent accuracy and scope to describe different parts of a com-
plex system, or phenomena occurring at different time and/or
length scales. Such approaches will require software packages
that can perform different kinds of computations on different
aspects of the same problem and/or different portions of the
same system, and that allow for interoperability or joint us-
age of the different modules. Different packages should at the
very least share the same input/output data formats; ideally
they should also share a number of mathematical and appli-
cation libraries, as well as the internal representation of some
of the data structures on which they operate. Individual re-
searchers or research groups ?nd it increasingly dif?cult to
meet all these requirements and to continue to develop and
maintain in-house software project of increasing complexity.
Thus, different and possibly collaborative solutions should be
sought.
A successful example comes from the software for sim-
ulations in quantum chemistry, that has often (but not al-
ways) evolved towards commercialization: the development
and maintenance of most well-known packages is devolved to
non-pro?t [5–8] or commercial companies [9–12]. The soft-
ware is released for a fee under some proprietary license that
may impose several restrictions to the availability of sources
(computer code in a high-level language) and to what can be
done with the software. This model has worked well, and is
also used by some of the leading development groups in the
condensed-matter electronic-structure community [13, 14],
while some proprietary projects allow for some free academic
usage of their products [14–19]. A commercial endeavor also
brings the distinctive advantage of a professional approach to
software development, maintenance, documentation, and sup-
port.
We believe however that a more interesting and fruitful al-
ternative can be pursued, and one that is closer to the spirit
of science and scienti?c endeavor, modeled on the experience
of the open-source software community. Under this model, a
large community of users has full access to the source code
and the development material, under the coordination of a
smaller group of core developers. In the long term, and in
the absence of entrenched monopolies, this strategy could be
more effective in providing good software solutions and in
nurturing a community engaged in providing those solutions,
as compared to the proprietary software strategy. In the case
of software for scienti?c usage, such an approach has the ad-
ditional, and by no means minor, advantage to be in line with
the tradition and best practice of science, that require repro-
ducibility of other people’s results, and where collaboration is
no less important than competition.
In this paper we will shortly describe our answer to the
above-mentioned problems, as embodied in our QUANTUM
ESPRESSO project (indeed, QUANTUM ESPRESSO stands
for opEn Source Package for Research in Electronic Structure,
Simulation, and Optimization). First, in Sec. 2, we describe
the guiding lines of our effort. In Sec. 3, we give an overview
of the current capabilities of QUANTUM ESPRESSO. In Sec.
4, we provide a short description of each software component
presently distributed within QUANTUM ESPRESSO. Finally,
Sec. 5 describes current developments and offers a perspec-
tive outlook. The Appendix sections discuss some of the more
speci?c technical details of the algorithms used, that have not
been documented elsewhere.
II. THE QUANTUM ESPRESSO PROJECT
QUANTUM ESPRESSO is an integrated suite of computer
codes for electronic-structure calculations and materials mod-
eling based on density-functional theory, plane waves basis
sets and pseudopotentials to represent electron-ion interac-
tions. QUANTUM ESPRESSO is free, open-source, software
distributed under the terms of the GNU General Public Li-
cense (GPL) [20].
The two main goals of this project are to foster method-
ological innovation in the ?eld of electronic-structure sim-
ulations and to provide a wide and diverse community of
end users with highly ef?cient, robust, and user-friendly soft-
ware implementing the most recent innovations in this ?eld.
Other open-source projects [21–25] exist, besides QUAN-
TUM ESPRESSO, that address electronic-structure calcula-
tions and various materials simulation techniques based on
them. Unlike some of these projects, QUANTUM ESPRESSO
does not aim at providing a single monolithic code able to per-
form several different tasks by specifying different input data
to a same executable. Our general philosophy is rather that of
an open distribution, i.e. an integrated suite of codes designed
to be interoperable, much in the spirit of a Linux distribution,
and thus built around a number of core components designed
and maintained by a small group of core developers, plus a
number of auxiliary/complementary codes designed, imple-
mented, and maintained by members of a wider community of
users. The distribution can even be redundant, with different
applications addressing the same problem in different ways;
at the end, the sole requirements to be added to the distribu-
tion are that: i) they are distributed under the same GPL li-
3
cence agreement [20] as other QUANTUM ESPRESSO com-
ponents; ii) they are fully interoperable with the other com-
ponents. Of course, they need to be scienti?cally sound, ver-
i?ed and validated. External contributors are encouraged to
join the QUANTUM ESPRESSO project, if they wish, while
maintaining their own individual distribution and advertise-
ment mode for their software (for instance, by maintaining
individual web sites with their own brand names [26]). To fa-
cilitate this, a web service called qe-forge [27], described
in the next subsection, has been recently put in place.
Interoperability of different components within QUANTUM
ESPRESSO is granted by the use of common formats for the
input, output, and work ?les. In addition, external contributors
are encouraged, but not by any means forced, to use the many
numerical and application libraries on which the core compo-
nents are built. Of course, this general philosophy must be
seen more as an objective to which a very complex software
project tends, rather than a starting point.
One of the main concerns that motivated the birth of the
QUANTUM ESPRESSO project is high performance, both
in serial and in parallel execution. High serial performance
across different architectures is achieved by the systematic
use of standardized mathematical libraries (BLAS, LAPACK
[28], and FFTW [29]) for which highly optimized implemen-
tations exist on many platforms; when proprietary optimiza-
tions of these libraries are not available, the user can compile
the library sources distributed with QUANTUM ESPRESSO.
Optimal performance in parallel execution is achieved through
the design of several parallelization levels, using sophisti-
cated communications algorithms, whose implementation of-
ten does not need to concern the developer, being embedded
in appropriate software layers. As a result the performance
of the key engines (PWscf and CP, see below) may scale on
massively parallel computers up to hundreds or thousands of
processors respectively.
The distribution is organized into a basic set of modules,
libraries, installation utilities, plus a number of directories,
each containing one or more executables, performing speci?c
tasks. The communications between the different executa-
bles take place via data ?les. We think that this kind of ap-
proach lowers the learning barrier for those who wish to con-
tribute to the project. The codes distributed with QUANTUM
ESPRESSO, including many auxiliary codes for the postpro-
cessing of the data generated by the simulations, are easy to
install and to use. The GNU configure and make util-
ities ensure a straightforward installation on many different
machines. Applications are run through text input ?les based
on Fortran namelists, that require the users to specify only an
essential but small subset of the many control variables avail-
able; a specialized graphical user interface (GUI) that is pro-
vided with the distribution facilitates this task for most com-
ponent programs. It is foreseen that in the near future the
design of special APIs (Application Programming Interfaces)
will make it easier to glue different components of the dis-
tribution together and with external applications, as well as
to interface them to other, custom-tailored, GUIs and/or com-
mand interpreters.
The QUANTUM ESPRESSO distribution is written,
mostly, in Fortran-95, with some parts in C or in Fortran-77.
Fortran-95 offers the possibility to introduce advanced pro-
gramming techniques without sacri?cing the performances.
Moreover Fortran is still the language of choice for high-
performance computing and it allows for easy integration of
legacy codes written in this language. A single source tree
is used for all architectures, with C preprocessor options se-
lecting a small subset of architecture-dependent code. Par-
allelization is achieved using the Message-Passing paradigm
and calls to standard MPI (Message Passing Interface) [30]
libraries. Most calls are hidden in a few routines that act as
an intermediate layer, accomplishing e.g. the tasks of sum-
ming a distributed quantity over processors, of collecting dis-
tributed arrays or distributing them across processors, and
to perform parallel three-dimensional Fast Fourier Transform
(FFT). This allows to develop straightforwardly and transpar-
ently new modules and functionalities that preserve the ef?-
cient parallelization backbone of the codes.
A. QE-forge
The ambition of the QUANTUM ESPRESSO project is not
limited to providing highly ef?cient and user-friendly soft-
ware for large-scale electronic-structure calculations and ma-
terials modeling. QUANTUM ESPRESSO aims at promot-
ing active cooperation between a vast and diverse commu-
nity of scientists developing new methods and algorithms in
electronic-structure theory and of end users interested in their
application to the numerical simulation of materials and de-
vices.
As mentioned, the main source of inspiration for the model
we want to promote is the successful cooperative experience
of the GNU/Linux developers’ and users’ community. One of
the main outcomes of this community has been the incorpora-
tion within the GNU/Linux operating system distributions of
third-party software packages, which, while being developed
and maintained by autonomous, and often very small, groups
of users, are put at the disposal of the entire community under
the terms of the GPL. The community, in turn, provides pos-
itive feedback and extensive validation by benchmarking new
developments, reporting bugs, and requesting new features.
These developments have largely bene?ted from the Source-
Forge code repository and software development service [31],
or by other similar services, such as RubyForge, Tigris.org,
BountySource, BerliOS, JavaForge, and GNU Savannah.
Inspired by this model, the QUANTUM ESPRESSO devel-
opers’ and users’ community has set up its own web portal,
named qe-forge [27]. The goal of qe-forge is to com-
plement the traditional web sites of individual scienti?c soft-
ware projects, which are passive instruments of information
retrieval, with a dynamical space for active content creation
and sharing. Its aim is to foster and simplify the coordination
and integration of the programming efforts of heterogeneous
groups and to ease the dissemination of the software tools thus
obtained.
qe-forge provides, through a user-friendly web inter-
face, an integrated development environment, whereby re-
4
searchers can freely upload, manage and maintain their own
software, while retaining full control over it, including the
right of not releasing it. The services so far available in-
clude source-code management software (CVS or SVN repos-
itory), mailing lists, public forums, bug tracking facilities,
download space, and wiki pages for projects’ documentation.
qe-forge is expected to be the main tool by which QUAN-
TUM ESPRESSO end users and external contributors can
maintain QUANTUM ESPRESSO-related projects and make
them available to the community.
III. SHORT DESCRIPTION OF QUANTUM ESPRESSO
QUANTUM ESPRESSO implements a variety of methods
and algorithms aimed at a chemically realistic modeling of
materials from the nanoscale upwards, based on the solution
of density-functional theory (DFT) [2, 3] problem, using a
plane waves (PW) basis set and pseudopotentials (PP) [32] to
represent electron-ion interactions.
The codes are constructed around the use of periodic
boundary conditions, which allows for a straightforward treat-
ment of in?nite crystalline systems, and an ef?cient con-
vergence to the thermodynamic limit for aperiodic but ex-
tended systems, such as liquids or amorphous materials. Fi-
nite systems are also treated using supercells; if required,
open-boundary conditions can be used through the use of the
density-countercharge method [33]. QUANTUM ESPRESSO
can thus be used for any crystal structure or supercell, and
for metals as well as for insulators. The atomic cores can be
described by separable [34] norm-conserving (NC) PPs [35],
ultra-soft (US) PPs [36], or by projector-augmented wave
(PAW) sets [37]. Many different exchange-correlation func-
tionals are available in the framework of the local-density
(LDA) or generalized-gradient approximation (GGA) [38],
plus advanced functionals like Hubbard U corrections and few
meta-GGA [39] and hybrid functionals [40–42]. The latter
is an area of very active development, and more details on
the implementation of hybrid functionals and related Fock-
exchange techniques are given in Appendix F.
The basic computations/simulations that can be performed
include:
• calculation of the Kohn-Sham orbitals and energies [43]
for isolated or extended/periodic systems, and of their
ground-state energies;
• complete structural optimizations of the microscopic
(atomic coordinates) and macroscopic (unit cell) de-
grees of freedom, using Hellmann-Feynman forces [44,
45] and stresses [46];
• ground states for magnetic or spin-polarized system, in-
cluding spin-orbit coupling [47] and non-collinear mag-
netism [48, 49];
• ab-initio molecular dynamics (MD), using either the
Car-Parrinello Lagrangian [50] or the Hellmann-
Feynman forces calculated on the Born-Oppenheimer
(BO) surface [51], in a variety of thermodynamical en-
sembles, including NPT variable-cell [52, 53] MD;
• density-functional perturbation theory (DFPT) [54–56],
to calculate second and third derivatives of the total en-
ergy at any arbitrary wavelength, providing phonon dis-
persions, electron-phonon and phonon-phonon interac-
tions, and static response functions (dielectric tensors,
Born effective charges, IR spectra, Raman tensors);
• location of saddle points and transition states via
transition-path optimization using the nudged elastic
band (NEB) method [57–59];
• ballistic conductance within the Landauer-Büttiker the-
ory using the scattering approach [60];
• generation of maximally localized Wannier functions
[61, 62] and related quantities;
• calculation of nuclear magnetic resonance (NMR) and
electronic paramagnetic resonance (EPR) parameters
[63, 64];
• calculation of K-edge X-ray absorption spectra [65].
Other more advanced or specialized capabilities are described
in the next sections, while ongoing projects (e.g. time-
dependent DFT, many-body perturbation theory) are men-
tioned in the last section. Selected applications were described
in Ref [66]. Several utilities for data postprocessing and inter-
facing to advanced graphic applications are available, allow-
ing e.g. to calculate scanning tunneling microscopy (STM)
images [67], the electron localization function (ELF) [68],
Löwdin charges [69], the density of states (DOS), and pla-
nar [70] or spherical averages of the charge and spin densities
and potentials.
A. Data ?le format
The interoperability of different software components
within a complex project such as QUANTUM ESPRESSO re-
lies on the careful design of ?le formats for data exchange. A
rational and open approach to data ?le formats is also essential
for interfacing applications within QUANTUM ESPRESSO
with third-party applications, and more generally to make the
results of lengthy and expensive computer simulations ac-
cessible to, and reproducible by, the scienti?c community at
large. The need for data ?le formats that make data exchange
easier than it is now is starting to be widely appreciated in the
electronic-structure community. This problem has many as-
pects and likely no simple, "one-size-?ts-all", solution. Data
?les should ideally be
• extensible: one should be able to add some more infor-
mation to a ?le without breaking all codes that read that
?le;
• self-documenting: it should be possible to understand
the contents of a ?le without too much effort;
5
• ef?cient: with data size in the order of GBytes for
large-scale calculations, slow or wasteful I/O should be
avoided.
The current trend in the electronic-structure community seems
to be the adoption of one of the following approaches:
• Structured ?le formats, notably Hierarchical Data For-
mat (HDF) [71] and network Common Data Form
(netCDF) [72], that are widely used since years in other
communities;
• ?le formats based on the Extensible Markup Language
(XML) [73].
It is unlikely that a common, standardized data formats will
ever prevail in our community. We feel that we should fo-
cus, rather than on standardization, on an approach that al-
lows an easy design and usage of simple and reliable con-
verters among different data formats. Prompted by these con-
siderations, QUANTUM ESPRESSO developers have opted
for a simple solution that tries to combine the advantages of
both the above-mentioned approaches. A single ?le contain-
ing all the data of a simulation is replaced by a data direc-
tory, containing several ?les and subdirectories, much in the
same way as it is done in the Mac OS X operating system.
The “head” ?le contains data written with ordinary Fortran
formatted I/O, identi?ed by XML tags. Only data of small
size, such as atomic positions, parameters used in the calcu-
lation, one-electron and total energies, are written in the head
?le. Data of potentially large size, such as PW coef?cients of
Kohn-Sham orbitals, charge density, potentials, are present as
links to separate ?les, written using unformatted Fortran I/O.
Data for each k-point is written to a separate subdirectory.
A lightweight library called iotk, standing for Input/Output
ToolKit [74], is used to read and write the data directory.
Another problem affecting interoperability of PW-PP codes
is the availability of data ?les containing atomic PP’s — one
of the basic ingredients of the calculation. There are many
different types of PP’s, many different codes generating PP’s
(see e.g. Ref [75–77]), each one with its own format. Again,
the choice has fallen on a simple solution that makes it easy
to write converters from and to the format used by QUANTUM
ESPRESSO. Each atomic PP is contained in a formatted ?le
(ef?ciency is not an issue here), described by a XML-like syn-
tax. The resulting format has been dubbed Uni?ed Pseudopo-
tential File (UPF). Several converters from other formats to
the UPF format are available in QUANTUM ESPRESSO.
IV. PARALLELIZATION
Keeping the pace with the evolution of high-end supercom-
puters is one of the guiding lines in the design of QUANTUM
ESPRESSO, with a signi?cant effort being dedicated to port-
ing it to the latest available architectures. This effort is moti-
vated not only by the need to stay at the forefront of architec-
tural innovation for large to very-large scale materials science
simulations, but also by the speed at which hardware features
speci?cally designed for supercomputers ?nd their way into
commodity computers.
The architecture of today’s supercomputers is character-
ized by multiple levels and layers of parallelism: the bot-
tom layer is the one affecting the instruction set of a single
core (simultaneous multithreading, hyperthreading); then one
has parallel processing at processor level (many CPU cores
inside a single processor sharing caches) and at node level
(many processors sharing the same memory inside the node);
at the top level, many nodes are ?nally interconnected with
a high-performance network. The main components of the
QUANTUM ESPRESSO distribution are designed to exploit
this highly structured hardware hierarchy. High performance
on massively parallel architectures is achieved by distributing
both data and computations in a hierarchical way across avail-
able processors, ending up with multiple parallelization levels
[78] that can be tuned to the speci?c application and to the
speci?c architecture. This remarkable characteristic makes it
possible for the main codes of the distribution to run in parallel
on most or all parallel machines with very good performance
in all cases.
More in detail, the various parallelization levels are geared
into a hierarchy of processor groups, identi?ed by different
MPI communicators. In this hierarchy, groups implement-
ing coarser-grained parallel tasks are split into groups im-
plementing ?ner-grained parallel tasks. The ?rst level is im-
age parallelization, implemented by dividing processors into
n
image
groups, each taking care of one or more images (i.e.
a point in the con?guration space, used by the NEB method).
The second level is pool parallelization, implemented by fur-
ther dividing each group of processors into n
pool
pools of
processors, each taking care of one or more k-points (see
Appendix D). The third level is plane-wave parallelization,
implemented by distributing real- and reciprocal-space grids
across the n
PW
processors of each pool. The ?nal level is task
group parallelization [79], in which processors are divided
into n
task
task groups of n
FFT
= n
PW
/n
task
processors,
each one taking care of different groups of electron states to
be Fourier-transformed, while each FFT is parallelized inside
a task group. A further paralellization level, linear-algebra,
coexists side-to-side with plane-wave parallelization, i.e. they
take care of different sets of operations, with different data dis-
tribution. Linear-algebra parallelization is implemented both
with custom algorithms and using ScaLAPACK [80], which
on massively parallel machines yield much superior perfor-
mances. Table I contains a summary of the ?ve levels cur-
rently implemented. With the recent addition of the two last
levels, most parallelization bottlenecks have been removed,
both computations and data structures are fully distributed,
scalability on parallel machines is only limited by the phys-
ical sizes of the system being simulated. Scalability does not
yet extend to tens of thousands of processors as in specially-
crafted code QBox[81], but excellent scalability on up to 4800
processors has been demonstrated (see Fig. 1) even for cases
where coarse-grained parallelization doesn’t help.
Scalability is thus guaranteed for large-scale simulations.
This being said, it is obvious that the size and speci?c nature
of the speci?c application sets quite natural limits to the max-
6
TABLE I: Summary of parallelization levels in QUANTUM ESPRESSO.
group distributed quantities communications performance
image NEB images very low linear CPU scaling,
fair to good load balancing;
does not distribute RAM
pool k-points low almost linear CPU scaling,
fair to good load balancing;
does not distribute RAM
plane-wave plane waves, G-vector coef?cients, high good CPU scaling,
R-space FFT arrays good load balancing,
distributes most RAM
task FFT on electron states high improves load balancing
linear algebra subspace Hamiltonians very high improves scaling,
and constraints matrices distributes more RAM
1000
1500
2000
2500
3000
3500
4000
4500
5000
0 1000 2000 3000 4000 5000
C
P
U
t
i
m
e
(
s
)
Number of CPUs
Scalability for > 1000 processors
PSIWAT
CNT (1)
CNT (2)
FIG. 1: Scalability as a function of the number of processors for
large-scale calculations. PSIWAT: gold surface covered by thiols in
interaction with water, 10.59 × 20.53 × 32.66 Å
3
cell, 587 atoms,
2552 electrons, PWscf code with n
pool
= 4 and n
tg
= 4 on A Cray
XT 4. CNT: porphyrin-functionalized nanotube, 1532 atoms, 5232
electrons; case (1), PWscf on a Cray XT 4; case (2), CP on a Cray
XT3(2) codes. CPU times for PSIWAT and CNT (2) were rescaled
so that they to fall in the same range as for CNT (1).
imum number of processors up to which the performances of
the various codes are expected to scale. For instance, the num-
ber of images in a NEB calculation sets a natural limit to the
level of image groups, or the number of electronic bands sets
a limit for the parallelization of the linear algebra operations.
Moreover some numerical algorithms scale better than others.
For example, the use of norm-conserving pseudopotentials al-
lows for a better scaling than ultrasoft pseudopotentials for
a same system, because a larger plane wave basis set and a
larger real- and reciprocal-space grids are required in the for-
mer case. On the other hand, using ultrasoft pseudopotentials
is generally faster because the use of a smaller basis set is
obviously more ef?cient, even though the overall parallel per-
formance may be not as good.
The efforts of the QUANTUM ESPRESSO developers’
50
100
150
200
250
300
350
0 100 200 300 400 500
C
P
U
t
i
m
e
(
s
)
Number of CPUs
CP scalability
BG/P 1 task
BG/P 4 tasks
BG/P 8 tasks
Altix 1 task
Altix 4 tasks
Altix 8 tasks
FIG. 2: CPU time t
cpu
(s) per electronic time step (CP code) on
an IBM Blue Gene/P (BG/P) and a SGI Altix for a fragment of an
A??peptide in water containing 838 atoms and 2311 electrons in a
22.1×22.9×19.9Å
3
cell, as a function of the number of processors,
for different numbers n
task
of task groups.
team are not limited to the performance on massively paral-
lel architectures. Simulations on systems containing several
hundreds of atoms are by now quite standard (see Fig. 2
for an example). Special attention is also paid to optimize
the performances for simulations of intermediate size (on sys-
tems comprising from several tens to a few hundreds inequiv-
alent atoms), to be performed on medium-size clusters, read-
ily available to many groups [82]. In particular, the QUAN-
TUM ESPRESSO developers’ team is now working to better
exploit new hardware trends, particularly in the ?eld of mul-
ticore architectures. By the time this paper will be printed,
mixed OpenMP-MPI parallelization will be available for en-
hanced performance on modern multicore CPU’s. Looking
ahead, future developments will likely focus on hybrid system
with hardware accelerators (GPUs and cell co-processors).
7
V. QUANTUM ESPRESSO PACKAGES
The complete QUANTUM ESPRESSO distribution is rel-
ative large: about 240,000 lines of code, excluding copies of
the external libraries. With such a sizable code basis, mod-
ularization becomes necessary. QUANTUM ESPRESSO is
presently divided into several executables, performing differ-
ent types of calculations, although some of them have over-
lapping functionalities. Typically there is a single set of func-
tions/subroutines or a single Fortran 90 module that performs
each speci?c task (e.g. matrix diagonalizations, or potential
updates), but there are still important exceptions to this rule,
re?ecting the different origin and different styles of the orig-
inal components. QUANTUM ESPRESSO has in fact been
built out of the merge and re-engineering of different pack-
ages, that were independently developed for several years. In
the following, the main components are brie?y described.
A. PWscf
PWscf implements an iterative approach to reach self-
consistency, using at each step iterative diagonalization tech-
niques, in the framework of the plane-wave pseudopotential
method. An early version of PWscf is described in Ref [83].
Both separable NC-PPs and US-PPs are implemented; re-
cently, also the projector augmented-wave method [37] has
been added, largely following the lines of Ref [84] for its
implementation. In the case of US-PPs, the electronic wave
functions can be made smoother at the price of having to aug-
ment their square modulus with additional contributions to re-
cover the actual physical charge densities. For this reason,
the charge density is more structured than the square of the
wavefunctions, and requires a larger energy cutoff for its plane
wave expansion (typically, 6 to 12 times larger; for a NC-PP,
a factor of 4 would be mathematically suf?cient). Hence, dif-
ferent real-space Fourier grids are introduced - a "soft" one
that represents the square of electronic wave functions, and a
"hard" one that represents the charge density [82, 85]. The
augmentation terms can be added either in reciprocal space
(using an exact but expensive algorithm) or directly in real
space (using an approximate but faster algorithm that exploits
the local character of the augmentation charges).
PWscf can use the established LDA and GGA exchange-
correlation functionals, including spin-polarization within the
scheme proposed in Ref [86] and can treat non-collinear
magnetism[48, 49] as e.g. induced by relativistic effects (spin-
orbit interactions) [87, 88] or by complex magnetic interac-
tions (e.g. in the presence of frustration). DFT + Hubbard U
calculations [89] are implemented for a simpli?ed (“no-J”)
rotationally invariant form [90] of the Hubbard term. Other
advanced functionals include TPSS meta-GGA[39], function-
als with ?nite-size corrections [91], and the PBE0 [40] and
B3LYP [41, 42] hybrids.
Self-consistency is achieved via the modi?ed Broyden
method of Ref [92], with some further re?nements that are
detailed in the Appendix. The sampling of the Brillouin Zone
(BZ) can be performed using either special [93, 94] k-points
provided in input or those automatically calculated starting
froma uniform grid. Crystal symmetries are automatically de-
tected and exploited to reduce computational costs, by restrict-
ing the sampling of the BZ to the irreducible wedge alone.
When only the ? point (k = 0) is used, advantage is taken
of the real character of the orbitals, allowing to store just half
of the Fourier components. BZ integrations in metallic sys-
tems are dealt with smearing/broadening techniques, such as
Fermi-Dirac, Gaussian, Methfessel-Paxton [95], and Marzari-
Vanderbilt cold smearing [96]; the tetrahedron method [97] is
also implemented. Note that ?nite-temperature effects on the
electronic properties can be easily accounted for by using the
Fermi-Dirac smearing not as a mathematical device, but as a
practical way of implementing the Mermin density-functional
approach [98].
Structural optimizations are performed using the Broyden-
Fletcher-Goldfarb-Shanno (BFGS) algorithm [99–101] or
damped dynamics; these can involve both the internal, micro-
scopic degrees of freedom (i.e. the atomic coordinates) and/or
the macroscopic ones (shape and size of the unit cell). The
calculation of minimum-energy paths, activation energies, and
transition states uses the Nudged Elastic Band (NEB) method
[57]. Potential energy surfaces as a function of suitably cho-
sen collective variables can be studied using Laio-Parrinello
metadynamics [102].
Microcanonical (NVE) MD is performed on the BO sur-
face, i.e. achieving electron self-consistency at each time step,
using the Verlet algorithm[103]. Canonical (NVT) dynamics
can be performed using velocity rescaling, or Anderson’s or
Berendsen’s thermostats [104]. Constant-pressure (NPT) MD
is performed by adding additional degrees of freedom for the
cell size and volume, using either the Parrinello-Rahman La-
grangian [105] or the invariant Lagrangian of Wentzcovitch
[53].
The effects of ?nite macroscopic electric ?elds on the elec-
tronic structure of the ground state can be accounted for ei-
ther through the method of Ref [106, 107] based on the Berry
phase, or (for slab geometries only) through a sawtooth ex-
ternal potential [108, 109]. A quantum fragment can be em-
bedded in a complex electrostatic environment that includes a
model solvent [110] and a counterion distribution [111], as is
typical of electrochemical systems.
B. CP
The CP code is the massively-parallel module for Car-
Parrinello ab-initio MD. CP can use both NC PPs [112] and
US PPs [85, 113]. In the latter case, the electron density
is augmented through a Fourier interpolation scheme in real
space (“box grid”) [82, 85] that is particular ef?cient for large
scale calculations. CP implements the same functionals as in
PWscf, with the exception of hybrid functionals; a simpli-
?ed one-electron self-interaction correction (SIC)[114] is also
available. The Car-Parrinello Lagrangian can be augmented
with Hubbard U corrections [115], or Hubbard-based penalty
functionals to impose arbitrary oxidation states [116].
Since the main applications of CP are for large systems
8
without translational symmetry (e.g. liquids, amorphous ma-
terials), Brillouin zone sampling is restricted to the ? point of
the supercell, allowing for real instead of complex wavefunc-
tions. Metallic systems can be treated in the framework of
“ensemble DFT” [117].
In the Car-Parrinello algorithm, microcanonical (NVE) MD
is performed on both electronic and nuclear degrees of free-
dom, treated on the same footing, using the Verlet algorithm.
The electronic equations of motion are accelerated through a
preconditioning scheme [118]. Constant-pressure (NPT) MD
is performed using the Parrinello-Rahman Lagrangian [105]
and additional degrees of freedom for the cell. Nosé-Hoover
thermostats [119] and Nosé-Hoover chains [120] allow to per-
form simulations in the different canonical ensembles.
CP can also be used to directly minimize the energy func-
tional to self-consistency while keeping the nuclei ?xed, or to
perform structural minimizations of nuclear positions, using
the “global minimization” approaches of in Refs. [121, 122],
and damped dynamics or conjugate-gradients on the elec-
tronic or ionic degrees of freedom. It can also perform NEB
and metadynamics calculations.
Finite homogeneous electric ?elds can be accounted for us-
ing the Berry’s phase method, adapted to systems with the
? point only [106]. This advanced feature can be used in
combination with MD to obtain the infrared spectra of liq-
uids [106, 123], the low- and high-frequency dielectric con-
stants [106, 124] and the coupling factors required for the cal-
culation of vibrational properties, including infrared, Raman
[125–127], and hyper-Raman [128] spectra.
C. PHonon
The PHonon code implements density-functional pertur-
bation theory (DFPT) [54–56] for the calculation of second-
and third-order derivatives of the energy with respect to
atomic displacements and to electric ?elds. The global min-
imization approach [129, 130] is used for the special case
of normal modes in ?nite (molecular) systems, where no BZ
sampling is required, while a self-consistent procedure [55] is
used in the general case, with the distinctive advantage that
the response to a perturbation of any arbitrary wavelength
can be calculated with a computational cost that is propor-
tional to that of the unperturbed system. Thus, response at
any wavevector, including very long-wavelength ones, can be
inexpensively calculated. This latter approach, and the tech-
nicalities involved in the calculation of effective charges and
interatomic force constants, are described in detail in Refs.
[55, 131].
Symmetry is fully exploited in order to reduce the amount
of computation. Lattice distortions transforming according
to irreducible representations of small dimensions are gener-
ated ?rst. The charge-density response to these lattice distor-
tions is then sampled at a number of discrete k-points in the
BZ, which is reduced according to the symmetry of the small
group of the phonon wavevector q. The grid of q needed for
the calculation of interatomic force constants reduces to one
wavevector per star: the dynamical matrices at the other q
vectors in the star are generated using the symmetry opera-
tions of the crystal. This approach allows us to speed up the
calculation without the need to store too much data for sym-
metrization.
The calculation of second-order derivatives of the energy
works also for US PP [132, 133] and for all GGA ?avors
[134, 135] used in PWscf and in CP. The extension of
PHonon to PAW [136], to noncollinear magnetism and to
fully relativistic US PPs which include spin-orbit coupling
[137] will be available by the time this paper will be printed.
Advanced features of the PHonon code include the calcu-
lation of third order derivatives of the energy, such as electron-
phonon or phonon-phonon interaction coef?cients. Electron-
phonon interactions are straightforwardly calculated from the
response of the self-consistent potential to a lattice distortion.
This involves a numerically-sensitive “double-delta” integra-
tion at the Fermi energy, that is performed using interpolations
on a dense k-point grid. Interpolation techniques based on
Wannier functions [138] will speed up considerably these cal-
culations. The calculation of the anharmonic force constants
from third-order derivatives of the electronic ground-state en-
ergy is described in Ref. [139] and is performed by a separate
code called d3. Static Raman coef?cients are calculated us-
ing the second-order response approach of Refs. [140, 141].
Both third-order derivatives and Raman-coef?cients calcula-
tions are currently implemented only for NC PPs.
D. atomic
The atomic code performs three different tasks: i) solu-
tion of the self-consistent all-electron radial Kohn-Shamequa-
tions (with a Coulomb nuclear potential and spherically sym-
metric charge density); ii) generation of NC PPs, of US PPs,
or of PAW data-sets; iii) test of the above PPs and data-sets.
The three tasks can be either separately executed or performed
in a single run. Three different all-electron equations are avail-
able: i) the non relativistic radial Kohn and Sham equations,
ii) the scalar relativistic approximation to the radial Dirac
equations [142], iii) the radial Dirac-like equations derived
within relativistic density functional theory [143, 144]. For i)
and ii) atomic magnetism is dealt within the local spin density
approximation, i.e. assuming an axis of magnetization. The
atomic code uses the same exchange and correlation energy
routines of PWscf and can deal with the same functionals.
The code is able to generate NC PPs directly in separable
form (also with multiple projectors per angular momentum
channel) via the Troullier-Martins [145] or the Rappe-Rabe-
Kaxiras-Joannopoulos [146] pseudization. US PPs can be
generated by a two-step pseudization process, starting from
a NC PPs, as described in Ref. [147], or using the solutions
of the all-electron equation and pseudizing the augmentation
functions [85]. The latter method is used also for the PAW
data-set generation. The generation of fully relativistic NC
and US PPs including spin-orbit coupling effects is also avail-
able. Converters are available to translate pseudopotentials
encoded in different formats (e.g. according to the Fritz-
Haber [75] or Vanderbilt [76] conventions) into the UPF for-
9
mat adopted by QUANTUM ESPRESSO.
Transferability tests can be made simultaneously for sev-
eral atomic con?gurations with or without spin-polarization,
by solving the non relativistic radial Kohn-Sham equations
generalized for separable nonlocal PPs and for the presence
of an overlap matrix.
E. PWcond
The PWcond code implements the scattering approach pro-
posed by Choi and Ihm [60] for the study of coherent electron
transport in atomic-sized nanocontacts within the Landauer-
Büttiker theory. Within this scheme the linear response bal-
listic conductance is proportional to the quantum-mechanical
electron transmission at the Fermi energy for an open quan-
tum system consisting of a scattering region (e.g., an atomic
chain or a molecule with some portions of left and right leads)
connected ideally from both sides to semi-in?nite metallic
leads. The transmission is evaluated by solving the Kohn-
Sham equations with the boundary conditions that an electron
coming fromthe left lead and propagating rightwards gets par-
tially re?ected and partially transmitted by the scattering re-
gion. The total transmission is obtained by summing all trans-
mission probabilities for all the propagating channels in the
left lead. As a byproduct of the method, the PWcond code
provides the complex band structures of the leads, that is the
Bloch states with complex k
z
in the direction of transport, de-
scribing wave functions exponentially growing or decaying in
the z direction. The original method formulated with NC PPs
has been generalized to US PPs both in the scalar relativistic
[148] and in the fully relativistic forms [149].
F. GIPAW
The GIPAW code allows for the calculation of physical pa-
rameters measured by i) nuclear magnetic resonance (NMR)
in insulators (the electric ?eld gradient (EFG) tensors and the
chemical shift tensors), and by ii) electronic paramagnetic res-
onance (EPR) spectroscopy for paramagnetic defects in solids
or in radicals (the hyper?ne tensors and the g-tensor). The
code also computes the magnetic susceptibility of nonmag-
netic insulators. The code is based on the PW-PP method,
and uses many subroutines of PWscf and of PHonon. The
code is currently restricted to NC PPs. All the NMR and
EPR parameters depend on the detailed shape of the elec-
tronic wave-functions near the nuclei and thus require the re-
construction of the all-electron wave-functions from the PP
wave-functions. For the properties de?ned at zero external
magnetic ?eld, namely the EFG and the hyper?ne tensors,
such reconstruction is performed as a post-processing step of
a self-consistent calculation using the PAW reconstruction, as
described for the EFG in Ref.[150] and for the hyper?ne ten-
sor in Ref.[151]. The g-tensor, the NMR chemical shifts and
the magnetic susceptibility are obtained from the orbital lin-
ear response to an external uniform magnetic ?eld. In pres-
ence of a magnetic ?eld the PAW method is no more gauge-
and translationally invariant. Gauge and translational invari-
ances are restored by using the gauge including projector aug-
mented wave (GIPAW) method [63, 64] both i) to describe
in the PP Hamiltonian the coupling of orbital degrees of free-
dom with the external magnetic ?eld, and ii) to reconstruct the
all-electron wave-functions, in presence of the external mag-
netic ?eld. In addition, the description of a uniform magnetic
?eld within periodic boundary condition is achieved by con-
sidering the long wave-length limit of a sinusoidally modu-
lated ?eld in real space [152, 153]. The NMR chemical shifts
are computed following the method described in Ref.[63], the
g-tensor following Ref.[154] and the magnetic susceptibility
following Refs.[63, 152]. Recently, a “converse” approach
to calculate chemical shifts has also been introduced [155],
based on recent developments on the Berry-phase theory of or-
bital magnetization; since it does not require a linear-response
calculation, it can be straightforwardly applied to arbitrarily
complex exchange-correlation functionals, and to very large
systems, albeit at computational costs that are proportional to
the number of chemical shifts that need to be calculated.
G. XSPECTRA
The XSPECTRA code allows for the calculation of K-edge
X-ray absorption spectra. The code calculates both the dipo-
lar and the quadrupolar matrix elements. The code uses the
self-consistent charge density produced by PWscf and acts
as a post-processing tool. The all-electron wavefunction is
reconstructed using the PAW method and its implementation
in the GIPAW code. The presence of a core-hole in the ?nal
state of the X-ray absorption process is simulated by using a
pseudopotential for the absorbing atom with a hole in the 1s
state. The calculation of the charge density is performed on a
supercell with one absorbing atom. From the self-consistent
charge density, the X-ray absorption spectra are obtained us-
ing the Lanczos method and a continuous fraction expansion
[65]. The advantage of this approach is that it is not necessary
to calculate many empty states even if features of the spec-
trum at very high energy are required. Correlation effects can
be simulated in a mean-?eld way using the Hubbard U cor-
rection [90] that has been included in the XSPECTRA code in
Ref. [156]. Currently the code works only with NC PPs and is
limited to collinear magnetism. Its extension to noncollinear
magnetismand to fully relativistic US PPs which include spin-
orbit coupling is under development.
H. Wannier90
Wannier90 [26, 157] is a code that calculates maximally-
localized Wannier functions in insulators or metals—
according to the algorithms described in Refs. [61, 62]—and
a number of properties that can be conveniently expressed in
a Wannier basis. The code is developed and maintained inde-
pendently by a Wannier development group [157] and can be
taken as a representative example of the philosophy described
earlier, where a project maintains its own individual distribu-
10
tion but provides full interoperability with the core compo-
nents of QUANTUM ESPRESSO in this case PWscf or CP.
These codes are in fact used as “quantum engines” to produce
the data onto which Wannier90 operates. The need to pro-
vide transparent protocols for interoperability has in turn fa-
cilitated the interfacing of wannier90 with other quantum
engines [14, 21], fostering a collaborative engagement with
the broader electronic-structure community that is also in the
spirit of QUANTUM ESPRESSO .
Wannier90 requires as input the scalar products between
wavefunctions at neighboring k-points, where these latter
form uniform meshes in the Brillouin zone. Often, it is
also convenient to provide scalar products between wavefunc-
tions and trial, localized real-space orbitals—these are used to
guide the localization procedure towards a desired, physical
minimum. As such, the code is not tied to a representation
of the wavefunctions in any particular basis—for PWscf and
CP a postprocessing utility is in charge of calculating these
scalar products using the plane-wave basis set of QUANTUM
ESPRESSO and either NC-PPs or US-PPs. Whenever ?
sampling is used, the simpli?ed algorithm of Ref. [158] is
adopted.
Besides calculating maximally localized Wannier func-
tions, the code is able to construct the Hamiltonian matrix
in this localized basis, providing a chemically accurate, and
transferable, tight-binding representation of the electronic-
structure of the system. This, in turn, can be used to construct
Green’s functions and self-energies for ballistic transport cal-
culations [159, 160], to determine the electronic-structure and
density-of-states of very large scale structures [160], to inter-
polate accurately the electronic band structure (i.e. the Hamil-
tonian) across the Brillouin zone [160, 161], or to interpolate
any other operator [161]. These latter capabilities are espe-
cially useful for the calculation of integrals that depend sen-
sitively on a submanifold of states; common examples come
from properties that depend sensitively on the Fermi surface,
such as electronic conductivity, electron-phonon couplings
Knight shifts, or the anomalous Hall effect. A related by-
product of Wannier90 is the capability of downfolding a
selected, physically signi?cant manifold of bands into a mini-
mal but accurate basis, to be used for model Hamiltonians that
can be treated with complex many-body approaches.
I. PostProc
The PostProc module contains a number of codes for
postprocessing and analysis of data ?les produced by PWscf
and CP. The following operations can be performed:
• Interfacing to graphical and molecular graphics appli-
cations. Charge and spin density, potentials, ELF [68]
and STM images [67] are extracted or calculated and
written to ?les that can be directly read by most com-
mon plotting programs, like xcrysden [162] and VMD
[163].
• Interfaces to other codes that use DFT results from
QUANTUM ESPRESSO for further calculations, such
as e.g.: pw2wannier90, an interface to the
wannier90 library and code [26, 157] (also in-
cluded in the QUANTUM ESPRESSO distribution);
pw2casino.f90, an interface to the casino quan-
tum Monte Carlo code [164]; wannier_ham.f90, a
tool to build a tight-binding representation of the KS
Hamiltonian to be used by the dmft code [165] (avail-
able at the qe-forge site); pw_export.f90, an in-
terface to the GW code SaX [166]; pw2gw.f90, an
interface to code DP [167] for dielectric property calcu-
lations, and to code EXC [168] for excited-state proper-
ties.
• Calculation of various quantities that are useful for the
analysis of the results. In addition to the already men-
tioned ELF and STM, one can calculate projections
over atomic states (e.g. Löwdin charges [69]), DOS and
Projected DOS (PDOS), planar and spherical averages,
and the complex macroscopic dielectric function in the
random-phase approximation (RPA).
J. PWgui
PWgui is the graphical user interface (GUI) for the
PWscf, PHonon, d3 and atomic codes as well as for
some of the main post-processing routines (e.g. pp.x and
projwfc.x). PWgui is an input ?le builder whose main
goal is to lower the learning barrier for the newcomer, who
has to struggle with the input syntax. Its event-driven mecha-
nism automatically adjusts the display of required input ?elds
(i.e. enables certain sets of widgets and disables others) to the
speci?c cases selected. It enables a preview of the format of
the (required) input ?le records for a given type of calculation.
The input ?les created by PWgui are syntactically correct (al-
though they can still be physically meaningless). It is possible
to upload previously generated input ?les for syntax checking
and/or to modify them. It is also possible to run calculations
from within the PWgui. In addition, PWgui can also use
the external xcrysden program [162] for the visualization
of molecular and/or crystal structures from the speci?ed in-
put data and for the visualization of properties (e.g. charge
densities or STM images).
As the QUANTUM ESPRESSO codes evolve, the input
?le syntax expands as well. This implies that PWgui has
to be continuously adapted. To effectively deal with such is-
sue, PWgui uses the GUIB concept [169]. GUIB builds on
the consideration that the input ?les for numerical simulation
codes have a rather simple structure and it exploits this sim-
plicity by de?ning a special meta-language with two purposes:
The ?rst is to de?ne the input-?le syntax, and the second is to
simultaneously automate the construction of the GUI on the
basis of such a de?nition.
A similar strategy has been recently adopted for the de-
scription of the QUANTUM ESPRESSO input ?le formats.
A single de?nition/description of a given input ?le serves as i)
a documentation per-se, ii) as a PWgui help documentation,
11
FIG. 3: Snapshot of the PWgui application. Left: PWgui’s main window; right: preview of speci?ed input data in text mode.
and iii) as a utility to synchronize the PWgui with up-to-date
input ?le formats.
VI. PERSPECTIVES AND OUTLOOK
Further developments and extensions of QUANTUM
ESPRESSO will be driven by the needs of the community
using it and working on it. Many of the soon-to-come ad-
ditions will deal with excited-state calculations within time-
dependent DFT (TDDFT [170, 171]) and/or many-body per-
turbation theory [172]. A new approach to the calculation of
optical spectra within TDDFT has been recently developed
[173], based on a ?nite-frequency generalization of density-
functional perturbation theory [54, 55], and implemented
in QUANTUM ESPRESSO. Another important development
presently under way is an ef?cient implementation of GW cal-
culations for large systems (whose size is of the order of a few
hundreds inequivalent atoms) [174]. The implementation of
ef?cient algorithms for calculating correlation energies at the
RPA level is also presently under way [175–177]. It is fore-
seen that by the time this paper will appear, many of these
developments will be publicly released.
It is hoped that many newfunctionalities will be made avail-
able to QUANTUM ESPRESSO users by external groups who
will make their own software compatible/interfaceable with
QUANTUM ESPRESSO. At the time of the writing of the
present paper, third-party scienti?c software available to the
QUANTUM ESPRESSO users’ community include: yambo,
a general-purpose code for excited-state calculations within
many-body perturbation theory [178]; casino, a code for
electronic-structure quantum Monte Carlo simulations [164];
wannier90, a code and a library for generating maximally
localized Wannier functions [26, 61]; want, a code for the
simulation of ballistic transport in nanostructures, based on
Wannier functions [179]; xcrysden, a molecular graphics
application, specially suited for periodic structures [162]. The
qe-forge portal is expected to burst the production and
availability of third-party software compatible with QUAN-
TUM ESPRESSO. Among the projects already available, or
soon-to-be available, on qe-forge, we mention: SaX [166],
an open-source project implementing state-of-the-art many-
body perturbation theory methods for excited states; dmft
[165], a code to perform Dynamical Mean-Field Theory cal-
culations on top of a tight-binding representation of the DFT
band structure; qha, a set of codes for calculating thermal
properties of materials within the quasi-harmonic approxima-
tion [180]; pwtk, a fully functional Tcl scripting interface to
PWscf [181].
Efforts towards better interoperability with third-party soft-
ware will be geared towards releasing accurate speci?cations
for data structures and data ?le formats and providing inter-
faces to and from other codes and packages used by the sci-
enti?c community. Further work will be also devoted to the
extension to the US-PPs and PAW schemes of the parts of
QUANTUM ESPRESSO that are now limited to NC-PPs.
The increasing availability of massively parallel machines
will likely lead to an increased interest towards large-scale cal-
culations. The ongoing effort in this ?eld will continue. A
special attention will be paid to the requirements imposed by
the architecture of the new machines, in particular multicore
CPUs, for which a mixed OpenMP-MPI approach seems to
be the only viable solution yielding maximum performances.
Grid computing and the commoditization of computer cluster
will also lead to great improvements in high-throughput cal-
culations for materials design and discovery.
12
VII. ACKNOWLEDGMENTS
The QUANTUM ESPRESSO project is an initiative of the
CNR-INFM DEMOCRITOS National Simulation Center in
Trieste (Italy) and its partners, in collaboration with MIT,
Princeton University, the University of Minnesota, the Ecole
Polytechnique Fédérale de Lausanne, the Université Pierre et
Marie Curie in Paris, the Jožef Stefan Institute in Ljubljana,
and the S3 research center in Modena. Many of the ideas
embodied in the QUANTUM ESPRESSO codes have ?our-
ished in the very stimulating environment of the International
School of Advanced Studies, SISSA, where Democritos is
hosted, and have bene?tted from the ingenuity of generations
of students and young postdocs. SISSA and the National Su-
percomputing Center CINECA are currently providing valu-
able support to the QUANTUM ESPRESSO project.
VIII. APPENDIX
This appendix contains the description of some algorithms
used in QUANTUM ESPRESSO that have not been docu-
mented elsewhere.
A. Self-consistency
The problem of ?nding a self-consistent solution to the
Kohn-Sham equations can be recast into the solution of a non-
linear problem
x = F[x], x = (x
1
, x
2
, . . . , x
N
), (1)
where vector x contains the N Fourier components or real-
space values of the charge density ? or the Kohn-Sham po-
tential V (the sum of Hartree and exchange-correlation poten-
tials); F[x
(in)
] is a functional of the input charge density or
potential x
(in)
, yielding the output vector x
(out)
via the so-
lution of Kohn-Sham equations. A solution can be found via
an iterative procedure. PWscf uses an algorithm based on
the modi?ed Broyden method [92] in which x contains the
components of the charge density in reciprocal space. Mixing
algorithms typically ?nd the optimal linear combination of a
few x
(in)
from previous iterations, that minimizes some suit-
ably de?ned norm||x
(out)
?x
(in)
||, vanishing at convergence,
that we will call in the following “scf norm”.
Ideally, the scf norm is a measure of the self-consistency
error on the total energy. Let us write an estimate of the latter
for the simplest case: an insulator with NC PPs and simple
LDA or GGA. At a given iteration we have
_
?
¯h
2
2m
?
2
+ V
ext
(r) + V
(in)
(r)
_
?
i
(r) =
i
?
i
(r), (2)
where
i
and ?
i
are Kohn-Sham energies and orbitals re-
spectively, i labels the occupied states, V
ext
is the sum of
the PPs of atomic cores (written for simplicity as a local po-
tential), the input Hartree and exchange-correlation potential
V
(in)
(r) = V
Hxc
[?
(in)
(r)] is a functional of the input charge
density ?
(in)
. The output charge density is given by
?
(out)
(r) =
i
|?
i
(r)|
2
. (3)
Let us compare the DFT energy calculated in the standard
way:
E =
i
_
?
?
i
(r)
_
?
¯h
2
2m
?
2
+ V
ext
(r)
_
?
i
(r)dr
+E
Hxc
[?
(out)
], (4)
where E
Hxc
is the Hartree and exchange-correlation en-
ergy, with the Harris-Weinert-Foulkes functional form, which
doesn’t use ?
(out)
:
E
=
i
_
?
?
i
(r)
_
?
¯h
2
2m
?
2
+ V
ext
(r) + V
(in)
(r)
_
?
i
(r)dr
?
_
?
(in)
V
(in)
(r) + E
Hxc
[?
(in)
] (5)
Both forms are variational, i.e. the ?rst-order variation of the
energy with respect to the charge density vanish, and both
converge to the same result when self-consistency is achieved.
Their difference can be approximated by the following expres-
sion, in which only the dominant Hartree term is considered:
E ?E
1
2
_
??(r)??(r
)
|r ?r
|
drdr
=
1
2
_
??(r)?V
H
(r
)dr (6)
where ?? = ?
(out)
??
(in)
and ?V
H
is the Hartree potential
generated by ??. Moreover it can be shown that, when ex-
change and correlation contributions to the electronic screen-
ing do not dominate over the electrostatic ones, this quantity
is an upper bound to the self-consistent error incurred when
using the standard form for the DFT energy. We therefore
take this term, which can be trivially calculated in reciprocal
space, as our squared scf norm:
||?
(out)
??
(in)
||
2
=
4?e
2
?
G
|??(G)|
2
G
2
, (7)
where G are the vectors in reciprocal space and ? is the vol-
ume of the unit cell.
Once the optimal linear combination of ?
(in)
from previous
iterations (typically 4 to 8) is determined, one adds a step in
the new search direction that is, in the simplest case, a fraction
of the optimal ?? or, taking advantage of some approximate
electronic screening[182], a preconditioned ??. In particular,
the simple, Thomas-Fermi, and local Thomas-Fermi mixing
described in Ref.[182] are implemented and used.
The above algorithm has been extended to more sophisti-
cated calculations, in which the x vector introduced above
may contain additional quantities: for DFT+U, occupancies of
13
atomic correlated states; for meta-GGA, kinetic energy den-
sity; for PAW, the quantities
i
?
i
|?
n
?
m
|?
i
, where the
? functions are the atomic-based projectors appearing in the
PAWformalism. The scf normis modi?ed accordingly in such
a way to include the additional variables in the estimated self-
consistency error.
B. Iterative diagonalization
During self-consistency one has to solve the generalized
eigenvalue problem for all N occupied states
H?
i
=
i
S?
i
, i = 1, . . . , N (8)
in which both H (the Hamiltonian) and S (the overlap ma-
trix) are available as operators (i.e. H? and S? products
can be calculated for a generic state ? ). Eigenvectors are
normalized according to the generalized orthonormality con-
straints ?
i
|S|?
j
= ?
ij
. This problem is solved using itera-
tive methods. Currently PWscf implements a block Davidson
algorithm and an alternative algorithm based on band-by-band
minimization using conjugate gradient.
1. Davidson
One starts from an initial set of orthonormalized trial or-
bitals ?
(0)
i
and of trial eigenvalues
(0)
i
= ?
(0)
i
|H|?
(0)
i
. The
starting set is typically obtained from the previous scf itera-
tion, if available, and if not, from the previous time step, or
optimization step, or from a superposition of atomic orbitals.
We introduce the residual vectors
g
(0)
i
= (H ?
(0)
i
S)?
(0)
i
, (9)
a measure of the error on the trial solution, and the correction
vectors ??
(0)
i
= Dg
(0)
i
, where D is a suitable approximation
to (H ?
(0)
i
S)
?1
. The eigenvalue problem is then solved in
the 2N-dimensional subspace spanned by the reduced basis
set ?
(0)
, formed by ?
(0)
i
= ?
(0)
i
and ?
(0)
i+N
= ??
(0)
i
:
2N
k=1
(H
jk
?
i
S
jk
)c
(i)
k
= 0, (10)
where
H
jk
= ?
(0)
j
|H|?
(0)
k
, S
jk
= ?
(0)
j
|S|?
(0)
k
. (11)
Conventional algorithms for matrix diagonalization are used
in this step. A new set of trial eigenvectors and eigenvalues is
obtained:
?
(1)
i
=
2N
j=1
c
(i)
j
?
(0)
j
,
(1)
i
= ?
(1)
i
|H|?
(1)
i
(12)
and the procedure is iterated until a satisfactory convergence is
achieved. Alternatively, one may enlarge the reduced basis set
with the new correction vectors ??
(1)
i
= Dg
(1)
i
, solve a 3N-
dimensional problem, and so on, until a pre?xed size of the
reduced basis set is reached. The latter approach is typically
slightly faster at the expenses of a larger memory usage.
The operator D must be easy to estimate. A natural choice
in the PWbasis set is a diagonal matrix, obtained keeping only
the diagonal term of the Hamiltonian:
k +G|D|k +G
=
?
GG
k +G|H ?S|k +G
(13)
where k is the Bloch vector of the electronic states under con-
sideration, |k+G
denotes PWs, an estimate of the highest
occupied eigenvalue. Since the Hamiltonian is a diagonally
dominant operator and the kinetic energy of PWs is the domi-
nant part at high G, this simple form is very effective.
2. Conjugate-Gradient
The eigenvalue problem of Eq.(8) can be recast into a se-
quence of constrained minimization problems:
min
_
_
?
i
|H|?
i
?
j?i
?
j
(?
i
|S|?
j
??
ij
)
_
_
, (14)
where the ?
j
are Lagrange multipliers. This can be solved
using preconditioned conjugate gradient algorithm with mi-
nor modi?cations to ensure constraint enforcement. The algo-
rithm here described was inspired by the conjugate-gradient
algorithm of Ref.[183], and is similar to one of the variants
described in Ref.[184].
Let us assume that eigenvectors ?
j
up to j = i ? 1 have
already been calculated. We start from an initial guess ?
(0)
for the i-th eigenvector, such that ?
(0)
|S|?
(0)
= 1 and
?
(0)
|S|?
j
= 0. We introduce a diagonal precondition ma-
trix P and auxiliary functions y = P
?1
? and solve the equiv-
alent problem
min
_
y|
¯
H|y ??
_
y|
¯
S|y ?1
__
, (15)
where
¯
H = PHP,
¯
S = PSP, under the additional orthonor-
mality constraints y|PS|?
j
= 0. The starting gradient of
Eq.(15)) is given by
g
(0)
= (
¯
H ??
¯
S)y
(0)
. (16)
By imposing that the gradient is orthonormal to the starting
vector: g
(0)
|
¯
S|y
(0)
= 0, one determines the value of the
Lagrange multiplier:
? =
y
(0)
|
¯
S
¯
H|y
(0)
y
(0)
|
¯
S
2
|y
(0)
. (17)
The remaining orthonormality constraints are imposed on
Pg
(0)
by explicit orthonormalization (e.g. Gram-Schmid) to
the ?
j
. We introduce the conjugate gradient h
(0)
, which for
14
the ?rst step is set equal to g
(0)
(after orthonormalization), and
the normalized direction n
(0)
= h
(0)
/h
(0)
|
¯
S|h
(0)
1/2
. We
search for the minimum of y
(1)
|
¯
H|y
(1)
along the direction
y
(1)
, de?ned as: [183]
y
(1)
= y
(0)
cos ? + n
(0)
sin?. (18)
This form ensures that the constraint on the norm is correctly
enforced. The calculation of the minimum can be analytically
performed and yields
? =
1
2
atan
_
a
(0)
(0)
?b
(0)
_
, (19)
where a
(0)
= 2y
(0)
|
¯
H|n
(0)
, b
(0)
= n
(0)
|
¯
H|n
(0)
, and
(0)
= y
(0)
|
¯
H|y
(0)
. The procedure is then iterated; at each
step the conjugate gradient is calculated from the gradient and
the conjugate gradient at the previous step, using the Polak-
Ribière formula:
h
= g
+ ?
(n?1)
h
(n?1)
, (20)
?
(n?1)
=
g
?g
(n?1)
|
¯
S|g
g
(n?1)
|
¯
S|g
(n?1)
. (21)
h
is subsequently re-orthogonalized to y
. We remarks
that in the practical implementation only Pg and Ph need to
be calculated and that only P
2
– the analogous of the D matrix
of Davidson algorithm – is actually used. A kinetic-only form
of P
2
has proved satisfactory:
k +G|P
2
|k +G
=
2m
¯h
2
(k +G)
2
?
GG
. (22)
C. Wavefunction extrapolation
In molecular dynamics runs and in structural relaxations,
extrapolations are employed to generate good initial guesses
for the wavefunctions at time t + dt from wavefunctions at
previous time steps. The extrapolation algorithms used are
similar to those described in Ref.[183]. The alignment pro-
cedure, needed when wavefunctions are the results of a self-
consistent calculation is as follows. The overlap matrix O
ij
between wavefunctions at consecutive time steps:
O
ij
= ?
i
(t + dt)|S(t + dt)|?
j
(t), (23)
can be used to generate the unitary transformation U [185]
that aligns ?(t +dt) to ?(t): ?
i
(t +dt) =
j
U
ij
?
j
(t +dt).
Since O is not unitary, it needs to be made unitary via e.g. the
unitarization procedure
U = (O
†
O)
?1/2
O
†
. (24)
The operation above is performed using a singular value de-
composition: Let the overlap matrix be O = vDw, where v
and w are unitary matrix and Dis a diagonal non-negative def-
inite matrix, whose eigenvalues are close to 1 if the two sets of
wavefunctions are very similar. The needed unitary transfor-
mation is then simply given by U w
†
v
†
. This procedure is
simpler than the original proposal and prevents the alignment
algorithm to break in the occasional situation where, due to
level crossing in the band structure between subsequent time
steps, one or more of the eigenvalues of the D matrix vanish.
D. Symmetry
Symmetry is exploited almost everywhere, with the notable
exception of CP. The latter is devised to study aperiodic sys-
tems or large supercells where symmetry is either absent or of
little use even if present.
In addition to lattice translations, the space group of a crys-
tal contains symmetry operations
ˆ
S combining rotations and
translations that leave the crystal unchanged:
ˆ
S ? {R|f },
where R is a 3 × 3 orthogonal matrix, f is a vector (called
fractional translation) and symmetry requires that any atomic
position, ?
s
is transformed into an equivalent one,
ˆ
S?
s
?
R(?
s
+ f ) = ?
ˆ
S(s)
. The rotational part of these operations
de?nes the crystal point group.
As a consequence of symmetry, roto-translated Kohn-Sham
orbitals are Kohn-Sham orbitals with the rotated Bloch vec-
tor:
ˆ
S?
i,k
(r) ? ?
i,k
(R
?1
r ?f ) = ?
i,Rk
(r). Where, strictly
speaking, the resulting wave-function at Rk does not neces-
sarily have the same band index as the original one but could
be some unitary transformation of states at Rk that share
with it the same single-particle eigenvalue. Since quantities
of physical interest are invariant for unitary rotations among
degenerate states this additional complication has no effect on
the ?nal result.
This is the basis for the symmetrization procedure used in
PWscf. One introduces a non-symmetrized charge density
(labeled by superscript
(ns)
) calculated on the irreducible BZ
(IBZ):
?
(ns)
(r) =
i
k?IBZ
w
k
|?
i,k
(r)|
2
. (25)
The factors w
k
(“weights”) are proportional to the number of
vectors in the star (i.e. inequivalent k vectors among all the
{Rk} vectors generated by the point-group rotations) and are
normalized to 1:
k?IBZ
w
k
= 1. Weights can either be
calculated or deduced from the literature on the special-point
technique[93, 94]. The charge density is then symmetrized as:
?(r) =
1
N
s
ˆ
S
ˆ
S?
(ns)
(r) =
1
N
s
ˆ
S
?
(ns)
(R
?1
r ?f ) (26)
where the sum runs over all N
s
symmetry operations.
The symmetrization technique can be extended to all quan-
tities that are expressed as sums over the BZ. Hellmann-
Feynman forces F
s
on atom s are thus calculated as follows:
F
s
=
1
N
s
ˆ
S
ˆ
SF
(ns)
s
=
1
N
s
ˆ
S
RF
(ns)
ˆ
S
?1
(s)
, (27)
15
where
ˆ
S
?1
(s) labels the atom into which the s?th atom trans-
forms (modulo a lattice translation vector) after application of
ˆ
S
?1
, the symmetry operation inverse of
ˆ
S. In a similar way
one determines the symmetrized stress, using the rule for ma-
trix transformation under a rotation:
?
??
=
1
N
s
ˆ
S
3
?,?=1
R
??
R
??
?
(ns)
??
. (28)
The PHonon package supplements the above technique
with a further strategy. Given the phonon wave-vector q, the
small group of q (the subgroup
ˆ
S
q
of crystal symmetry opera-
tions that leave q invariant) is identi?ed and the reducible rep-
resentation de?ned by the 3N
at
atomic displacements along
cartesian axis is decomposed into n
irr
irreducible representa-
tions (irreps) ?
(q)
j
, j = 1, . . . , n
irr
. The dimensions of the ir-
reducible representations are small, with ?
j
? 3 in most cases,
up to 6 in some special cases (zone-boundary wave-vectors q
in nonsymmorphic groups). Each irrep, j, is therefore de?ned
by a set of ?
j
linear combinations of atomic displacements
that transform into each other under the symmetry operations
of the small group of q. In the self-consistent solution of the
linear response equations, only perturbations associated to a
given irrep need to be treated together and different irreps can
be solved independently. This feature is exploited to reduce
the amount of memory required by the calculation and is suit-
able for coarse-grained parallelization and for execution on a
GRID infrastructure [186].
The wavefunction response, ??
(j,?)
k+q,i
(r), to displacements
along irrep j, ?
(q)
j,?
(where ? = 1, . . . , ?
j
labels different part-
ners of the given irrep), is then calculated. The lattice-periodic
unsymmetrized charge response, ??
(ns)
q,j,?
(r), has the form:
??
(ns)
q,j,?
(r) = e
?iq·r
4
i
k?IBZ(q)
w
k
?
?
k,i
(r)??
(j,?)
k+q,i
(r),
(29)
where the notation IBZ(q) indicates the IBZ calculated as-
suming the small group of q as symmetry group, and the
weights w
k
are calculated accordingly. The symmetrized
charge response is calculated as
??
q,j,?
(r) =
1
N
s
(q)
ˆ
S
q
e
?iqf
?
j
?=1
D(
ˆ
S
q
)
??
??
(ns)
q,j,?
(R
?1
r?f )
(30)
where D(
ˆ
S
q
) is the matrix representation of the action of the
symmetry operation
ˆ
S
q
? {R|f } for the j?th irrep ?
(q)
j
. At
the end of the self-consistent procedure, the force constant
matrix C
s?,t?
(q) (where s, t label atoms, ?, ? cartesian co-
ordinates) is calculated. Force constants at all vectors in the
star of q are then obtained using symmetry:
C
s?,t?
(Rq) =
?,?
R
??
R
??
C
ˆ
S
?1
(s)?,
ˆ
S
?1
(t)?
(q), (31)
where
ˆ
S ? {R|f } is a symmetry operation of the crystal
group but not of the small group of q.
E. Fock exchange
Hybrid functionals are characterized by the inclusion of a
fraction of exact (i.e. non-local) Fock exchange in the de?-
nition of the exchange-correlation functional. For a periodic
system, the Fock exchange energy per unit cell is given by:
E
x
= ?
e
2
N
kv
k
v
_
?
?
kv
(r)?
k
v
(r)?
?
k
v
(r
)?
kv
(r
)
|r ?r
|
drdr
,
(32)
where an insulating and non magnetic system is assumed for
simplicity. Integrals and wave-function normalizations are de-
?ned over the whole crystal volume, V = N? (? being the
unit cell volume), and the summations run over all occupied
bands and all N k-points de?ned in the BZ by Born-von Kár-
mán boundary conditions. The calculation of this term is per-
formed exploiting the dual-space formalism: auxiliary coden-
sities, ?k
,v
k,v
(r) = ?
?
k
,v
(r)?
k,v
(r) are computed in real space
and transformed to reciprocal space by FFT, where the associ-
ated electrostatic energies are accumulated. The application of
the Fock exchange operator to a wavefunction involves addi-
tional FFTs and real-space array multiplications. These basic
operations need to be repeated for all the occupied bands and
all the points in the BZ grid. For this reason the computational
cost of the exact exchange calculation is very high, at least an
order of magnitude larger than for non-hybrid functional cal-
culations.
In order to limit the computational cost, an auxiliary grid of
q-points in the BZ, centered at the ? point, can be introduced
and the summation over k
be limited to the subset k
= k+q.
Of course convergence with respect to this additional param-
eter needs to be checked, but often a grid coarser than the one
used for computing densities and potentials is suf?cient.
The direct evaluation of the Fock energy on regular grids in
the BZ is however problematic due to an integrable divergence
that appears in the q ?0 limit. This problem is addressed re-
sorting to a procedure, ?rst proposed by Gygi and Baldereschi
[187], where an integrable term that displays the same diver-
gence is subtracted from the expression for the exchange en-
ergy and its analytic integral over the BZ is separately added
back to it. Some care must still be paid [175] in order to es-
timate the contribution of the q = 0 term in the sum, which
contains a 0/0 limit that cannot be calculated from informa-
tion at q = 0 only. This term is estimated [175] assuming that
the grid of q-points used for evaluating the exchange integrals
is dense enough that a coarser grid, including only every sec-
ond point in each direction, would also be equally accurate.
Since the limiting term contributes to the integral with differ-
ent weights in the two grids, one can extract its value from the
condition that the two integral give the same result. This pro-
cedure removes an error proportional to the inverse of the unit
cell volume ? that would otherwise appear if this term were
simply neglected.
16
[1] N. Marzari, MRS BULLETIN 31, 681 (2006).
[2] R. G. Parr and W. Yang, Density Functional Theory of Atoms
and Molecules (Oxford University Press, 1989).
[3] R. M. Dreizler and E. K. U. Gross, Density Functional Theory
(Springer-Verlag, 1990).
[4] R. Martin, Electronic Structure: Basic Theory and Practi-
cal Methods (Cambridge University Press, Cambridge, UK,
2004).
[5] E. J. Baerends, J. Autschbach, A. Bérces, F. Bickelhaupt,
C. Bo, P. M. Boerrigter, L. Cavallo, D. P. Chong, L. Deng,
R. M. Dickson, et al., ADF2008.01, SCM, Theoretical Chem-
istry, Vrije Universiteit, Amsterdam, The Netherlands, URL
http://www.scm.com.
[6] R. Dovesi, B. Civalleri, R. Orlando, C. Roetti, and V. R. Saun-
ders, Ab initio quantum simulation in solid state chemistry, in
Reviews in Computational Chemistry, Ch. 1, Vol. 21, K. B.
Lipkowitz, R. Larter, T. R. Cundari editors, John Wiley and
Sons, New York (2005), URL http://www.crystal.
unito.it.
[7] G. Karlström, R. Lindh, P.-A. Malmqvist, B. O. Roos,
U. Ryde, V. Veryazov, P.-O. Widmark, M. Cossi, B. Schim-
melpfennig, P. Neogrady, et al., Comp. Mater. Sci. 28, 222
(2003), URL http://www.teokem.lu.se/molcas.
[8] R. Ahlrichs, F. Furche, C. Hättig, W. M. Klopper, M. Sierka,
and F. Weigend, TURBOMOLE, URL http://www.
turbomole-gmbh.com.
[9] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuse-
ria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr.,
T. Vreven, K. N. Kudin, J. C. Burant, et al., Gaussian 03,
Revision C.02, Gaussian, Inc., Wallingford, CT, 2004, URL
http://www.gaussian.com.
[10] Materials Studio, URL http://accelrys.com/
products/materials-studio.
[11] Jaguar: rapid ab initio electronic structure package,
Schrödinger, LLC., URL http://www.schrodinger.
com/.
[12] Spartan, Wavefunction, Inc., URL http://www.
wavefun.com/products/spartan.html.
[13] G. Kresse and J. Furthmuller, Comp. Mater. Sci. 6, 15 (1996),
URL http://cmp.univie.ac.at/vasp.
[14] URL http://www.castep.org.
[15] The CPMD consortium page, coordinated by M. Parrinello
and W. Andreoni, Copyright IBM Corp 1990-2008, Copy-
right MPI für Festkörperforschung Stuttgart 1997-2001, URL
http://www.cpmd.org.
[16] E. J. Bylaska, W. A. de Jong, N. Govind, K. Kowalski,
T. P. Straatsma, M. Valiev, D. Wang, E. Apra, T. L. Win-
dus, J. Hammond, et al., Nwchem, a computational chem-
istry package for parallel computers, version 5.1, Paci?c
Northwest National Laboratory, Richland, Washington 99352-
0999, USA, URL http://www.emsl.pnl.gov/docs/
nwchem/nwchem.html.
[17] R. A. Kendall, E. Apra, D. E. Bernholdt, E. J. Bylaska,
M. Dupuis, G. I. Fann, R. J. Harrison, J. Ju, J. A. Nichols,
J. Nieplocha, et al., Comput. Phys. Commun. 128, 260 (2000).
[18] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert,
M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A.
Nguyen, S. J. Su, et al., J. Comput. Chem. 14, 1347 (1993).
[19] M. S. Gordon and M. W. Schmidt, Advances in elec-
tronic structure theory: Gamess a decade later, in The-
ory and Applications of Computational Chemistry, the ?rst
forty years, C. E. Dykstra, G. Frenking, K. S. Kim, G.
E. Scuseria editors, ch. 41, pp 1167-1189, Elsevier, Ams-
terdam, 2005, URL http://www.msg.chem.iastate.
edu/GAMESS/GAMESS.html.
[20] URL http://www.gnu.org/licenses/.
[21] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs,
G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jol-
let, et al., Comp. Mater. Sci. 25, 478 (2002), URL http:
//www.abinit.org.
[22] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello,
T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103
(2005), URL http://cp2k.berlios.de.
[23] S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng. 4, 56
(2002), URL http://www.fysik.dtu.dk/campos.
[24] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys.
Rev. B 71 (2005), URL https://wiki.fysik.dtu.
dk/gpaw.
[25] T. D. Crawford, C. D. Sherrill, E. F. Valeev, J. T. Fermann,
R. A. King, M. L. Leininger, S. T. Brown, C. L. Janssen, E. T.
Seidl, J. P. Kenny, et al., J. Comput. Chem. 28, 1610 (2007),
URL http://www.psicode.org.
[26] A. A. Mosto?, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt,
and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
[27] URL http://qe-forge.org.
[28] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel,
J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling,
A. McKenney, et al., LAPACK Users’ Guide (Society for In-
dustrial and Applied Mathematics, Philadelphia, PA, 1999),
3rd ed., ISBN 0-89871-447-8 (paperback).
[29] M. Frigo and S. G. Johnson, Proceedings of the IEEE 93, 216
(2005), special issue on "Program Generation, Optimization,
and Platform Adaptation".
[30] M. P. I. Forum, Int. J. Supercomputer Applications 8 (3/4)
(1994).
[31] URL http://en.wikipedia.org/wiki/
SourceForge.
[32] W. E. Pickett, Comput. Phys. Rep, 9, 115 (1989).
[33] I. Dabo, B. Kozinsky, N. E. Singh-Miller, and N. Marzari,
Phys. Rev. B 77, 115139 (2008).
[34] L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425
(1982).
[35] D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett.
43, 1494 (1979).
[36] D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
[37] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
[38] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, D. J.
Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992).
[39] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys.
Rev. Lett. 91, 146401 (2003).
[40] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105,
9982 (1996).
[41] A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
[42] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J.
Frisch, J. Phys. Chem. 98, 11623 (1994).
[43] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[44] H. Hellmann, Einführung in die Quantenchemie (Deutliche,
1937).
[45] R. P. Feynman, Phys. Rev. 56, 340 (1939).
[46] O. H. Nielsen and R. M. Martin, Phys. Rev. B 32, 3780 (1985).
[47] P. Pyykkö, Chem. Rev. 88, 563 (1988).
[48] T. Oda, A. Pasquarello, and R. Car, Phys. Rev. Lett. 80, 3622
17
(1998).
[49] R. Gebauer and S. Baroni, Phys. Rev. B 61, R6459 (2000).
[50] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
[51] R. M. Wentzcovitch and J. L. Martins, Sol. St. Commun. 78,
831 (1991).
[52] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196
(1980).
[53] R. M. Wentzcovitch, Phys. Rev. B 44, 2358 (1991).
[54] S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58,
1861 (1987).
[55] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi,
Rev. Mod. Phys. 73, 515 (2001).
[56] X. Gonze, Phys. Rev. A 52, 1096 (1995).
[57] G. Henkelman and H. Jónsson, J. Chem. Phys. 111, 7010
(1999).
[58] W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66,
052301 (2002).
[59] K. J. Caspersen and E. A. Carter, P. Natl. Acad. Sci. USA 102,
6738 (2005).
[60] H. J. Choi and J. Ihm, Phys. Rev. B 59, 2267 (1999).
[61] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
[62] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65,
035109 (2001).
[63] C. J. Pickard and F. Mauri, Phys. Rev. B 63, 245101 (2001).
[64] C. J. Pickard and F. Mauri, Phys. Rev. Lett. 91, 196401 (2003).
[65] M. Taillefumier, D. Cabaret, A. M. Flank, and F. Mauri, Phys.
Rev. B 66, 195107 (2002).
[66] S. Scandolo, P. Giannozzi, C. Cavazzoni, S. de Gironcoli,
A. Pasquarello, and S. Baroni, Z. Kristallogr. 220, 574 (2005).
[67] J. Tersoff and D. R. Hamann, Phys. Rev. B 31, 805 (1985).
[68] A. D. Becke and K. E. Edgecombe, J. Chem. Phys. 92, 5397
(1990).
[69] A. Szabo and N. Ostlund, Modern Quantum Chemistry
(Dover, 1996).
[70] A. Baldereschi, S. Baroni, and R. Resta, Phys. Rev. Lett. 61,
734 (1988).
[71] T. H. Group, URL http://www.hdfgroup.org.
[72] URL http://www.unidata.ucar.edu/software/
netcdf.
[73] URL http://www.quantum-simulation.org.
[74] G. Bussi, URL http://www.s3.infm.it/iotk.
[75] M. Fuchs and M. Schef?er, Comp. Phys. Comm. 119, 67
(1999).
[76] URL http://www.physics.rutgers.edu/~dhv/
uspp/.
[77] OPIUM – pseudopotential generation project, URL http:
//opium.sourceforge.net.
[78] P. Giannozzi and C. Cavazzoni, Nuovo Cimento C (in press,
2009).
[79] J. Hutter and A. Curioni, ChemPhysChem 6, 1788 (2005).
[80] L. S. Blackford, J. Choi, A. Cleary, J. Demmel, I. Dhillon,
J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley,
et al., SC Conference 0, 5 (1996).
[81] F. Gygi, E. W. Draeger, M. Schulz, B. R. de Supin-
ski, J. A. Gunnels, V. Austel, J. C. Sexton, F. Franchetti,
S. Kral, C. W. Ueberhuber, et al., IBM J. Res. Dev.
52, 137 (2008), URL http://eslab.ucdavis.edu/
software/qbox/index.htm.
[82] P. Giannozzi, F. de Angelis, and R. Car, J. Chem. Phys. 120,
5903 (2004).
[83] A. Dal Corso, A pseudopotential plane waves program
(pwscf) and some case studies, in Lecture Notes in Chem-
istry, Vol. 67, C. Pisani editor, Springer Verlag, Berlin (1996).
[84] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
[85] K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vander-
bilt, Phys. Rev. B 47, 10142 (1993).
[86] J. A. White and D. M. Bird, Phys. Rev. B 50, 4954 (1994).
[87] A. Dal Corso and A. Mosca Conte, Phys. Rev. B 71, 115106
(2005).
[88] A. Mosca Conte, Ph.D. thesis, SISSA/ISAS, Trieste Italy
(2007), URL http://www.sissa.it/cm/thesis/
2007/moscaconte.pdf.
[89] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B
44, 943 (1991).
[90] M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105
(2005).
[91] H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100,
126404 (2008).
[92] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
[93] D. J. Chadi and M. L. Cohen, Phys. Rev. B 8, 5747 (1973).
[94] H. J. Monkhorst and J. Pack, Phys. Rev. B 13, 5188 (1976).
[95] M. Methfessel and A. Paxton, Phys. Rev. B 40, 3616 (1989).
[96] N. Marzari, D. Vanderbilt, A. de Vita, and M. C. Payne, Phys.
Rev. Lett. 82, 3296 (1999).
[97] P. E. Blöchl, O. Jepsen, and O. Andersen, Phys. Rev. B 49,
16223 (1994).
[98] N. D. Mermin, Phys. Rev. 137, A1441 (1965).
[99] R. Fletcher, Practical Methods of Optimization (John Wiley
and Sons, 1987).
[100] S. R. Billeter, A. J. Turner, and W. Thiel, Phys. Chem. Chem.
Phys. 2, 2177 (2000).
[101] S. R. Billeter, A. Curioni, and W. Andreoni, Comp. Mater. Sci.
27, 437 (2003).
[102] C. Micheletti, A. Laio, and M. Parrinello, Phys. Rev. Lett. 92,
170601 (2004).
[103] L. Verlet, Phys. Rev. 159, 98 (1967).
[104] M. P. Allen and D. J. Tildesley, Computer Simulations of Liq-
uids (Clarendon Press, 1986).
[105] M. Bernasconi, G. L. Chiarotti, P. Focher, S. Scandolo,
E. Tosatti, and M. Parrinello, J. Phys. Chem. Solids 56, 501
(1995).
[106] P. Umari and A. Pasquarello, Phys. Rev. Lett. 89, 157602
(2002).
[107] I. Souza, J. Íñiguez, and D. Vanderbilt, Phys. Rev. Lett. 89,
117602 (2002).
[108] K. Kunc and R. Resta, Phys. Rev. Lett. 51, 686 (1983).
[109] J. Tobik and A. Dal Corso, J. Chem. Phys. 120, 9934 (2004).
[110] D. A. Scherlis, J.-L. Fattebert, F. Gygi, M. Cococcioni, and
N. Marzari, J. Chem. Phys. 124, 074103 (2006).
[111] I. Dabo, E. Cances, Y. Li, and N. Marzari (2009), URL http:
//arxiv.org/abs/0901.0096v2.
[112] C. Cavazzoni and G. L. Chiarotti, Comput. Phys. Commun.
123, 56 (1999).
[113] A. Pasquarello, K. Laasonen, R. Car, C. Lee, and D. Vander-
bilt, Phys. Rev. Lett. 69, 1982 (1992).
[114] M. d’Avezac, M. Calandra, and F. Mauri, Phys. Rev. B 71,
205210 (2005).
[115] H.-L. Sit, M. Cococcioni, and N. Marzari, J. of Electroanalyt-
ical Chemistry 607, 107 (2007).
[116] H.-L. Sit, M. Cococcioni, and N. Marzari, Phys. Rev. Lett. 97,
028303 (2006).
[117] N. Marzari, D. Vanderbilt, and M. C. Payne, Phys. Rev. Lett.
79, 1337 (1997).
[118] F. Tassone, F. Mauri, and R. Car, Phys. Rev. B 50, 10561
(1994).
[119] D. J. Tobias, G. J. Martyna, and M. L. Klein, J. Phys. Chem.
97, 12959 (1993).
[120] G. J. Martyna, M. L. Klein, and M. Tuckerman, The Journal
18
of Chemical Physics 97, 2635 (1992).
[121] M. C. Payne, M. P. Teter, D. C. Allen, T. A. Arias, and J. D.
Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992).
[122] D. Marx and J. Hutter, Modern methods and algorithms of
quantum chemistry, John von Neumann Institute for Comput-
ing, FZ Jülich, pp. 301–449 (2000).
[123] V. Dubois, P. Umari, and A. Pasquarello, Chem. Phys. Lett.
390, 193 (2004).
[124] F. Giustino and A. Pasquarello, Phys. Rev. Lett. 95, 187402
(2005).
[125] P. Umari and A. Pasquarello, Diam. Relat. Mater. 14, 1255
(2005).
[126] P. Umari and A. Pasquarello, Phys. Rev. Lett. 95, 137401
(2005).
[127] L. Giacomazzi, P. Umari, and A. Pasquarello, Phys. Rev. Lett.
95, 075505 (2005).
[128] P. Umari and A. Pasquarello, Phys. Rev. Lett. 98, 176402
(2007).
[129] P. Giannozzi and S. Baroni, J. Chem. Phys. 100, 8537 (1994).
[130] X. Gonze, Phys. Rev. B 55, 10337 (1997).
[131] X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
[132] A. Dal Corso, A. Pasquarello, and A. Baldereschi, Phys. Rev.
B 56, R11369 (1997).
[133] A. Dal Corso, Phys. Rev. B 64, 235118 (2001).
[134] F. Favot and A. Dal Corso, Phys. Rev. B 60, 11427 (1999).
[135] A. Dal Corso and S. de Gironcoli, Phys. Rev. B62, 273 (2000).
[136] A. Dal Corso, submitted (2009).
[137] A. Dal Corso, Phys. Rev. B 76, 054308 (2007).
[138] F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B 76,
165108 (2007).
[139] M. Lazzeri and S. de Gironcoli, Phys. Rev. B 65, 245402
(2002).
[140] M. Lazzeri and F. Mauri, Phys. Rev. Lett. 90, 036401 (2003).
[141] M. Lazzeri and F. Mauri, Phys. Rev. B 68, 161101(R) (2003).
[142] D. D. Koelling and B. N. Harmon, J. Phys. C: Solid State Phys.
10, 3107 (1977).
[143] A. H. MacDonald and S. H. Vosko, J. Phys. C: Solid State
Phys. 12, 2977 (1979).
[144] A. K. Rajagopal and J. Callaway, Phys. Rev. B 7, 1912 (1973).
[145] N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
[146] A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos,
Phys. Rev. B 41, 1227 (1990).
[147] G. Kresse and J. Hafner, J. Phys.-Condens. Mat. 6, 8245
(1994).
[148] A. Smogunov, A. Dal Corso, and E. Tosatti, Phys. Rev. B 70,
045417 (2004).
[149] A. Dal Corso, A. Smogunov, and E. Tosatti, Phys. Rev. B 74,
045429 (2006).
[150] M. Profeta, F. Mauri, and C. J. Pickard, J. Am. Chem. Soc.
125, 541 (2003).
[151] C. G. van de Walle and P. E. Blöchl, Phys. Rev. B 47, 4244
(1993).
[152] F. Mauri and S. G. Louie, Phys. Rev. Lett. 76, 4246 (1996).
[153] F. Mauri, B. Pfrommer, and S. G. Louie, Phys. Rev. Lett. 77,
5300 (1996).
[154] C. J. Pickard and F. Mauri, Phys. Rev. Lett. 88, 086403 (2002).
[155] T. Thonhauser, D. Ceresoli, A. Mosto?, N. Marzari, R. Resta,
and D. Vanderbilt, submitted to Phys. Rev. Lett. (2007), URL
http://arxiv.org/abs/0709.4429v1.
[156] C. Gougoussis, M. Calandra, A. Seitsonen, C. Brouder,
A. Shukla, and F. Mauri (2008), arXiv/0806.4706v1.
[157] URL http://www.wannier.org/.
[158] P. Silvestrelli, N. Marzari, D. Vanderbilt, and M. Parrinello,
Solid State Commun. 107, 7 (1998).
[159] A. Calzolari, I. Souza, N. Marzari, and M. B. Nardelli, Phys.
Rev. B 69, 035108 (2004).
[160] Y.-S. Leee, M. B. Nardelli, and N. Marzari, Phys. Rev. Lett.
95, 076804 (2005).
[161] J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, Phys. Rev.
B 75, 195121 (2007).
[162] A. Kokalj, J. Mol. Graph. Model. 17, 176 (1999), URL http:
//www.xcrysden.org/.
[163] W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph.
Model. 14, 33 (1996).
[164] URL http://www.tcm.phy.cam.ac.uk/~mdt26/
casino2.html.
[165] D. Korotin, A. V. Kozhevnikov, S. L. Skornyakov, I. Leonov,
N. Binggeli, V. I. Anisimov, and G. Trimarchi, Euro. Phys. J.
B 65, 91 (2008), URL http://[email protected].
[166] L. Martin-Samos and G. Bussi, J. Comp. Phys. Comm (2009),
URL http://sax-project.org.
[167] URL http://dp-code.org.
[168] URL http://www.bethe-salpeter.org.
[169] A. Kokalj, Comp. Mater. Sci. 28, 155 (2003), URL http:
//www-k3.ijs.si/kokalj/guib/.
[170] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
[171] M. A. L. Marques, C. L. Ullrich, F. Nogueira, A. Rubio,
K. Burke, and E. K. U. Gross, eds., Time-Dependent Den-
sity Functional Theory, vol. 706 of Lecture notes in Physics
(Springer-Verlag, Berlin, Heidelberg, 2006), DOI-10.1007/3-
540-35426-3-17.
[172] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601
(2002).
[173] D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, J. Chem. Phys.
128, 154105 (2008).
[174] P. Umari, G. Stenuit, and S. Baroni (2008), arXiv/0811.1453.
[175] H.-V. Nguyen and S. de Gironcoli, submitted to Phys. Rev. B
(2009), URL http://arxiv.org/abs/0902.0889.
[176] H.-V. Nguyen and S. de Gironcoli, submitted to Phys. Rev. B
(2009), URL http://arxiv.org/abs/0902.0883.
[177] H.-V. Nguyen, Ph.D. thesis, SISSA (2008), URL
http://www.sissa.it/cm/thesis/2008/
VietHuyNguyen_PhDthesis.pd%f.
[178] A. Marini, C. Hogan, M. Grüning, and D. Varsano, Comp.
Phys. Commun. p. in press (2009), URL http://www.
yambo-code.org.
[179] A. Calzolari, I. Souza, N. Marzari, and M. B. Nardelli,
Phys. Rev. B 69, 035108 (2004), URL http://www.
wannier-transport.org.
[180] URL http://[email protected].
[181] A. Kokalj, pwtk: a Tcl scripting interface to PWscf, soon to
be available on http://qe-forge.net.
[182] D. Raczkowski, A. Canning, and L. W. Wang, Phys. Rev. B
64, R121101 (2001).
[183] T. A. Arias, M. C. Payne, and J. D. Joannopoulos, Phys. Rev.
B 45, 1538 (1992).
[184] A. Qteish, Phys. Rev. B 52, 14497 (1995).
[185] C. A. Mead, Rev. Mod. Phys. 64, 51 (1992).
[186] R. di Meo, A. Dal Corso, P. Giannozzi, and S. Cozzini, Pro-
ceedings of the COST School, Trieste (2009), to be published.
[187] F. Gygi and A. Baldereschi, Phys. Rev. B 34, 4405 (1986).
doc_420872609.pdf
Quantum technology is a new field of physics and engineering, which transitions some of the stranger features of quantum mechanics, especially quantum entanglement and most recently quantum tunneling, into practical applications such as quantum computing, quantum cryptography, quantum simulation, quantum metrology, quantum sensing, and quantum imaging.
QUANTUM ESPRESSO: a modular and open-source software project
for quantum simulations of materials
Paolo Giannozzi,
1, 2
Stefano Baroni,
1, 3
Nicola Bonini,
4
Matteo Calandra,
5
Roberto Car,
6
Carlo Cavazzoni,
7, 8
Davide Ceresoli,
4
Guido L. Chiarotti,
9
Matteo Cococcioni,
10
Ismaila Dabo,
11
Andrea Dal Corso,
1, 3
Stefano de Gironcoli,
1, 3
Ralph Gebauer,
1, 12
Uwe Gerstmann,
13
Christos Gougoussis,
5
Anton Kokalj,
14
Michele Lazzeri,
5
Layla Martin Samos Colomer,
1
Nicola
Marzari,
4
Francesco Mauri,
5
Stefano Paolini,
3, 9
Alfredo Pasquarello,
15
Lorenzo Paulatto,
1, 3
Carlo Sbraccia
?
,
1
Sandro
Scandolo,
1, 12
Gabriele Sclauzero,
1, 3
Ari P. Seitsonen,
5
Alexander Smogunov,
12
Paolo Umari,
1
and Renata M. Wentzcovitch
10
1
CNR-INFM Democritos National Simulation Center, 34100 Trieste, Italy
2
Dipartimento di Fisica, Università degli Studi di Udine, via delle Scienze 208, 33100 Udine, Italy
3
SISSA – Scuola Internazionale Superiore di Studi Avanzati, via Beirut 2-4, 34014 Trieste Grignano, Italy
4
Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139 USA
5
Institut de Minéralogie et de Physique des Milieux Condensés,
Université Pierre et Marie Curie, CNRS, IPGP, 140 rue de Lourmel, 75015 Paris, France
6
Department of Chemistry, Princeton University, Princeton NJ 08544, USA
7
CINECA National Supercomputing Center, Casalecchio di Reno, 40033 Bologna, Italy
8
CNR-INFM S3 Research Center, 41100 Modena, Italy
9
SPIN s.r.l. via del Follatoio 12, 34148 Trieste, Italy
10
Department of Chemical Engineering and Materials Science, University of Minnesota,
151 Amundson Hall, 421 Washington Avenue SE, Minneapolis MN 55455, USA
11
Université Paris-Est, CERMICS, Projet Micmac ENPC-INRIA,
6-8 avenue Blaise Pascal, 77455 Marne-la-Vallée Cedex 2, France
12
The Abdus Salam International Centre for Theoretical Physics, Strada Costiera 11, 34014 Trieste Grignano, Italy
13
Theoretische Physik, Universität Paderborn, D-33098 Paderborn, Germany
14
Jožef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia
15
Ecole Polytechnique Fédérale de Lausanne (EPFL), Institute of Theoretical Physics,
and Institut Romand de Recherche Numérique en Physique des Matériaux (IRRMA), CH-1015 Lausanne, Switzerland
QUANTUM ESPRESSO is an integrated suite of computer codes for electronic-structure calculations and
materials modeling, based on density-functional theory, plane waves, and pseudopotentials (norm-conserving,
ultrasoft, and projector-augmented wave). QUANTUM ESPRESSO stands for opEn Source Package for Re-
search in Electronic Structure, Simulation, and Optimization. It is freely available to researchers around the
world under the terms of the GNU General Public License. QUANTUM ESPRESSO builds upon newly-
restructured electronic-structure codes that have been developed and tested by some of the original authors
of novel electronic-structure algorithms and applied in the last twenty years by some of the leading materials
modeling groups worldwide. Innovation and ef?ciency are still its main focus, with special attention paid to
massively-parallel architectures, and a great effort being devoted to user friendliness. QUANTUM ESPRESSO
is evolving towards a distribution of independent and inter-operable codes in the spirit of an open-source project,
where researchers active in the ?eld of electronic-structure calculations are encouraged to participate in the
project by contributing their own codes or by implementing their own ideas into existing codes.
I. INTRODUCTION
The combination of methodological and algorithmic in-
novations and ever-increasing computer power is deliver-
ing a simulation revolution in materials modeling, starting
from the nanoscale up to bulk materials and devices [1].
Electronic-structure simulations based on density-functional
theory (DFT) [2–4] have been instrumental to this revolu-
tion, and their application has now spread outside a restricted
core of researchers in condensed-matter theory and quantum
chemistry, involving a vast community of end users with very
diverse scienti?c backgrounds and research interests. Sus-
taining this revolution and extending its bene?cial effects to
the many ?elds of science and technology that can capital-
?
Present address: Constellation Energy Commodities Group 7th Floor, 61
Aldwich, London, WC2B 4AE, United Kingdom
ize on it represents a multifold challenge. In our view it
is also a most urgent, fascinating and fruitful endeavor, able
to deliver new forms for scienti?c exploration and discovery,
where a very complex infrastructure—made of software rather
than hardware—can be made available to any researcher, and
whose capabilities continue to increase thanks to the method-
ological innovations and computing power scalability alluded
to above.
Over the past few decades, innovation in materials simula-
tion and modeling has resulted from the concerted efforts of
many individuals and groups worldwide, often of small size.
Their success has been made possible by a combination of
competences, ranging from the ability to address meaningful
and challenging problems, to a rigorous insight into theoret-
ical methods, ending with a marked sensibility to matters of
numerical accuracy and algorithmic ef?ciency. The readiness
to implement new algorithms that utilize novel ideas requires
total control over the software being used–for this reason,
the physics community has long relied on in-house computer
2
codes to develop and implement new ideas and algorithms.
Transitioning these development codes to production tools is
nevertheless essential, both to extensively validate new meth-
ods and to speed up their acceptance by the scienti?c com-
munity. At the same time, the dissemination of codes has to
be substantial, to justify the learning efforts of PhD students
and young postdocs who would soon be confronted with the
necessity of deploying their competences in different research
groups. In order to sustain innovation in numerical simulation,
we believe there should be little, if any, distinction between
development and production codes; computer codes should
be easy to maintain, to understand by different generations of
young researchers, to modify, and extend; they should be easy
to use by the layman, as well as general and ?exible enough
to be enticing for a vast and diverse community of end users.
One easily understands that such con?icting requirements can
only be tempered, if anything, within organized and modular
software projects.
Software modularity also comes as a necessity when com-
plex problems in complex materials need to be tackled with
an array of different methods and techniques. Multiscale ap-
proaches, in particular, strive to combine methods with differ-
ent accuracy and scope to describe different parts of a com-
plex system, or phenomena occurring at different time and/or
length scales. Such approaches will require software packages
that can perform different kinds of computations on different
aspects of the same problem and/or different portions of the
same system, and that allow for interoperability or joint us-
age of the different modules. Different packages should at the
very least share the same input/output data formats; ideally
they should also share a number of mathematical and appli-
cation libraries, as well as the internal representation of some
of the data structures on which they operate. Individual re-
searchers or research groups ?nd it increasingly dif?cult to
meet all these requirements and to continue to develop and
maintain in-house software project of increasing complexity.
Thus, different and possibly collaborative solutions should be
sought.
A successful example comes from the software for sim-
ulations in quantum chemistry, that has often (but not al-
ways) evolved towards commercialization: the development
and maintenance of most well-known packages is devolved to
non-pro?t [5–8] or commercial companies [9–12]. The soft-
ware is released for a fee under some proprietary license that
may impose several restrictions to the availability of sources
(computer code in a high-level language) and to what can be
done with the software. This model has worked well, and is
also used by some of the leading development groups in the
condensed-matter electronic-structure community [13, 14],
while some proprietary projects allow for some free academic
usage of their products [14–19]. A commercial endeavor also
brings the distinctive advantage of a professional approach to
software development, maintenance, documentation, and sup-
port.
We believe however that a more interesting and fruitful al-
ternative can be pursued, and one that is closer to the spirit
of science and scienti?c endeavor, modeled on the experience
of the open-source software community. Under this model, a
large community of users has full access to the source code
and the development material, under the coordination of a
smaller group of core developers. In the long term, and in
the absence of entrenched monopolies, this strategy could be
more effective in providing good software solutions and in
nurturing a community engaged in providing those solutions,
as compared to the proprietary software strategy. In the case
of software for scienti?c usage, such an approach has the ad-
ditional, and by no means minor, advantage to be in line with
the tradition and best practice of science, that require repro-
ducibility of other people’s results, and where collaboration is
no less important than competition.
In this paper we will shortly describe our answer to the
above-mentioned problems, as embodied in our QUANTUM
ESPRESSO project (indeed, QUANTUM ESPRESSO stands
for opEn Source Package for Research in Electronic Structure,
Simulation, and Optimization). First, in Sec. 2, we describe
the guiding lines of our effort. In Sec. 3, we give an overview
of the current capabilities of QUANTUM ESPRESSO. In Sec.
4, we provide a short description of each software component
presently distributed within QUANTUM ESPRESSO. Finally,
Sec. 5 describes current developments and offers a perspec-
tive outlook. The Appendix sections discuss some of the more
speci?c technical details of the algorithms used, that have not
been documented elsewhere.
II. THE QUANTUM ESPRESSO PROJECT
QUANTUM ESPRESSO is an integrated suite of computer
codes for electronic-structure calculations and materials mod-
eling based on density-functional theory, plane waves basis
sets and pseudopotentials to represent electron-ion interac-
tions. QUANTUM ESPRESSO is free, open-source, software
distributed under the terms of the GNU General Public Li-
cense (GPL) [20].
The two main goals of this project are to foster method-
ological innovation in the ?eld of electronic-structure sim-
ulations and to provide a wide and diverse community of
end users with highly ef?cient, robust, and user-friendly soft-
ware implementing the most recent innovations in this ?eld.
Other open-source projects [21–25] exist, besides QUAN-
TUM ESPRESSO, that address electronic-structure calcula-
tions and various materials simulation techniques based on
them. Unlike some of these projects, QUANTUM ESPRESSO
does not aim at providing a single monolithic code able to per-
form several different tasks by specifying different input data
to a same executable. Our general philosophy is rather that of
an open distribution, i.e. an integrated suite of codes designed
to be interoperable, much in the spirit of a Linux distribution,
and thus built around a number of core components designed
and maintained by a small group of core developers, plus a
number of auxiliary/complementary codes designed, imple-
mented, and maintained by members of a wider community of
users. The distribution can even be redundant, with different
applications addressing the same problem in different ways;
at the end, the sole requirements to be added to the distribu-
tion are that: i) they are distributed under the same GPL li-
3
cence agreement [20] as other QUANTUM ESPRESSO com-
ponents; ii) they are fully interoperable with the other com-
ponents. Of course, they need to be scienti?cally sound, ver-
i?ed and validated. External contributors are encouraged to
join the QUANTUM ESPRESSO project, if they wish, while
maintaining their own individual distribution and advertise-
ment mode for their software (for instance, by maintaining
individual web sites with their own brand names [26]). To fa-
cilitate this, a web service called qe-forge [27], described
in the next subsection, has been recently put in place.
Interoperability of different components within QUANTUM
ESPRESSO is granted by the use of common formats for the
input, output, and work ?les. In addition, external contributors
are encouraged, but not by any means forced, to use the many
numerical and application libraries on which the core compo-
nents are built. Of course, this general philosophy must be
seen more as an objective to which a very complex software
project tends, rather than a starting point.
One of the main concerns that motivated the birth of the
QUANTUM ESPRESSO project is high performance, both
in serial and in parallel execution. High serial performance
across different architectures is achieved by the systematic
use of standardized mathematical libraries (BLAS, LAPACK
[28], and FFTW [29]) for which highly optimized implemen-
tations exist on many platforms; when proprietary optimiza-
tions of these libraries are not available, the user can compile
the library sources distributed with QUANTUM ESPRESSO.
Optimal performance in parallel execution is achieved through
the design of several parallelization levels, using sophisti-
cated communications algorithms, whose implementation of-
ten does not need to concern the developer, being embedded
in appropriate software layers. As a result the performance
of the key engines (PWscf and CP, see below) may scale on
massively parallel computers up to hundreds or thousands of
processors respectively.
The distribution is organized into a basic set of modules,
libraries, installation utilities, plus a number of directories,
each containing one or more executables, performing speci?c
tasks. The communications between the different executa-
bles take place via data ?les. We think that this kind of ap-
proach lowers the learning barrier for those who wish to con-
tribute to the project. The codes distributed with QUANTUM
ESPRESSO, including many auxiliary codes for the postpro-
cessing of the data generated by the simulations, are easy to
install and to use. The GNU configure and make util-
ities ensure a straightforward installation on many different
machines. Applications are run through text input ?les based
on Fortran namelists, that require the users to specify only an
essential but small subset of the many control variables avail-
able; a specialized graphical user interface (GUI) that is pro-
vided with the distribution facilitates this task for most com-
ponent programs. It is foreseen that in the near future the
design of special APIs (Application Programming Interfaces)
will make it easier to glue different components of the dis-
tribution together and with external applications, as well as
to interface them to other, custom-tailored, GUIs and/or com-
mand interpreters.
The QUANTUM ESPRESSO distribution is written,
mostly, in Fortran-95, with some parts in C or in Fortran-77.
Fortran-95 offers the possibility to introduce advanced pro-
gramming techniques without sacri?cing the performances.
Moreover Fortran is still the language of choice for high-
performance computing and it allows for easy integration of
legacy codes written in this language. A single source tree
is used for all architectures, with C preprocessor options se-
lecting a small subset of architecture-dependent code. Par-
allelization is achieved using the Message-Passing paradigm
and calls to standard MPI (Message Passing Interface) [30]
libraries. Most calls are hidden in a few routines that act as
an intermediate layer, accomplishing e.g. the tasks of sum-
ming a distributed quantity over processors, of collecting dis-
tributed arrays or distributing them across processors, and
to perform parallel three-dimensional Fast Fourier Transform
(FFT). This allows to develop straightforwardly and transpar-
ently new modules and functionalities that preserve the ef?-
cient parallelization backbone of the codes.
A. QE-forge
The ambition of the QUANTUM ESPRESSO project is not
limited to providing highly ef?cient and user-friendly soft-
ware for large-scale electronic-structure calculations and ma-
terials modeling. QUANTUM ESPRESSO aims at promot-
ing active cooperation between a vast and diverse commu-
nity of scientists developing new methods and algorithms in
electronic-structure theory and of end users interested in their
application to the numerical simulation of materials and de-
vices.
As mentioned, the main source of inspiration for the model
we want to promote is the successful cooperative experience
of the GNU/Linux developers’ and users’ community. One of
the main outcomes of this community has been the incorpora-
tion within the GNU/Linux operating system distributions of
third-party software packages, which, while being developed
and maintained by autonomous, and often very small, groups
of users, are put at the disposal of the entire community under
the terms of the GPL. The community, in turn, provides pos-
itive feedback and extensive validation by benchmarking new
developments, reporting bugs, and requesting new features.
These developments have largely bene?ted from the Source-
Forge code repository and software development service [31],
or by other similar services, such as RubyForge, Tigris.org,
BountySource, BerliOS, JavaForge, and GNU Savannah.
Inspired by this model, the QUANTUM ESPRESSO devel-
opers’ and users’ community has set up its own web portal,
named qe-forge [27]. The goal of qe-forge is to com-
plement the traditional web sites of individual scienti?c soft-
ware projects, which are passive instruments of information
retrieval, with a dynamical space for active content creation
and sharing. Its aim is to foster and simplify the coordination
and integration of the programming efforts of heterogeneous
groups and to ease the dissemination of the software tools thus
obtained.
qe-forge provides, through a user-friendly web inter-
face, an integrated development environment, whereby re-
4
searchers can freely upload, manage and maintain their own
software, while retaining full control over it, including the
right of not releasing it. The services so far available in-
clude source-code management software (CVS or SVN repos-
itory), mailing lists, public forums, bug tracking facilities,
download space, and wiki pages for projects’ documentation.
qe-forge is expected to be the main tool by which QUAN-
TUM ESPRESSO end users and external contributors can
maintain QUANTUM ESPRESSO-related projects and make
them available to the community.
III. SHORT DESCRIPTION OF QUANTUM ESPRESSO
QUANTUM ESPRESSO implements a variety of methods
and algorithms aimed at a chemically realistic modeling of
materials from the nanoscale upwards, based on the solution
of density-functional theory (DFT) [2, 3] problem, using a
plane waves (PW) basis set and pseudopotentials (PP) [32] to
represent electron-ion interactions.
The codes are constructed around the use of periodic
boundary conditions, which allows for a straightforward treat-
ment of in?nite crystalline systems, and an ef?cient con-
vergence to the thermodynamic limit for aperiodic but ex-
tended systems, such as liquids or amorphous materials. Fi-
nite systems are also treated using supercells; if required,
open-boundary conditions can be used through the use of the
density-countercharge method [33]. QUANTUM ESPRESSO
can thus be used for any crystal structure or supercell, and
for metals as well as for insulators. The atomic cores can be
described by separable [34] norm-conserving (NC) PPs [35],
ultra-soft (US) PPs [36], or by projector-augmented wave
(PAW) sets [37]. Many different exchange-correlation func-
tionals are available in the framework of the local-density
(LDA) or generalized-gradient approximation (GGA) [38],
plus advanced functionals like Hubbard U corrections and few
meta-GGA [39] and hybrid functionals [40–42]. The latter
is an area of very active development, and more details on
the implementation of hybrid functionals and related Fock-
exchange techniques are given in Appendix F.
The basic computations/simulations that can be performed
include:
• calculation of the Kohn-Sham orbitals and energies [43]
for isolated or extended/periodic systems, and of their
ground-state energies;
• complete structural optimizations of the microscopic
(atomic coordinates) and macroscopic (unit cell) de-
grees of freedom, using Hellmann-Feynman forces [44,
45] and stresses [46];
• ground states for magnetic or spin-polarized system, in-
cluding spin-orbit coupling [47] and non-collinear mag-
netism [48, 49];
• ab-initio molecular dynamics (MD), using either the
Car-Parrinello Lagrangian [50] or the Hellmann-
Feynman forces calculated on the Born-Oppenheimer
(BO) surface [51], in a variety of thermodynamical en-
sembles, including NPT variable-cell [52, 53] MD;
• density-functional perturbation theory (DFPT) [54–56],
to calculate second and third derivatives of the total en-
ergy at any arbitrary wavelength, providing phonon dis-
persions, electron-phonon and phonon-phonon interac-
tions, and static response functions (dielectric tensors,
Born effective charges, IR spectra, Raman tensors);
• location of saddle points and transition states via
transition-path optimization using the nudged elastic
band (NEB) method [57–59];
• ballistic conductance within the Landauer-Büttiker the-
ory using the scattering approach [60];
• generation of maximally localized Wannier functions
[61, 62] and related quantities;
• calculation of nuclear magnetic resonance (NMR) and
electronic paramagnetic resonance (EPR) parameters
[63, 64];
• calculation of K-edge X-ray absorption spectra [65].
Other more advanced or specialized capabilities are described
in the next sections, while ongoing projects (e.g. time-
dependent DFT, many-body perturbation theory) are men-
tioned in the last section. Selected applications were described
in Ref [66]. Several utilities for data postprocessing and inter-
facing to advanced graphic applications are available, allow-
ing e.g. to calculate scanning tunneling microscopy (STM)
images [67], the electron localization function (ELF) [68],
Löwdin charges [69], the density of states (DOS), and pla-
nar [70] or spherical averages of the charge and spin densities
and potentials.
A. Data ?le format
The interoperability of different software components
within a complex project such as QUANTUM ESPRESSO re-
lies on the careful design of ?le formats for data exchange. A
rational and open approach to data ?le formats is also essential
for interfacing applications within QUANTUM ESPRESSO
with third-party applications, and more generally to make the
results of lengthy and expensive computer simulations ac-
cessible to, and reproducible by, the scienti?c community at
large. The need for data ?le formats that make data exchange
easier than it is now is starting to be widely appreciated in the
electronic-structure community. This problem has many as-
pects and likely no simple, "one-size-?ts-all", solution. Data
?les should ideally be
• extensible: one should be able to add some more infor-
mation to a ?le without breaking all codes that read that
?le;
• self-documenting: it should be possible to understand
the contents of a ?le without too much effort;
5
• ef?cient: with data size in the order of GBytes for
large-scale calculations, slow or wasteful I/O should be
avoided.
The current trend in the electronic-structure community seems
to be the adoption of one of the following approaches:
• Structured ?le formats, notably Hierarchical Data For-
mat (HDF) [71] and network Common Data Form
(netCDF) [72], that are widely used since years in other
communities;
• ?le formats based on the Extensible Markup Language
(XML) [73].
It is unlikely that a common, standardized data formats will
ever prevail in our community. We feel that we should fo-
cus, rather than on standardization, on an approach that al-
lows an easy design and usage of simple and reliable con-
verters among different data formats. Prompted by these con-
siderations, QUANTUM ESPRESSO developers have opted
for a simple solution that tries to combine the advantages of
both the above-mentioned approaches. A single ?le contain-
ing all the data of a simulation is replaced by a data direc-
tory, containing several ?les and subdirectories, much in the
same way as it is done in the Mac OS X operating system.
The “head” ?le contains data written with ordinary Fortran
formatted I/O, identi?ed by XML tags. Only data of small
size, such as atomic positions, parameters used in the calcu-
lation, one-electron and total energies, are written in the head
?le. Data of potentially large size, such as PW coef?cients of
Kohn-Sham orbitals, charge density, potentials, are present as
links to separate ?les, written using unformatted Fortran I/O.
Data for each k-point is written to a separate subdirectory.
A lightweight library called iotk, standing for Input/Output
ToolKit [74], is used to read and write the data directory.
Another problem affecting interoperability of PW-PP codes
is the availability of data ?les containing atomic PP’s — one
of the basic ingredients of the calculation. There are many
different types of PP’s, many different codes generating PP’s
(see e.g. Ref [75–77]), each one with its own format. Again,
the choice has fallen on a simple solution that makes it easy
to write converters from and to the format used by QUANTUM
ESPRESSO. Each atomic PP is contained in a formatted ?le
(ef?ciency is not an issue here), described by a XML-like syn-
tax. The resulting format has been dubbed Uni?ed Pseudopo-
tential File (UPF). Several converters from other formats to
the UPF format are available in QUANTUM ESPRESSO.
IV. PARALLELIZATION
Keeping the pace with the evolution of high-end supercom-
puters is one of the guiding lines in the design of QUANTUM
ESPRESSO, with a signi?cant effort being dedicated to port-
ing it to the latest available architectures. This effort is moti-
vated not only by the need to stay at the forefront of architec-
tural innovation for large to very-large scale materials science
simulations, but also by the speed at which hardware features
speci?cally designed for supercomputers ?nd their way into
commodity computers.
The architecture of today’s supercomputers is character-
ized by multiple levels and layers of parallelism: the bot-
tom layer is the one affecting the instruction set of a single
core (simultaneous multithreading, hyperthreading); then one
has parallel processing at processor level (many CPU cores
inside a single processor sharing caches) and at node level
(many processors sharing the same memory inside the node);
at the top level, many nodes are ?nally interconnected with
a high-performance network. The main components of the
QUANTUM ESPRESSO distribution are designed to exploit
this highly structured hardware hierarchy. High performance
on massively parallel architectures is achieved by distributing
both data and computations in a hierarchical way across avail-
able processors, ending up with multiple parallelization levels
[78] that can be tuned to the speci?c application and to the
speci?c architecture. This remarkable characteristic makes it
possible for the main codes of the distribution to run in parallel
on most or all parallel machines with very good performance
in all cases.
More in detail, the various parallelization levels are geared
into a hierarchy of processor groups, identi?ed by different
MPI communicators. In this hierarchy, groups implement-
ing coarser-grained parallel tasks are split into groups im-
plementing ?ner-grained parallel tasks. The ?rst level is im-
age parallelization, implemented by dividing processors into
n
image
groups, each taking care of one or more images (i.e.
a point in the con?guration space, used by the NEB method).
The second level is pool parallelization, implemented by fur-
ther dividing each group of processors into n
pool
pools of
processors, each taking care of one or more k-points (see
Appendix D). The third level is plane-wave parallelization,
implemented by distributing real- and reciprocal-space grids
across the n
PW
processors of each pool. The ?nal level is task
group parallelization [79], in which processors are divided
into n
task
task groups of n
FFT
= n
PW
/n
task
processors,
each one taking care of different groups of electron states to
be Fourier-transformed, while each FFT is parallelized inside
a task group. A further paralellization level, linear-algebra,
coexists side-to-side with plane-wave parallelization, i.e. they
take care of different sets of operations, with different data dis-
tribution. Linear-algebra parallelization is implemented both
with custom algorithms and using ScaLAPACK [80], which
on massively parallel machines yield much superior perfor-
mances. Table I contains a summary of the ?ve levels cur-
rently implemented. With the recent addition of the two last
levels, most parallelization bottlenecks have been removed,
both computations and data structures are fully distributed,
scalability on parallel machines is only limited by the phys-
ical sizes of the system being simulated. Scalability does not
yet extend to tens of thousands of processors as in specially-
crafted code QBox[81], but excellent scalability on up to 4800
processors has been demonstrated (see Fig. 1) even for cases
where coarse-grained parallelization doesn’t help.
Scalability is thus guaranteed for large-scale simulations.
This being said, it is obvious that the size and speci?c nature
of the speci?c application sets quite natural limits to the max-
6
TABLE I: Summary of parallelization levels in QUANTUM ESPRESSO.
group distributed quantities communications performance
image NEB images very low linear CPU scaling,
fair to good load balancing;
does not distribute RAM
pool k-points low almost linear CPU scaling,
fair to good load balancing;
does not distribute RAM
plane-wave plane waves, G-vector coef?cients, high good CPU scaling,
R-space FFT arrays good load balancing,
distributes most RAM
task FFT on electron states high improves load balancing
linear algebra subspace Hamiltonians very high improves scaling,
and constraints matrices distributes more RAM
1000
1500
2000
2500
3000
3500
4000
4500
5000
0 1000 2000 3000 4000 5000
C
P
U
t
i
m
e
(
s
)
Number of CPUs
Scalability for > 1000 processors
PSIWAT
CNT (1)
CNT (2)
FIG. 1: Scalability as a function of the number of processors for
large-scale calculations. PSIWAT: gold surface covered by thiols in
interaction with water, 10.59 × 20.53 × 32.66 Å
3
cell, 587 atoms,
2552 electrons, PWscf code with n
pool
= 4 and n
tg
= 4 on A Cray
XT 4. CNT: porphyrin-functionalized nanotube, 1532 atoms, 5232
electrons; case (1), PWscf on a Cray XT 4; case (2), CP on a Cray
XT3(2) codes. CPU times for PSIWAT and CNT (2) were rescaled
so that they to fall in the same range as for CNT (1).
imum number of processors up to which the performances of
the various codes are expected to scale. For instance, the num-
ber of images in a NEB calculation sets a natural limit to the
level of image groups, or the number of electronic bands sets
a limit for the parallelization of the linear algebra operations.
Moreover some numerical algorithms scale better than others.
For example, the use of norm-conserving pseudopotentials al-
lows for a better scaling than ultrasoft pseudopotentials for
a same system, because a larger plane wave basis set and a
larger real- and reciprocal-space grids are required in the for-
mer case. On the other hand, using ultrasoft pseudopotentials
is generally faster because the use of a smaller basis set is
obviously more ef?cient, even though the overall parallel per-
formance may be not as good.
The efforts of the QUANTUM ESPRESSO developers’
50
100
150
200
250
300
350
0 100 200 300 400 500
C
P
U
t
i
m
e
(
s
)
Number of CPUs
CP scalability
BG/P 1 task
BG/P 4 tasks
BG/P 8 tasks
Altix 1 task
Altix 4 tasks
Altix 8 tasks
FIG. 2: CPU time t
cpu
(s) per electronic time step (CP code) on
an IBM Blue Gene/P (BG/P) and a SGI Altix for a fragment of an
A??peptide in water containing 838 atoms and 2311 electrons in a
22.1×22.9×19.9Å
3
cell, as a function of the number of processors,
for different numbers n
task
of task groups.
team are not limited to the performance on massively paral-
lel architectures. Simulations on systems containing several
hundreds of atoms are by now quite standard (see Fig. 2
for an example). Special attention is also paid to optimize
the performances for simulations of intermediate size (on sys-
tems comprising from several tens to a few hundreds inequiv-
alent atoms), to be performed on medium-size clusters, read-
ily available to many groups [82]. In particular, the QUAN-
TUM ESPRESSO developers’ team is now working to better
exploit new hardware trends, particularly in the ?eld of mul-
ticore architectures. By the time this paper will be printed,
mixed OpenMP-MPI parallelization will be available for en-
hanced performance on modern multicore CPU’s. Looking
ahead, future developments will likely focus on hybrid system
with hardware accelerators (GPUs and cell co-processors).
7
V. QUANTUM ESPRESSO PACKAGES
The complete QUANTUM ESPRESSO distribution is rel-
ative large: about 240,000 lines of code, excluding copies of
the external libraries. With such a sizable code basis, mod-
ularization becomes necessary. QUANTUM ESPRESSO is
presently divided into several executables, performing differ-
ent types of calculations, although some of them have over-
lapping functionalities. Typically there is a single set of func-
tions/subroutines or a single Fortran 90 module that performs
each speci?c task (e.g. matrix diagonalizations, or potential
updates), but there are still important exceptions to this rule,
re?ecting the different origin and different styles of the orig-
inal components. QUANTUM ESPRESSO has in fact been
built out of the merge and re-engineering of different pack-
ages, that were independently developed for several years. In
the following, the main components are brie?y described.
A. PWscf
PWscf implements an iterative approach to reach self-
consistency, using at each step iterative diagonalization tech-
niques, in the framework of the plane-wave pseudopotential
method. An early version of PWscf is described in Ref [83].
Both separable NC-PPs and US-PPs are implemented; re-
cently, also the projector augmented-wave method [37] has
been added, largely following the lines of Ref [84] for its
implementation. In the case of US-PPs, the electronic wave
functions can be made smoother at the price of having to aug-
ment their square modulus with additional contributions to re-
cover the actual physical charge densities. For this reason,
the charge density is more structured than the square of the
wavefunctions, and requires a larger energy cutoff for its plane
wave expansion (typically, 6 to 12 times larger; for a NC-PP,
a factor of 4 would be mathematically suf?cient). Hence, dif-
ferent real-space Fourier grids are introduced - a "soft" one
that represents the square of electronic wave functions, and a
"hard" one that represents the charge density [82, 85]. The
augmentation terms can be added either in reciprocal space
(using an exact but expensive algorithm) or directly in real
space (using an approximate but faster algorithm that exploits
the local character of the augmentation charges).
PWscf can use the established LDA and GGA exchange-
correlation functionals, including spin-polarization within the
scheme proposed in Ref [86] and can treat non-collinear
magnetism[48, 49] as e.g. induced by relativistic effects (spin-
orbit interactions) [87, 88] or by complex magnetic interac-
tions (e.g. in the presence of frustration). DFT + Hubbard U
calculations [89] are implemented for a simpli?ed (“no-J”)
rotationally invariant form [90] of the Hubbard term. Other
advanced functionals include TPSS meta-GGA[39], function-
als with ?nite-size corrections [91], and the PBE0 [40] and
B3LYP [41, 42] hybrids.
Self-consistency is achieved via the modi?ed Broyden
method of Ref [92], with some further re?nements that are
detailed in the Appendix. The sampling of the Brillouin Zone
(BZ) can be performed using either special [93, 94] k-points
provided in input or those automatically calculated starting
froma uniform grid. Crystal symmetries are automatically de-
tected and exploited to reduce computational costs, by restrict-
ing the sampling of the BZ to the irreducible wedge alone.
When only the ? point (k = 0) is used, advantage is taken
of the real character of the orbitals, allowing to store just half
of the Fourier components. BZ integrations in metallic sys-
tems are dealt with smearing/broadening techniques, such as
Fermi-Dirac, Gaussian, Methfessel-Paxton [95], and Marzari-
Vanderbilt cold smearing [96]; the tetrahedron method [97] is
also implemented. Note that ?nite-temperature effects on the
electronic properties can be easily accounted for by using the
Fermi-Dirac smearing not as a mathematical device, but as a
practical way of implementing the Mermin density-functional
approach [98].
Structural optimizations are performed using the Broyden-
Fletcher-Goldfarb-Shanno (BFGS) algorithm [99–101] or
damped dynamics; these can involve both the internal, micro-
scopic degrees of freedom (i.e. the atomic coordinates) and/or
the macroscopic ones (shape and size of the unit cell). The
calculation of minimum-energy paths, activation energies, and
transition states uses the Nudged Elastic Band (NEB) method
[57]. Potential energy surfaces as a function of suitably cho-
sen collective variables can be studied using Laio-Parrinello
metadynamics [102].
Microcanonical (NVE) MD is performed on the BO sur-
face, i.e. achieving electron self-consistency at each time step,
using the Verlet algorithm[103]. Canonical (NVT) dynamics
can be performed using velocity rescaling, or Anderson’s or
Berendsen’s thermostats [104]. Constant-pressure (NPT) MD
is performed by adding additional degrees of freedom for the
cell size and volume, using either the Parrinello-Rahman La-
grangian [105] or the invariant Lagrangian of Wentzcovitch
[53].
The effects of ?nite macroscopic electric ?elds on the elec-
tronic structure of the ground state can be accounted for ei-
ther through the method of Ref [106, 107] based on the Berry
phase, or (for slab geometries only) through a sawtooth ex-
ternal potential [108, 109]. A quantum fragment can be em-
bedded in a complex electrostatic environment that includes a
model solvent [110] and a counterion distribution [111], as is
typical of electrochemical systems.
B. CP
The CP code is the massively-parallel module for Car-
Parrinello ab-initio MD. CP can use both NC PPs [112] and
US PPs [85, 113]. In the latter case, the electron density
is augmented through a Fourier interpolation scheme in real
space (“box grid”) [82, 85] that is particular ef?cient for large
scale calculations. CP implements the same functionals as in
PWscf, with the exception of hybrid functionals; a simpli-
?ed one-electron self-interaction correction (SIC)[114] is also
available. The Car-Parrinello Lagrangian can be augmented
with Hubbard U corrections [115], or Hubbard-based penalty
functionals to impose arbitrary oxidation states [116].
Since the main applications of CP are for large systems
8
without translational symmetry (e.g. liquids, amorphous ma-
terials), Brillouin zone sampling is restricted to the ? point of
the supercell, allowing for real instead of complex wavefunc-
tions. Metallic systems can be treated in the framework of
“ensemble DFT” [117].
In the Car-Parrinello algorithm, microcanonical (NVE) MD
is performed on both electronic and nuclear degrees of free-
dom, treated on the same footing, using the Verlet algorithm.
The electronic equations of motion are accelerated through a
preconditioning scheme [118]. Constant-pressure (NPT) MD
is performed using the Parrinello-Rahman Lagrangian [105]
and additional degrees of freedom for the cell. Nosé-Hoover
thermostats [119] and Nosé-Hoover chains [120] allow to per-
form simulations in the different canonical ensembles.
CP can also be used to directly minimize the energy func-
tional to self-consistency while keeping the nuclei ?xed, or to
perform structural minimizations of nuclear positions, using
the “global minimization” approaches of in Refs. [121, 122],
and damped dynamics or conjugate-gradients on the elec-
tronic or ionic degrees of freedom. It can also perform NEB
and metadynamics calculations.
Finite homogeneous electric ?elds can be accounted for us-
ing the Berry’s phase method, adapted to systems with the
? point only [106]. This advanced feature can be used in
combination with MD to obtain the infrared spectra of liq-
uids [106, 123], the low- and high-frequency dielectric con-
stants [106, 124] and the coupling factors required for the cal-
culation of vibrational properties, including infrared, Raman
[125–127], and hyper-Raman [128] spectra.
C. PHonon
The PHonon code implements density-functional pertur-
bation theory (DFPT) [54–56] for the calculation of second-
and third-order derivatives of the energy with respect to
atomic displacements and to electric ?elds. The global min-
imization approach [129, 130] is used for the special case
of normal modes in ?nite (molecular) systems, where no BZ
sampling is required, while a self-consistent procedure [55] is
used in the general case, with the distinctive advantage that
the response to a perturbation of any arbitrary wavelength
can be calculated with a computational cost that is propor-
tional to that of the unperturbed system. Thus, response at
any wavevector, including very long-wavelength ones, can be
inexpensively calculated. This latter approach, and the tech-
nicalities involved in the calculation of effective charges and
interatomic force constants, are described in detail in Refs.
[55, 131].
Symmetry is fully exploited in order to reduce the amount
of computation. Lattice distortions transforming according
to irreducible representations of small dimensions are gener-
ated ?rst. The charge-density response to these lattice distor-
tions is then sampled at a number of discrete k-points in the
BZ, which is reduced according to the symmetry of the small
group of the phonon wavevector q. The grid of q needed for
the calculation of interatomic force constants reduces to one
wavevector per star: the dynamical matrices at the other q
vectors in the star are generated using the symmetry opera-
tions of the crystal. This approach allows us to speed up the
calculation without the need to store too much data for sym-
metrization.
The calculation of second-order derivatives of the energy
works also for US PP [132, 133] and for all GGA ?avors
[134, 135] used in PWscf and in CP. The extension of
PHonon to PAW [136], to noncollinear magnetism and to
fully relativistic US PPs which include spin-orbit coupling
[137] will be available by the time this paper will be printed.
Advanced features of the PHonon code include the calcu-
lation of third order derivatives of the energy, such as electron-
phonon or phonon-phonon interaction coef?cients. Electron-
phonon interactions are straightforwardly calculated from the
response of the self-consistent potential to a lattice distortion.
This involves a numerically-sensitive “double-delta” integra-
tion at the Fermi energy, that is performed using interpolations
on a dense k-point grid. Interpolation techniques based on
Wannier functions [138] will speed up considerably these cal-
culations. The calculation of the anharmonic force constants
from third-order derivatives of the electronic ground-state en-
ergy is described in Ref. [139] and is performed by a separate
code called d3. Static Raman coef?cients are calculated us-
ing the second-order response approach of Refs. [140, 141].
Both third-order derivatives and Raman-coef?cients calcula-
tions are currently implemented only for NC PPs.
D. atomic
The atomic code performs three different tasks: i) solu-
tion of the self-consistent all-electron radial Kohn-Shamequa-
tions (with a Coulomb nuclear potential and spherically sym-
metric charge density); ii) generation of NC PPs, of US PPs,
or of PAW data-sets; iii) test of the above PPs and data-sets.
The three tasks can be either separately executed or performed
in a single run. Three different all-electron equations are avail-
able: i) the non relativistic radial Kohn and Sham equations,
ii) the scalar relativistic approximation to the radial Dirac
equations [142], iii) the radial Dirac-like equations derived
within relativistic density functional theory [143, 144]. For i)
and ii) atomic magnetism is dealt within the local spin density
approximation, i.e. assuming an axis of magnetization. The
atomic code uses the same exchange and correlation energy
routines of PWscf and can deal with the same functionals.
The code is able to generate NC PPs directly in separable
form (also with multiple projectors per angular momentum
channel) via the Troullier-Martins [145] or the Rappe-Rabe-
Kaxiras-Joannopoulos [146] pseudization. US PPs can be
generated by a two-step pseudization process, starting from
a NC PPs, as described in Ref. [147], or using the solutions
of the all-electron equation and pseudizing the augmentation
functions [85]. The latter method is used also for the PAW
data-set generation. The generation of fully relativistic NC
and US PPs including spin-orbit coupling effects is also avail-
able. Converters are available to translate pseudopotentials
encoded in different formats (e.g. according to the Fritz-
Haber [75] or Vanderbilt [76] conventions) into the UPF for-
9
mat adopted by QUANTUM ESPRESSO.
Transferability tests can be made simultaneously for sev-
eral atomic con?gurations with or without spin-polarization,
by solving the non relativistic radial Kohn-Sham equations
generalized for separable nonlocal PPs and for the presence
of an overlap matrix.
E. PWcond
The PWcond code implements the scattering approach pro-
posed by Choi and Ihm [60] for the study of coherent electron
transport in atomic-sized nanocontacts within the Landauer-
Büttiker theory. Within this scheme the linear response bal-
listic conductance is proportional to the quantum-mechanical
electron transmission at the Fermi energy for an open quan-
tum system consisting of a scattering region (e.g., an atomic
chain or a molecule with some portions of left and right leads)
connected ideally from both sides to semi-in?nite metallic
leads. The transmission is evaluated by solving the Kohn-
Sham equations with the boundary conditions that an electron
coming fromthe left lead and propagating rightwards gets par-
tially re?ected and partially transmitted by the scattering re-
gion. The total transmission is obtained by summing all trans-
mission probabilities for all the propagating channels in the
left lead. As a byproduct of the method, the PWcond code
provides the complex band structures of the leads, that is the
Bloch states with complex k
z
in the direction of transport, de-
scribing wave functions exponentially growing or decaying in
the z direction. The original method formulated with NC PPs
has been generalized to US PPs both in the scalar relativistic
[148] and in the fully relativistic forms [149].
F. GIPAW
The GIPAW code allows for the calculation of physical pa-
rameters measured by i) nuclear magnetic resonance (NMR)
in insulators (the electric ?eld gradient (EFG) tensors and the
chemical shift tensors), and by ii) electronic paramagnetic res-
onance (EPR) spectroscopy for paramagnetic defects in solids
or in radicals (the hyper?ne tensors and the g-tensor). The
code also computes the magnetic susceptibility of nonmag-
netic insulators. The code is based on the PW-PP method,
and uses many subroutines of PWscf and of PHonon. The
code is currently restricted to NC PPs. All the NMR and
EPR parameters depend on the detailed shape of the elec-
tronic wave-functions near the nuclei and thus require the re-
construction of the all-electron wave-functions from the PP
wave-functions. For the properties de?ned at zero external
magnetic ?eld, namely the EFG and the hyper?ne tensors,
such reconstruction is performed as a post-processing step of
a self-consistent calculation using the PAW reconstruction, as
described for the EFG in Ref.[150] and for the hyper?ne ten-
sor in Ref.[151]. The g-tensor, the NMR chemical shifts and
the magnetic susceptibility are obtained from the orbital lin-
ear response to an external uniform magnetic ?eld. In pres-
ence of a magnetic ?eld the PAW method is no more gauge-
and translationally invariant. Gauge and translational invari-
ances are restored by using the gauge including projector aug-
mented wave (GIPAW) method [63, 64] both i) to describe
in the PP Hamiltonian the coupling of orbital degrees of free-
dom with the external magnetic ?eld, and ii) to reconstruct the
all-electron wave-functions, in presence of the external mag-
netic ?eld. In addition, the description of a uniform magnetic
?eld within periodic boundary condition is achieved by con-
sidering the long wave-length limit of a sinusoidally modu-
lated ?eld in real space [152, 153]. The NMR chemical shifts
are computed following the method described in Ref.[63], the
g-tensor following Ref.[154] and the magnetic susceptibility
following Refs.[63, 152]. Recently, a “converse” approach
to calculate chemical shifts has also been introduced [155],
based on recent developments on the Berry-phase theory of or-
bital magnetization; since it does not require a linear-response
calculation, it can be straightforwardly applied to arbitrarily
complex exchange-correlation functionals, and to very large
systems, albeit at computational costs that are proportional to
the number of chemical shifts that need to be calculated.
G. XSPECTRA
The XSPECTRA code allows for the calculation of K-edge
X-ray absorption spectra. The code calculates both the dipo-
lar and the quadrupolar matrix elements. The code uses the
self-consistent charge density produced by PWscf and acts
as a post-processing tool. The all-electron wavefunction is
reconstructed using the PAW method and its implementation
in the GIPAW code. The presence of a core-hole in the ?nal
state of the X-ray absorption process is simulated by using a
pseudopotential for the absorbing atom with a hole in the 1s
state. The calculation of the charge density is performed on a
supercell with one absorbing atom. From the self-consistent
charge density, the X-ray absorption spectra are obtained us-
ing the Lanczos method and a continuous fraction expansion
[65]. The advantage of this approach is that it is not necessary
to calculate many empty states even if features of the spec-
trum at very high energy are required. Correlation effects can
be simulated in a mean-?eld way using the Hubbard U cor-
rection [90] that has been included in the XSPECTRA code in
Ref. [156]. Currently the code works only with NC PPs and is
limited to collinear magnetism. Its extension to noncollinear
magnetismand to fully relativistic US PPs which include spin-
orbit coupling is under development.
H. Wannier90
Wannier90 [26, 157] is a code that calculates maximally-
localized Wannier functions in insulators or metals—
according to the algorithms described in Refs. [61, 62]—and
a number of properties that can be conveniently expressed in
a Wannier basis. The code is developed and maintained inde-
pendently by a Wannier development group [157] and can be
taken as a representative example of the philosophy described
earlier, where a project maintains its own individual distribu-
10
tion but provides full interoperability with the core compo-
nents of QUANTUM ESPRESSO in this case PWscf or CP.
These codes are in fact used as “quantum engines” to produce
the data onto which Wannier90 operates. The need to pro-
vide transparent protocols for interoperability has in turn fa-
cilitated the interfacing of wannier90 with other quantum
engines [14, 21], fostering a collaborative engagement with
the broader electronic-structure community that is also in the
spirit of QUANTUM ESPRESSO .
Wannier90 requires as input the scalar products between
wavefunctions at neighboring k-points, where these latter
form uniform meshes in the Brillouin zone. Often, it is
also convenient to provide scalar products between wavefunc-
tions and trial, localized real-space orbitals—these are used to
guide the localization procedure towards a desired, physical
minimum. As such, the code is not tied to a representation
of the wavefunctions in any particular basis—for PWscf and
CP a postprocessing utility is in charge of calculating these
scalar products using the plane-wave basis set of QUANTUM
ESPRESSO and either NC-PPs or US-PPs. Whenever ?
sampling is used, the simpli?ed algorithm of Ref. [158] is
adopted.
Besides calculating maximally localized Wannier func-
tions, the code is able to construct the Hamiltonian matrix
in this localized basis, providing a chemically accurate, and
transferable, tight-binding representation of the electronic-
structure of the system. This, in turn, can be used to construct
Green’s functions and self-energies for ballistic transport cal-
culations [159, 160], to determine the electronic-structure and
density-of-states of very large scale structures [160], to inter-
polate accurately the electronic band structure (i.e. the Hamil-
tonian) across the Brillouin zone [160, 161], or to interpolate
any other operator [161]. These latter capabilities are espe-
cially useful for the calculation of integrals that depend sen-
sitively on a submanifold of states; common examples come
from properties that depend sensitively on the Fermi surface,
such as electronic conductivity, electron-phonon couplings
Knight shifts, or the anomalous Hall effect. A related by-
product of Wannier90 is the capability of downfolding a
selected, physically signi?cant manifold of bands into a mini-
mal but accurate basis, to be used for model Hamiltonians that
can be treated with complex many-body approaches.
I. PostProc
The PostProc module contains a number of codes for
postprocessing and analysis of data ?les produced by PWscf
and CP. The following operations can be performed:
• Interfacing to graphical and molecular graphics appli-
cations. Charge and spin density, potentials, ELF [68]
and STM images [67] are extracted or calculated and
written to ?les that can be directly read by most com-
mon plotting programs, like xcrysden [162] and VMD
[163].
• Interfaces to other codes that use DFT results from
QUANTUM ESPRESSO for further calculations, such
as e.g.: pw2wannier90, an interface to the
wannier90 library and code [26, 157] (also in-
cluded in the QUANTUM ESPRESSO distribution);
pw2casino.f90, an interface to the casino quan-
tum Monte Carlo code [164]; wannier_ham.f90, a
tool to build a tight-binding representation of the KS
Hamiltonian to be used by the dmft code [165] (avail-
able at the qe-forge site); pw_export.f90, an in-
terface to the GW code SaX [166]; pw2gw.f90, an
interface to code DP [167] for dielectric property calcu-
lations, and to code EXC [168] for excited-state proper-
ties.
• Calculation of various quantities that are useful for the
analysis of the results. In addition to the already men-
tioned ELF and STM, one can calculate projections
over atomic states (e.g. Löwdin charges [69]), DOS and
Projected DOS (PDOS), planar and spherical averages,
and the complex macroscopic dielectric function in the
random-phase approximation (RPA).
J. PWgui
PWgui is the graphical user interface (GUI) for the
PWscf, PHonon, d3 and atomic codes as well as for
some of the main post-processing routines (e.g. pp.x and
projwfc.x). PWgui is an input ?le builder whose main
goal is to lower the learning barrier for the newcomer, who
has to struggle with the input syntax. Its event-driven mecha-
nism automatically adjusts the display of required input ?elds
(i.e. enables certain sets of widgets and disables others) to the
speci?c cases selected. It enables a preview of the format of
the (required) input ?le records for a given type of calculation.
The input ?les created by PWgui are syntactically correct (al-
though they can still be physically meaningless). It is possible
to upload previously generated input ?les for syntax checking
and/or to modify them. It is also possible to run calculations
from within the PWgui. In addition, PWgui can also use
the external xcrysden program [162] for the visualization
of molecular and/or crystal structures from the speci?ed in-
put data and for the visualization of properties (e.g. charge
densities or STM images).
As the QUANTUM ESPRESSO codes evolve, the input
?le syntax expands as well. This implies that PWgui has
to be continuously adapted. To effectively deal with such is-
sue, PWgui uses the GUIB concept [169]. GUIB builds on
the consideration that the input ?les for numerical simulation
codes have a rather simple structure and it exploits this sim-
plicity by de?ning a special meta-language with two purposes:
The ?rst is to de?ne the input-?le syntax, and the second is to
simultaneously automate the construction of the GUI on the
basis of such a de?nition.
A similar strategy has been recently adopted for the de-
scription of the QUANTUM ESPRESSO input ?le formats.
A single de?nition/description of a given input ?le serves as i)
a documentation per-se, ii) as a PWgui help documentation,
11
FIG. 3: Snapshot of the PWgui application. Left: PWgui’s main window; right: preview of speci?ed input data in text mode.
and iii) as a utility to synchronize the PWgui with up-to-date
input ?le formats.
VI. PERSPECTIVES AND OUTLOOK
Further developments and extensions of QUANTUM
ESPRESSO will be driven by the needs of the community
using it and working on it. Many of the soon-to-come ad-
ditions will deal with excited-state calculations within time-
dependent DFT (TDDFT [170, 171]) and/or many-body per-
turbation theory [172]. A new approach to the calculation of
optical spectra within TDDFT has been recently developed
[173], based on a ?nite-frequency generalization of density-
functional perturbation theory [54, 55], and implemented
in QUANTUM ESPRESSO. Another important development
presently under way is an ef?cient implementation of GW cal-
culations for large systems (whose size is of the order of a few
hundreds inequivalent atoms) [174]. The implementation of
ef?cient algorithms for calculating correlation energies at the
RPA level is also presently under way [175–177]. It is fore-
seen that by the time this paper will appear, many of these
developments will be publicly released.
It is hoped that many newfunctionalities will be made avail-
able to QUANTUM ESPRESSO users by external groups who
will make their own software compatible/interfaceable with
QUANTUM ESPRESSO. At the time of the writing of the
present paper, third-party scienti?c software available to the
QUANTUM ESPRESSO users’ community include: yambo,
a general-purpose code for excited-state calculations within
many-body perturbation theory [178]; casino, a code for
electronic-structure quantum Monte Carlo simulations [164];
wannier90, a code and a library for generating maximally
localized Wannier functions [26, 61]; want, a code for the
simulation of ballistic transport in nanostructures, based on
Wannier functions [179]; xcrysden, a molecular graphics
application, specially suited for periodic structures [162]. The
qe-forge portal is expected to burst the production and
availability of third-party software compatible with QUAN-
TUM ESPRESSO. Among the projects already available, or
soon-to-be available, on qe-forge, we mention: SaX [166],
an open-source project implementing state-of-the-art many-
body perturbation theory methods for excited states; dmft
[165], a code to perform Dynamical Mean-Field Theory cal-
culations on top of a tight-binding representation of the DFT
band structure; qha, a set of codes for calculating thermal
properties of materials within the quasi-harmonic approxima-
tion [180]; pwtk, a fully functional Tcl scripting interface to
PWscf [181].
Efforts towards better interoperability with third-party soft-
ware will be geared towards releasing accurate speci?cations
for data structures and data ?le formats and providing inter-
faces to and from other codes and packages used by the sci-
enti?c community. Further work will be also devoted to the
extension to the US-PPs and PAW schemes of the parts of
QUANTUM ESPRESSO that are now limited to NC-PPs.
The increasing availability of massively parallel machines
will likely lead to an increased interest towards large-scale cal-
culations. The ongoing effort in this ?eld will continue. A
special attention will be paid to the requirements imposed by
the architecture of the new machines, in particular multicore
CPUs, for which a mixed OpenMP-MPI approach seems to
be the only viable solution yielding maximum performances.
Grid computing and the commoditization of computer cluster
will also lead to great improvements in high-throughput cal-
culations for materials design and discovery.
12
VII. ACKNOWLEDGMENTS
The QUANTUM ESPRESSO project is an initiative of the
CNR-INFM DEMOCRITOS National Simulation Center in
Trieste (Italy) and its partners, in collaboration with MIT,
Princeton University, the University of Minnesota, the Ecole
Polytechnique Fédérale de Lausanne, the Université Pierre et
Marie Curie in Paris, the Jožef Stefan Institute in Ljubljana,
and the S3 research center in Modena. Many of the ideas
embodied in the QUANTUM ESPRESSO codes have ?our-
ished in the very stimulating environment of the International
School of Advanced Studies, SISSA, where Democritos is
hosted, and have bene?tted from the ingenuity of generations
of students and young postdocs. SISSA and the National Su-
percomputing Center CINECA are currently providing valu-
able support to the QUANTUM ESPRESSO project.
VIII. APPENDIX
This appendix contains the description of some algorithms
used in QUANTUM ESPRESSO that have not been docu-
mented elsewhere.
A. Self-consistency
The problem of ?nding a self-consistent solution to the
Kohn-Sham equations can be recast into the solution of a non-
linear problem
x = F[x], x = (x
1
, x
2
, . . . , x
N
), (1)
where vector x contains the N Fourier components or real-
space values of the charge density ? or the Kohn-Sham po-
tential V (the sum of Hartree and exchange-correlation poten-
tials); F[x
(in)
] is a functional of the input charge density or
potential x
(in)
, yielding the output vector x
(out)
via the so-
lution of Kohn-Sham equations. A solution can be found via
an iterative procedure. PWscf uses an algorithm based on
the modi?ed Broyden method [92] in which x contains the
components of the charge density in reciprocal space. Mixing
algorithms typically ?nd the optimal linear combination of a
few x
(in)
from previous iterations, that minimizes some suit-
ably de?ned norm||x
(out)
?x
(in)
||, vanishing at convergence,
that we will call in the following “scf norm”.
Ideally, the scf norm is a measure of the self-consistency
error on the total energy. Let us write an estimate of the latter
for the simplest case: an insulator with NC PPs and simple
LDA or GGA. At a given iteration we have
_
?
¯h
2
2m
?
2
+ V
ext
(r) + V
(in)
(r)
_
?
i
(r) =
i
?
i
(r), (2)
where
i
and ?
i
are Kohn-Sham energies and orbitals re-
spectively, i labels the occupied states, V
ext
is the sum of
the PPs of atomic cores (written for simplicity as a local po-
tential), the input Hartree and exchange-correlation potential
V
(in)
(r) = V
Hxc
[?
(in)
(r)] is a functional of the input charge
density ?
(in)
. The output charge density is given by
?
(out)
(r) =
i
|?
i
(r)|
2
. (3)
Let us compare the DFT energy calculated in the standard
way:
E =
i
_
?
?
i
(r)
_
?
¯h
2
2m
?
2
+ V
ext
(r)
_
?
i
(r)dr
+E
Hxc
[?
(out)
], (4)
where E
Hxc
is the Hartree and exchange-correlation en-
ergy, with the Harris-Weinert-Foulkes functional form, which
doesn’t use ?
(out)
:
E
=
i
_
?
?
i
(r)
_
?
¯h
2
2m
?
2
+ V
ext
(r) + V
(in)
(r)
_
?
i
(r)dr
?
_
?
(in)
V
(in)
(r) + E
Hxc
[?
(in)
] (5)
Both forms are variational, i.e. the ?rst-order variation of the
energy with respect to the charge density vanish, and both
converge to the same result when self-consistency is achieved.
Their difference can be approximated by the following expres-
sion, in which only the dominant Hartree term is considered:
E ?E
1
2
_
??(r)??(r
)
|r ?r
|
drdr
=
1
2
_
??(r)?V
H
(r
)dr (6)
where ?? = ?
(out)
??
(in)
and ?V
H
is the Hartree potential
generated by ??. Moreover it can be shown that, when ex-
change and correlation contributions to the electronic screen-
ing do not dominate over the electrostatic ones, this quantity
is an upper bound to the self-consistent error incurred when
using the standard form for the DFT energy. We therefore
take this term, which can be trivially calculated in reciprocal
space, as our squared scf norm:
||?
(out)
??
(in)
||
2
=
4?e
2
?
G
|??(G)|
2
G
2
, (7)
where G are the vectors in reciprocal space and ? is the vol-
ume of the unit cell.
Once the optimal linear combination of ?
(in)
from previous
iterations (typically 4 to 8) is determined, one adds a step in
the new search direction that is, in the simplest case, a fraction
of the optimal ?? or, taking advantage of some approximate
electronic screening[182], a preconditioned ??. In particular,
the simple, Thomas-Fermi, and local Thomas-Fermi mixing
described in Ref.[182] are implemented and used.
The above algorithm has been extended to more sophisti-
cated calculations, in which the x vector introduced above
may contain additional quantities: for DFT+U, occupancies of
13
atomic correlated states; for meta-GGA, kinetic energy den-
sity; for PAW, the quantities
i
?
i
|?
n
?
m
|?
i
, where the
? functions are the atomic-based projectors appearing in the
PAWformalism. The scf normis modi?ed accordingly in such
a way to include the additional variables in the estimated self-
consistency error.
B. Iterative diagonalization
During self-consistency one has to solve the generalized
eigenvalue problem for all N occupied states
H?
i
=
i
S?
i
, i = 1, . . . , N (8)
in which both H (the Hamiltonian) and S (the overlap ma-
trix) are available as operators (i.e. H? and S? products
can be calculated for a generic state ? ). Eigenvectors are
normalized according to the generalized orthonormality con-
straints ?
i
|S|?
j
= ?
ij
. This problem is solved using itera-
tive methods. Currently PWscf implements a block Davidson
algorithm and an alternative algorithm based on band-by-band
minimization using conjugate gradient.
1. Davidson
One starts from an initial set of orthonormalized trial or-
bitals ?
(0)
i
and of trial eigenvalues
(0)
i
= ?
(0)
i
|H|?
(0)
i
. The
starting set is typically obtained from the previous scf itera-
tion, if available, and if not, from the previous time step, or
optimization step, or from a superposition of atomic orbitals.
We introduce the residual vectors
g
(0)
i
= (H ?
(0)
i
S)?
(0)
i
, (9)
a measure of the error on the trial solution, and the correction
vectors ??
(0)
i
= Dg
(0)
i
, where D is a suitable approximation
to (H ?
(0)
i
S)
?1
. The eigenvalue problem is then solved in
the 2N-dimensional subspace spanned by the reduced basis
set ?
(0)
, formed by ?
(0)
i
= ?
(0)
i
and ?
(0)
i+N
= ??
(0)
i
:
2N
k=1
(H
jk
?
i
S
jk
)c
(i)
k
= 0, (10)
where
H
jk
= ?
(0)
j
|H|?
(0)
k
, S
jk
= ?
(0)
j
|S|?
(0)
k
. (11)
Conventional algorithms for matrix diagonalization are used
in this step. A new set of trial eigenvectors and eigenvalues is
obtained:
?
(1)
i
=
2N
j=1
c
(i)
j
?
(0)
j
,
(1)
i
= ?
(1)
i
|H|?
(1)
i
(12)
and the procedure is iterated until a satisfactory convergence is
achieved. Alternatively, one may enlarge the reduced basis set
with the new correction vectors ??
(1)
i
= Dg
(1)
i
, solve a 3N-
dimensional problem, and so on, until a pre?xed size of the
reduced basis set is reached. The latter approach is typically
slightly faster at the expenses of a larger memory usage.
The operator D must be easy to estimate. A natural choice
in the PWbasis set is a diagonal matrix, obtained keeping only
the diagonal term of the Hamiltonian:
k +G|D|k +G
=
?
GG
k +G|H ?S|k +G
(13)
where k is the Bloch vector of the electronic states under con-
sideration, |k+G
denotes PWs, an estimate of the highest
occupied eigenvalue. Since the Hamiltonian is a diagonally
dominant operator and the kinetic energy of PWs is the domi-
nant part at high G, this simple form is very effective.
2. Conjugate-Gradient
The eigenvalue problem of Eq.(8) can be recast into a se-
quence of constrained minimization problems:
min
_
_
?
i
|H|?
i
?
j?i
?
j
(?
i
|S|?
j
??
ij
)
_
_
, (14)
where the ?
j
are Lagrange multipliers. This can be solved
using preconditioned conjugate gradient algorithm with mi-
nor modi?cations to ensure constraint enforcement. The algo-
rithm here described was inspired by the conjugate-gradient
algorithm of Ref.[183], and is similar to one of the variants
described in Ref.[184].
Let us assume that eigenvectors ?
j
up to j = i ? 1 have
already been calculated. We start from an initial guess ?
(0)
for the i-th eigenvector, such that ?
(0)
|S|?
(0)
= 1 and
?
(0)
|S|?
j
= 0. We introduce a diagonal precondition ma-
trix P and auxiliary functions y = P
?1
? and solve the equiv-
alent problem
min
_
y|
¯
H|y ??
_
y|
¯
S|y ?1
__
, (15)
where
¯
H = PHP,
¯
S = PSP, under the additional orthonor-
mality constraints y|PS|?
j
= 0. The starting gradient of
Eq.(15)) is given by
g
(0)
= (
¯
H ??
¯
S)y
(0)
. (16)
By imposing that the gradient is orthonormal to the starting
vector: g
(0)
|
¯
S|y
(0)
= 0, one determines the value of the
Lagrange multiplier:
? =
y
(0)
|
¯
S
¯
H|y
(0)
y
(0)
|
¯
S
2
|y
(0)
. (17)
The remaining orthonormality constraints are imposed on
Pg
(0)
by explicit orthonormalization (e.g. Gram-Schmid) to
the ?
j
. We introduce the conjugate gradient h
(0)
, which for
14
the ?rst step is set equal to g
(0)
(after orthonormalization), and
the normalized direction n
(0)
= h
(0)
/h
(0)
|
¯
S|h
(0)
1/2
. We
search for the minimum of y
(1)
|
¯
H|y
(1)
along the direction
y
(1)
, de?ned as: [183]
y
(1)
= y
(0)
cos ? + n
(0)
sin?. (18)
This form ensures that the constraint on the norm is correctly
enforced. The calculation of the minimum can be analytically
performed and yields
? =
1
2
atan
_
a
(0)
(0)
?b
(0)
_
, (19)
where a
(0)
= 2y
(0)
|
¯
H|n
(0)
, b
(0)
= n
(0)
|
¯
H|n
(0)
, and
(0)
= y
(0)
|
¯
H|y
(0)
. The procedure is then iterated; at each
step the conjugate gradient is calculated from the gradient and
the conjugate gradient at the previous step, using the Polak-
Ribière formula:
h
= g
+ ?
(n?1)
h
(n?1)
, (20)
?
(n?1)
=
g
?g
(n?1)
|
¯
S|g
g
(n?1)
|
¯
S|g
(n?1)
. (21)
h
is subsequently re-orthogonalized to y
. We remarks
that in the practical implementation only Pg and Ph need to
be calculated and that only P
2
– the analogous of the D matrix
of Davidson algorithm – is actually used. A kinetic-only form
of P
2
has proved satisfactory:
k +G|P
2
|k +G
=
2m
¯h
2
(k +G)
2
?
GG
. (22)
C. Wavefunction extrapolation
In molecular dynamics runs and in structural relaxations,
extrapolations are employed to generate good initial guesses
for the wavefunctions at time t + dt from wavefunctions at
previous time steps. The extrapolation algorithms used are
similar to those described in Ref.[183]. The alignment pro-
cedure, needed when wavefunctions are the results of a self-
consistent calculation is as follows. The overlap matrix O
ij
between wavefunctions at consecutive time steps:
O
ij
= ?
i
(t + dt)|S(t + dt)|?
j
(t), (23)
can be used to generate the unitary transformation U [185]
that aligns ?(t +dt) to ?(t): ?
i
(t +dt) =
j
U
ij
?
j
(t +dt).
Since O is not unitary, it needs to be made unitary via e.g. the
unitarization procedure
U = (O
†
O)
?1/2
O
†
. (24)
The operation above is performed using a singular value de-
composition: Let the overlap matrix be O = vDw, where v
and w are unitary matrix and Dis a diagonal non-negative def-
inite matrix, whose eigenvalues are close to 1 if the two sets of
wavefunctions are very similar. The needed unitary transfor-
mation is then simply given by U w
†
v
†
. This procedure is
simpler than the original proposal and prevents the alignment
algorithm to break in the occasional situation where, due to
level crossing in the band structure between subsequent time
steps, one or more of the eigenvalues of the D matrix vanish.
D. Symmetry
Symmetry is exploited almost everywhere, with the notable
exception of CP. The latter is devised to study aperiodic sys-
tems or large supercells where symmetry is either absent or of
little use even if present.
In addition to lattice translations, the space group of a crys-
tal contains symmetry operations
ˆ
S combining rotations and
translations that leave the crystal unchanged:
ˆ
S ? {R|f },
where R is a 3 × 3 orthogonal matrix, f is a vector (called
fractional translation) and symmetry requires that any atomic
position, ?
s
is transformed into an equivalent one,
ˆ
S?
s
?
R(?
s
+ f ) = ?
ˆ
S(s)
. The rotational part of these operations
de?nes the crystal point group.
As a consequence of symmetry, roto-translated Kohn-Sham
orbitals are Kohn-Sham orbitals with the rotated Bloch vec-
tor:
ˆ
S?
i,k
(r) ? ?
i,k
(R
?1
r ?f ) = ?
i,Rk
(r). Where, strictly
speaking, the resulting wave-function at Rk does not neces-
sarily have the same band index as the original one but could
be some unitary transformation of states at Rk that share
with it the same single-particle eigenvalue. Since quantities
of physical interest are invariant for unitary rotations among
degenerate states this additional complication has no effect on
the ?nal result.
This is the basis for the symmetrization procedure used in
PWscf. One introduces a non-symmetrized charge density
(labeled by superscript
(ns)
) calculated on the irreducible BZ
(IBZ):
?
(ns)
(r) =
i
k?IBZ
w
k
|?
i,k
(r)|
2
. (25)
The factors w
k
(“weights”) are proportional to the number of
vectors in the star (i.e. inequivalent k vectors among all the
{Rk} vectors generated by the point-group rotations) and are
normalized to 1:
k?IBZ
w
k
= 1. Weights can either be
calculated or deduced from the literature on the special-point
technique[93, 94]. The charge density is then symmetrized as:
?(r) =
1
N
s
ˆ
S
ˆ
S?
(ns)
(r) =
1
N
s
ˆ
S
?
(ns)
(R
?1
r ?f ) (26)
where the sum runs over all N
s
symmetry operations.
The symmetrization technique can be extended to all quan-
tities that are expressed as sums over the BZ. Hellmann-
Feynman forces F
s
on atom s are thus calculated as follows:
F
s
=
1
N
s
ˆ
S
ˆ
SF
(ns)
s
=
1
N
s
ˆ
S
RF
(ns)
ˆ
S
?1
(s)
, (27)
15
where
ˆ
S
?1
(s) labels the atom into which the s?th atom trans-
forms (modulo a lattice translation vector) after application of
ˆ
S
?1
, the symmetry operation inverse of
ˆ
S. In a similar way
one determines the symmetrized stress, using the rule for ma-
trix transformation under a rotation:
?
??
=
1
N
s
ˆ
S
3
?,?=1
R
??
R
??
?
(ns)
??
. (28)
The PHonon package supplements the above technique
with a further strategy. Given the phonon wave-vector q, the
small group of q (the subgroup
ˆ
S
q
of crystal symmetry opera-
tions that leave q invariant) is identi?ed and the reducible rep-
resentation de?ned by the 3N
at
atomic displacements along
cartesian axis is decomposed into n
irr
irreducible representa-
tions (irreps) ?
(q)
j
, j = 1, . . . , n
irr
. The dimensions of the ir-
reducible representations are small, with ?
j
? 3 in most cases,
up to 6 in some special cases (zone-boundary wave-vectors q
in nonsymmorphic groups). Each irrep, j, is therefore de?ned
by a set of ?
j
linear combinations of atomic displacements
that transform into each other under the symmetry operations
of the small group of q. In the self-consistent solution of the
linear response equations, only perturbations associated to a
given irrep need to be treated together and different irreps can
be solved independently. This feature is exploited to reduce
the amount of memory required by the calculation and is suit-
able for coarse-grained parallelization and for execution on a
GRID infrastructure [186].
The wavefunction response, ??
(j,?)
k+q,i
(r), to displacements
along irrep j, ?
(q)
j,?
(where ? = 1, . . . , ?
j
labels different part-
ners of the given irrep), is then calculated. The lattice-periodic
unsymmetrized charge response, ??
(ns)
q,j,?
(r), has the form:
??
(ns)
q,j,?
(r) = e
?iq·r
4
i
k?IBZ(q)
w
k
?
?
k,i
(r)??
(j,?)
k+q,i
(r),
(29)
where the notation IBZ(q) indicates the IBZ calculated as-
suming the small group of q as symmetry group, and the
weights w
k
are calculated accordingly. The symmetrized
charge response is calculated as
??
q,j,?
(r) =
1
N
s
(q)
ˆ
S
q
e
?iqf
?
j
?=1
D(
ˆ
S
q
)
??
??
(ns)
q,j,?
(R
?1
r?f )
(30)
where D(
ˆ
S
q
) is the matrix representation of the action of the
symmetry operation
ˆ
S
q
? {R|f } for the j?th irrep ?
(q)
j
. At
the end of the self-consistent procedure, the force constant
matrix C
s?,t?
(q) (where s, t label atoms, ?, ? cartesian co-
ordinates) is calculated. Force constants at all vectors in the
star of q are then obtained using symmetry:
C
s?,t?
(Rq) =
?,?
R
??
R
??
C
ˆ
S
?1
(s)?,
ˆ
S
?1
(t)?
(q), (31)
where
ˆ
S ? {R|f } is a symmetry operation of the crystal
group but not of the small group of q.
E. Fock exchange
Hybrid functionals are characterized by the inclusion of a
fraction of exact (i.e. non-local) Fock exchange in the de?-
nition of the exchange-correlation functional. For a periodic
system, the Fock exchange energy per unit cell is given by:
E
x
= ?
e
2
N
kv
k
v
_
?
?
kv
(r)?
k
v
(r)?
?
k
v
(r
)?
kv
(r
)
|r ?r
|
drdr
,
(32)
where an insulating and non magnetic system is assumed for
simplicity. Integrals and wave-function normalizations are de-
?ned over the whole crystal volume, V = N? (? being the
unit cell volume), and the summations run over all occupied
bands and all N k-points de?ned in the BZ by Born-von Kár-
mán boundary conditions. The calculation of this term is per-
formed exploiting the dual-space formalism: auxiliary coden-
sities, ?k
,v
k,v
(r) = ?
?
k
,v
(r)?
k,v
(r) are computed in real space
and transformed to reciprocal space by FFT, where the associ-
ated electrostatic energies are accumulated. The application of
the Fock exchange operator to a wavefunction involves addi-
tional FFTs and real-space array multiplications. These basic
operations need to be repeated for all the occupied bands and
all the points in the BZ grid. For this reason the computational
cost of the exact exchange calculation is very high, at least an
order of magnitude larger than for non-hybrid functional cal-
culations.
In order to limit the computational cost, an auxiliary grid of
q-points in the BZ, centered at the ? point, can be introduced
and the summation over k
be limited to the subset k
= k+q.
Of course convergence with respect to this additional param-
eter needs to be checked, but often a grid coarser than the one
used for computing densities and potentials is suf?cient.
The direct evaluation of the Fock energy on regular grids in
the BZ is however problematic due to an integrable divergence
that appears in the q ?0 limit. This problem is addressed re-
sorting to a procedure, ?rst proposed by Gygi and Baldereschi
[187], where an integrable term that displays the same diver-
gence is subtracted from the expression for the exchange en-
ergy and its analytic integral over the BZ is separately added
back to it. Some care must still be paid [175] in order to es-
timate the contribution of the q = 0 term in the sum, which
contains a 0/0 limit that cannot be calculated from informa-
tion at q = 0 only. This term is estimated [175] assuming that
the grid of q-points used for evaluating the exchange integrals
is dense enough that a coarser grid, including only every sec-
ond point in each direction, would also be equally accurate.
Since the limiting term contributes to the integral with differ-
ent weights in the two grids, one can extract its value from the
condition that the two integral give the same result. This pro-
cedure removes an error proportional to the inverse of the unit
cell volume ? that would otherwise appear if this term were
simply neglected.
16
[1] N. Marzari, MRS BULLETIN 31, 681 (2006).
[2] R. G. Parr and W. Yang, Density Functional Theory of Atoms
and Molecules (Oxford University Press, 1989).
[3] R. M. Dreizler and E. K. U. Gross, Density Functional Theory
(Springer-Verlag, 1990).
[4] R. Martin, Electronic Structure: Basic Theory and Practi-
cal Methods (Cambridge University Press, Cambridge, UK,
2004).
[5] E. J. Baerends, J. Autschbach, A. Bérces, F. Bickelhaupt,
C. Bo, P. M. Boerrigter, L. Cavallo, D. P. Chong, L. Deng,
R. M. Dickson, et al., ADF2008.01, SCM, Theoretical Chem-
istry, Vrije Universiteit, Amsterdam, The Netherlands, URL
http://www.scm.com.
[6] R. Dovesi, B. Civalleri, R. Orlando, C. Roetti, and V. R. Saun-
ders, Ab initio quantum simulation in solid state chemistry, in
Reviews in Computational Chemistry, Ch. 1, Vol. 21, K. B.
Lipkowitz, R. Larter, T. R. Cundari editors, John Wiley and
Sons, New York (2005), URL http://www.crystal.
unito.it.
[7] G. Karlström, R. Lindh, P.-A. Malmqvist, B. O. Roos,
U. Ryde, V. Veryazov, P.-O. Widmark, M. Cossi, B. Schim-
melpfennig, P. Neogrady, et al., Comp. Mater. Sci. 28, 222
(2003), URL http://www.teokem.lu.se/molcas.
[8] R. Ahlrichs, F. Furche, C. Hättig, W. M. Klopper, M. Sierka,
and F. Weigend, TURBOMOLE, URL http://www.
turbomole-gmbh.com.
[9] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuse-
ria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr.,
T. Vreven, K. N. Kudin, J. C. Burant, et al., Gaussian 03,
Revision C.02, Gaussian, Inc., Wallingford, CT, 2004, URL
http://www.gaussian.com.
[10] Materials Studio, URL http://accelrys.com/
products/materials-studio.
[11] Jaguar: rapid ab initio electronic structure package,
Schrödinger, LLC., URL http://www.schrodinger.
com/.
[12] Spartan, Wavefunction, Inc., URL http://www.
wavefun.com/products/spartan.html.
[13] G. Kresse and J. Furthmuller, Comp. Mater. Sci. 6, 15 (1996),
URL http://cmp.univie.ac.at/vasp.
[14] URL http://www.castep.org.
[15] The CPMD consortium page, coordinated by M. Parrinello
and W. Andreoni, Copyright IBM Corp 1990-2008, Copy-
right MPI für Festkörperforschung Stuttgart 1997-2001, URL
http://www.cpmd.org.
[16] E. J. Bylaska, W. A. de Jong, N. Govind, K. Kowalski,
T. P. Straatsma, M. Valiev, D. Wang, E. Apra, T. L. Win-
dus, J. Hammond, et al., Nwchem, a computational chem-
istry package for parallel computers, version 5.1, Paci?c
Northwest National Laboratory, Richland, Washington 99352-
0999, USA, URL http://www.emsl.pnl.gov/docs/
nwchem/nwchem.html.
[17] R. A. Kendall, E. Apra, D. E. Bernholdt, E. J. Bylaska,
M. Dupuis, G. I. Fann, R. J. Harrison, J. Ju, J. A. Nichols,
J. Nieplocha, et al., Comput. Phys. Commun. 128, 260 (2000).
[18] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert,
M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A.
Nguyen, S. J. Su, et al., J. Comput. Chem. 14, 1347 (1993).
[19] M. S. Gordon and M. W. Schmidt, Advances in elec-
tronic structure theory: Gamess a decade later, in The-
ory and Applications of Computational Chemistry, the ?rst
forty years, C. E. Dykstra, G. Frenking, K. S. Kim, G.
E. Scuseria editors, ch. 41, pp 1167-1189, Elsevier, Ams-
terdam, 2005, URL http://www.msg.chem.iastate.
edu/GAMESS/GAMESS.html.
[20] URL http://www.gnu.org/licenses/.
[21] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs,
G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jol-
let, et al., Comp. Mater. Sci. 25, 478 (2002), URL http:
//www.abinit.org.
[22] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello,
T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103
(2005), URL http://cp2k.berlios.de.
[23] S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng. 4, 56
(2002), URL http://www.fysik.dtu.dk/campos.
[24] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys.
Rev. B 71 (2005), URL https://wiki.fysik.dtu.
dk/gpaw.
[25] T. D. Crawford, C. D. Sherrill, E. F. Valeev, J. T. Fermann,
R. A. King, M. L. Leininger, S. T. Brown, C. L. Janssen, E. T.
Seidl, J. P. Kenny, et al., J. Comput. Chem. 28, 1610 (2007),
URL http://www.psicode.org.
[26] A. A. Mosto?, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt,
and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
[27] URL http://qe-forge.org.
[28] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel,
J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling,
A. McKenney, et al., LAPACK Users’ Guide (Society for In-
dustrial and Applied Mathematics, Philadelphia, PA, 1999),
3rd ed., ISBN 0-89871-447-8 (paperback).
[29] M. Frigo and S. G. Johnson, Proceedings of the IEEE 93, 216
(2005), special issue on "Program Generation, Optimization,
and Platform Adaptation".
[30] M. P. I. Forum, Int. J. Supercomputer Applications 8 (3/4)
(1994).
[31] URL http://en.wikipedia.org/wiki/
SourceForge.
[32] W. E. Pickett, Comput. Phys. Rep, 9, 115 (1989).
[33] I. Dabo, B. Kozinsky, N. E. Singh-Miller, and N. Marzari,
Phys. Rev. B 77, 115139 (2008).
[34] L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425
(1982).
[35] D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett.
43, 1494 (1979).
[36] D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
[37] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
[38] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, D. J.
Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992).
[39] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys.
Rev. Lett. 91, 146401 (2003).
[40] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105,
9982 (1996).
[41] A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
[42] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J.
Frisch, J. Phys. Chem. 98, 11623 (1994).
[43] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[44] H. Hellmann, Einführung in die Quantenchemie (Deutliche,
1937).
[45] R. P. Feynman, Phys. Rev. 56, 340 (1939).
[46] O. H. Nielsen and R. M. Martin, Phys. Rev. B 32, 3780 (1985).
[47] P. Pyykkö, Chem. Rev. 88, 563 (1988).
[48] T. Oda, A. Pasquarello, and R. Car, Phys. Rev. Lett. 80, 3622
17
(1998).
[49] R. Gebauer and S. Baroni, Phys. Rev. B 61, R6459 (2000).
[50] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
[51] R. M. Wentzcovitch and J. L. Martins, Sol. St. Commun. 78,
831 (1991).
[52] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196
(1980).
[53] R. M. Wentzcovitch, Phys. Rev. B 44, 2358 (1991).
[54] S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58,
1861 (1987).
[55] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi,
Rev. Mod. Phys. 73, 515 (2001).
[56] X. Gonze, Phys. Rev. A 52, 1096 (1995).
[57] G. Henkelman and H. Jónsson, J. Chem. Phys. 111, 7010
(1999).
[58] W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66,
052301 (2002).
[59] K. J. Caspersen and E. A. Carter, P. Natl. Acad. Sci. USA 102,
6738 (2005).
[60] H. J. Choi and J. Ihm, Phys. Rev. B 59, 2267 (1999).
[61] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
[62] I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65,
035109 (2001).
[63] C. J. Pickard and F. Mauri, Phys. Rev. B 63, 245101 (2001).
[64] C. J. Pickard and F. Mauri, Phys. Rev. Lett. 91, 196401 (2003).
[65] M. Taillefumier, D. Cabaret, A. M. Flank, and F. Mauri, Phys.
Rev. B 66, 195107 (2002).
[66] S. Scandolo, P. Giannozzi, C. Cavazzoni, S. de Gironcoli,
A. Pasquarello, and S. Baroni, Z. Kristallogr. 220, 574 (2005).
[67] J. Tersoff and D. R. Hamann, Phys. Rev. B 31, 805 (1985).
[68] A. D. Becke and K. E. Edgecombe, J. Chem. Phys. 92, 5397
(1990).
[69] A. Szabo and N. Ostlund, Modern Quantum Chemistry
(Dover, 1996).
[70] A. Baldereschi, S. Baroni, and R. Resta, Phys. Rev. Lett. 61,
734 (1988).
[71] T. H. Group, URL http://www.hdfgroup.org.
[72] URL http://www.unidata.ucar.edu/software/
netcdf.
[73] URL http://www.quantum-simulation.org.
[74] G. Bussi, URL http://www.s3.infm.it/iotk.
[75] M. Fuchs and M. Schef?er, Comp. Phys. Comm. 119, 67
(1999).
[76] URL http://www.physics.rutgers.edu/~dhv/
uspp/.
[77] OPIUM – pseudopotential generation project, URL http:
//opium.sourceforge.net.
[78] P. Giannozzi and C. Cavazzoni, Nuovo Cimento C (in press,
2009).
[79] J. Hutter and A. Curioni, ChemPhysChem 6, 1788 (2005).
[80] L. S. Blackford, J. Choi, A. Cleary, J. Demmel, I. Dhillon,
J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley,
et al., SC Conference 0, 5 (1996).
[81] F. Gygi, E. W. Draeger, M. Schulz, B. R. de Supin-
ski, J. A. Gunnels, V. Austel, J. C. Sexton, F. Franchetti,
S. Kral, C. W. Ueberhuber, et al., IBM J. Res. Dev.
52, 137 (2008), URL http://eslab.ucdavis.edu/
software/qbox/index.htm.
[82] P. Giannozzi, F. de Angelis, and R. Car, J. Chem. Phys. 120,
5903 (2004).
[83] A. Dal Corso, A pseudopotential plane waves program
(pwscf) and some case studies, in Lecture Notes in Chem-
istry, Vol. 67, C. Pisani editor, Springer Verlag, Berlin (1996).
[84] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
[85] K. Laasonen, A. Pasquarello, R. Car, C. Lee, and D. Vander-
bilt, Phys. Rev. B 47, 10142 (1993).
[86] J. A. White and D. M. Bird, Phys. Rev. B 50, 4954 (1994).
[87] A. Dal Corso and A. Mosca Conte, Phys. Rev. B 71, 115106
(2005).
[88] A. Mosca Conte, Ph.D. thesis, SISSA/ISAS, Trieste Italy
(2007), URL http://www.sissa.it/cm/thesis/
2007/moscaconte.pdf.
[89] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B
44, 943 (1991).
[90] M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71, 035105
(2005).
[91] H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100,
126404 (2008).
[92] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
[93] D. J. Chadi and M. L. Cohen, Phys. Rev. B 8, 5747 (1973).
[94] H. J. Monkhorst and J. Pack, Phys. Rev. B 13, 5188 (1976).
[95] M. Methfessel and A. Paxton, Phys. Rev. B 40, 3616 (1989).
[96] N. Marzari, D. Vanderbilt, A. de Vita, and M. C. Payne, Phys.
Rev. Lett. 82, 3296 (1999).
[97] P. E. Blöchl, O. Jepsen, and O. Andersen, Phys. Rev. B 49,
16223 (1994).
[98] N. D. Mermin, Phys. Rev. 137, A1441 (1965).
[99] R. Fletcher, Practical Methods of Optimization (John Wiley
and Sons, 1987).
[100] S. R. Billeter, A. J. Turner, and W. Thiel, Phys. Chem. Chem.
Phys. 2, 2177 (2000).
[101] S. R. Billeter, A. Curioni, and W. Andreoni, Comp. Mater. Sci.
27, 437 (2003).
[102] C. Micheletti, A. Laio, and M. Parrinello, Phys. Rev. Lett. 92,
170601 (2004).
[103] L. Verlet, Phys. Rev. 159, 98 (1967).
[104] M. P. Allen and D. J. Tildesley, Computer Simulations of Liq-
uids (Clarendon Press, 1986).
[105] M. Bernasconi, G. L. Chiarotti, P. Focher, S. Scandolo,
E. Tosatti, and M. Parrinello, J. Phys. Chem. Solids 56, 501
(1995).
[106] P. Umari and A. Pasquarello, Phys. Rev. Lett. 89, 157602
(2002).
[107] I. Souza, J. Íñiguez, and D. Vanderbilt, Phys. Rev. Lett. 89,
117602 (2002).
[108] K. Kunc and R. Resta, Phys. Rev. Lett. 51, 686 (1983).
[109] J. Tobik and A. Dal Corso, J. Chem. Phys. 120, 9934 (2004).
[110] D. A. Scherlis, J.-L. Fattebert, F. Gygi, M. Cococcioni, and
N. Marzari, J. Chem. Phys. 124, 074103 (2006).
[111] I. Dabo, E. Cances, Y. Li, and N. Marzari (2009), URL http:
//arxiv.org/abs/0901.0096v2.
[112] C. Cavazzoni and G. L. Chiarotti, Comput. Phys. Commun.
123, 56 (1999).
[113] A. Pasquarello, K. Laasonen, R. Car, C. Lee, and D. Vander-
bilt, Phys. Rev. Lett. 69, 1982 (1992).
[114] M. d’Avezac, M. Calandra, and F. Mauri, Phys. Rev. B 71,
205210 (2005).
[115] H.-L. Sit, M. Cococcioni, and N. Marzari, J. of Electroanalyt-
ical Chemistry 607, 107 (2007).
[116] H.-L. Sit, M. Cococcioni, and N. Marzari, Phys. Rev. Lett. 97,
028303 (2006).
[117] N. Marzari, D. Vanderbilt, and M. C. Payne, Phys. Rev. Lett.
79, 1337 (1997).
[118] F. Tassone, F. Mauri, and R. Car, Phys. Rev. B 50, 10561
(1994).
[119] D. J. Tobias, G. J. Martyna, and M. L. Klein, J. Phys. Chem.
97, 12959 (1993).
[120] G. J. Martyna, M. L. Klein, and M. Tuckerman, The Journal
18
of Chemical Physics 97, 2635 (1992).
[121] M. C. Payne, M. P. Teter, D. C. Allen, T. A. Arias, and J. D.
Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992).
[122] D. Marx and J. Hutter, Modern methods and algorithms of
quantum chemistry, John von Neumann Institute for Comput-
ing, FZ Jülich, pp. 301–449 (2000).
[123] V. Dubois, P. Umari, and A. Pasquarello, Chem. Phys. Lett.
390, 193 (2004).
[124] F. Giustino and A. Pasquarello, Phys. Rev. Lett. 95, 187402
(2005).
[125] P. Umari and A. Pasquarello, Diam. Relat. Mater. 14, 1255
(2005).
[126] P. Umari and A. Pasquarello, Phys. Rev. Lett. 95, 137401
(2005).
[127] L. Giacomazzi, P. Umari, and A. Pasquarello, Phys. Rev. Lett.
95, 075505 (2005).
[128] P. Umari and A. Pasquarello, Phys. Rev. Lett. 98, 176402
(2007).
[129] P. Giannozzi and S. Baroni, J. Chem. Phys. 100, 8537 (1994).
[130] X. Gonze, Phys. Rev. B 55, 10337 (1997).
[131] X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
[132] A. Dal Corso, A. Pasquarello, and A. Baldereschi, Phys. Rev.
B 56, R11369 (1997).
[133] A. Dal Corso, Phys. Rev. B 64, 235118 (2001).
[134] F. Favot and A. Dal Corso, Phys. Rev. B 60, 11427 (1999).
[135] A. Dal Corso and S. de Gironcoli, Phys. Rev. B62, 273 (2000).
[136] A. Dal Corso, submitted (2009).
[137] A. Dal Corso, Phys. Rev. B 76, 054308 (2007).
[138] F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B 76,
165108 (2007).
[139] M. Lazzeri and S. de Gironcoli, Phys. Rev. B 65, 245402
(2002).
[140] M. Lazzeri and F. Mauri, Phys. Rev. Lett. 90, 036401 (2003).
[141] M. Lazzeri and F. Mauri, Phys. Rev. B 68, 161101(R) (2003).
[142] D. D. Koelling and B. N. Harmon, J. Phys. C: Solid State Phys.
10, 3107 (1977).
[143] A. H. MacDonald and S. H. Vosko, J. Phys. C: Solid State
Phys. 12, 2977 (1979).
[144] A. K. Rajagopal and J. Callaway, Phys. Rev. B 7, 1912 (1973).
[145] N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
[146] A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos,
Phys. Rev. B 41, 1227 (1990).
[147] G. Kresse and J. Hafner, J. Phys.-Condens. Mat. 6, 8245
(1994).
[148] A. Smogunov, A. Dal Corso, and E. Tosatti, Phys. Rev. B 70,
045417 (2004).
[149] A. Dal Corso, A. Smogunov, and E. Tosatti, Phys. Rev. B 74,
045429 (2006).
[150] M. Profeta, F. Mauri, and C. J. Pickard, J. Am. Chem. Soc.
125, 541 (2003).
[151] C. G. van de Walle and P. E. Blöchl, Phys. Rev. B 47, 4244
(1993).
[152] F. Mauri and S. G. Louie, Phys. Rev. Lett. 76, 4246 (1996).
[153] F. Mauri, B. Pfrommer, and S. G. Louie, Phys. Rev. Lett. 77,
5300 (1996).
[154] C. J. Pickard and F. Mauri, Phys. Rev. Lett. 88, 086403 (2002).
[155] T. Thonhauser, D. Ceresoli, A. Mosto?, N. Marzari, R. Resta,
and D. Vanderbilt, submitted to Phys. Rev. Lett. (2007), URL
http://arxiv.org/abs/0709.4429v1.
[156] C. Gougoussis, M. Calandra, A. Seitsonen, C. Brouder,
A. Shukla, and F. Mauri (2008), arXiv/0806.4706v1.
[157] URL http://www.wannier.org/.
[158] P. Silvestrelli, N. Marzari, D. Vanderbilt, and M. Parrinello,
Solid State Commun. 107, 7 (1998).
[159] A. Calzolari, I. Souza, N. Marzari, and M. B. Nardelli, Phys.
Rev. B 69, 035108 (2004).
[160] Y.-S. Leee, M. B. Nardelli, and N. Marzari, Phys. Rev. Lett.
95, 076804 (2005).
[161] J. R. Yates, X. Wang, D. Vanderbilt, and I. Souza, Phys. Rev.
B 75, 195121 (2007).
[162] A. Kokalj, J. Mol. Graph. Model. 17, 176 (1999), URL http:
//www.xcrysden.org/.
[163] W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph.
Model. 14, 33 (1996).
[164] URL http://www.tcm.phy.cam.ac.uk/~mdt26/
casino2.html.
[165] D. Korotin, A. V. Kozhevnikov, S. L. Skornyakov, I. Leonov,
N. Binggeli, V. I. Anisimov, and G. Trimarchi, Euro. Phys. J.
B 65, 91 (2008), URL http://[email protected].
[166] L. Martin-Samos and G. Bussi, J. Comp. Phys. Comm (2009),
URL http://sax-project.org.
[167] URL http://dp-code.org.
[168] URL http://www.bethe-salpeter.org.
[169] A. Kokalj, Comp. Mater. Sci. 28, 155 (2003), URL http:
//www-k3.ijs.si/kokalj/guib/.
[170] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
[171] M. A. L. Marques, C. L. Ullrich, F. Nogueira, A. Rubio,
K. Burke, and E. K. U. Gross, eds., Time-Dependent Den-
sity Functional Theory, vol. 706 of Lecture notes in Physics
(Springer-Verlag, Berlin, Heidelberg, 2006), DOI-10.1007/3-
540-35426-3-17.
[172] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601
(2002).
[173] D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, J. Chem. Phys.
128, 154105 (2008).
[174] P. Umari, G. Stenuit, and S. Baroni (2008), arXiv/0811.1453.
[175] H.-V. Nguyen and S. de Gironcoli, submitted to Phys. Rev. B
(2009), URL http://arxiv.org/abs/0902.0889.
[176] H.-V. Nguyen and S. de Gironcoli, submitted to Phys. Rev. B
(2009), URL http://arxiv.org/abs/0902.0883.
[177] H.-V. Nguyen, Ph.D. thesis, SISSA (2008), URL
http://www.sissa.it/cm/thesis/2008/
VietHuyNguyen_PhDthesis.pd%f.
[178] A. Marini, C. Hogan, M. Grüning, and D. Varsano, Comp.
Phys. Commun. p. in press (2009), URL http://www.
yambo-code.org.
[179] A. Calzolari, I. Souza, N. Marzari, and M. B. Nardelli,
Phys. Rev. B 69, 035108 (2004), URL http://www.
wannier-transport.org.
[180] URL http://[email protected].
[181] A. Kokalj, pwtk: a Tcl scripting interface to PWscf, soon to
be available on http://qe-forge.net.
[182] D. Raczkowski, A. Canning, and L. W. Wang, Phys. Rev. B
64, R121101 (2001).
[183] T. A. Arias, M. C. Payne, and J. D. Joannopoulos, Phys. Rev.
B 45, 1538 (1992).
[184] A. Qteish, Phys. Rev. B 52, 14497 (1995).
[185] C. A. Mead, Rev. Mod. Phys. 64, 51 (1992).
[186] R. di Meo, A. Dal Corso, P. Giannozzi, and S. Cozzini, Pro-
ceedings of the COST School, Trieste (2009), to be published.
[187] F. Gygi and A. Baldereschi, Phys. Rev. B 34, 4405 (1986).
doc_420872609.pdf