Home Search Collections Journals About Contact us My IOPscience

The SIESTA method; developments and applicability

This article has been downloaded from IOPscience. Please scroll down to see the full text article.

2008 J. Phys.: Condens. Matter 20 064208

(http://iopscience.iop.org/0953-8984/20/6/064208)

View the table of contents for this issue, or go to the journal homepage for more

Download details:
IP Address: 138.26.31.3
The article was downloaded on 12/12/2012 at 20:51

Please note that terms and conditions apply.

IOP PUBLISHING JOURNAL OF PHYSICS: CONDENSED MATTER
J. Phys.: Condens. Matter 20 (2008) 064208 (6pp) doi:10.1088/0953-8984/20/6/064208

The SIESTA method; developments and

applicability
Emilio Artacho1,2,10, E Anglada3 , O Diéguez4 , J D Gale5 , A Garcı́a4 ,
J Junquera6, R M Martin7 , P Ordejón8, J M Pruneda4,
D Sánchez-Portal9 and J M Soler3
1
Department of Chemistry, University of California, Berkeley, CA 94720, USA
2
Department of Earth Sciences, University of Cambridge, Cambridge CB2 3EQ, UK
3
Departamento de Fı́sica de la Materia Condensada, C-III, Universidad Autónoma,
28049 Madrid, Spain
4
Instituto de Ciencia de Materiales de Barcelona (ICMAB-CSIC), Campus de Bellaterra,
08193 Barcelona, Spain
5
Nanochemistry Research Institute, Curtin University of Technology, Perth, Western Australia
6845, Australia
6
Departamento de Ciencias de la Tierra y Fı́sica de la Materia Condensada, Universidad de
Cantabria, Avenida de Los Castros s/n, 39005 Santander, Spain
7
Department of Physics, University of Illinois, Urbana, IL 61801, USA
8
Centre d’Investigació en Nanociència i Nanotecnologia (CSIC-ICN), Campus de Bellaterra,
08193 Barcelona, Spain
9
Centro Mixto CSIC-UPV/EHU, Facultad de Quı́mica, Apartado 1072, 20080 San Sebastián,
Spain

Received 4 November 2007

Published 24 January 2008
Online at stacks.iop.org/JPhysCM/20/064208

Abstract
Recent developments in and around the SIESTA method of first-principles simulation of
condensed matter are described and reviewed, with emphasis on (i) the applicability of the
method for large and varied systems, (ii) efficient basis sets for the standards of accuracy of
density-functional methods, (iii) new implementations, and (iv) extensions beyond ground-state
calculations.

1. Introduction and overview methods were proposed and tested on empirical Hamiltonians,
the attention soon shifted to the problem of how to perform
Within the last two decades, the first-principles simulation full DFT calculations with a linear-scaling effort, since the
of condensed matter systems has expanded spectacularly, calculation of the Kohn–Sham Hamiltonian also displayed
from physics and chemistry into life, earth, nanoscale and inefficient scaling. The SIESTA method started in 1995 [2]
materials sciences. This success has been based on both the by merging Sankey’s finite-support atomic orbitals [3], with
steady growth of computing power and the development of a 3D real-space grid representation of the density and the basis
methods based on density-functional theory (DFT). However, functions. The grid offered a natural way of dealing with
deeply rooted algorithms within the standard methods demand long-range electrostatics, which does not allow for localized
computer resources that grow too rapidly (scaling as N 3 ) treatments (localization being the key for linear scaling). This
when increasing the number of atoms in a simulation box, N . method was coupled to several linear-scaling solvers, and
Although cube-scaling algorithms are still used in most of the coded in the SIESTA program in 1996 [4], a pioneer in linear-
DFT calculations done today (up to about 1000 atoms on the scaling DFT.
best supercomputers), it is obvious that in the long run linear- Other linear-scaling efforts started around that time in both
scaling algorithms will be of advantage. the quantum chemistry (QC) and condensed matter physics
The early 1990s witnessed considerable activity in communities. When overviewing the status of linear-scaling
the search for algorithms that could solve one-particle-like DFT, the first point to emphasize is that we are referring to
Hamiltonians in linear-scaling operations [1]. After several Kohn–Sham based DFT, as opposed to the so-called orbital-
10 Permanent address: Department of Earth Sciences, University of free DFT [5], which computes the total energy directly from
Cambridge, Cambridge CB2 3EQ, UK. the electron density with approximations to the kinetic energy

0953-8984/08/064208+06$30.00 1 © 2008 IOP Publishing Ltd Printed in the UK

J. Phys.: Condens. Matter 20 (2008) 064208 E Artacho et al

functional in addition to exchange and correlation. This Finally, there are hybrid methods that use a mix
approach is not only linear-scaling, but extremely efficient, of technologies from both the physics and chemistry
and best suited for metals, where linear-scaling Kohn–Sham communities: atomic-orbital basis sets and integration grids,
approaches still fail. The main weakness of orbital-free DFT is or plane waves as auxiliary basis set. This is the case of
that the kinetic energy functionals proposed so far offer good the CP2K program (formerly QUICKSTEP) [18], which uses the
approximations only for systems close to the homogeneous Gaussian bases of the chemists, and calculates some integrals
electron liquid, i.e. mostly simple metals [5]. as in QC, some using plane waves [19]. SIESTA as a method also
Within Kohn–Sham approaches, several groups in QC falls within this family, but using numerical atomic orbitals
have obtained and implemented methods for building the of finite support instead of Gaussians. This implies several
Hamiltonian in linear-scaling operations, using their traditional differences. Gaussian basis sets have been developed over
Gaussian basis sets, while the linear-scaling Hamiltonian decades of QC work, and are tabulated, and thus can be
solvers have not been widely introduced so far. Head- used mainly off the shelf. This is not the case for numerical
Gordon and collaborators started the line within QC [6] orbitals of finite support, which are generated by other means
and their methods have been implemented in the QCHEM (see below). On the other hand, numerical orbitals offer
program [7]. Challacombe’s contributions [8] are implemented flexibility in the choice of radial shapes of the basis functions
in the MondoSCF package [9], and the work of Scuseria and at no cost, while varying the radial shapes in Gaussian bases
co-workers [10] have been incorporated into the GAUSSIAN requires adding more Gaussian primitives, which increases the
program [11]. A characteristic of the QC approach to linear computing effort needed for some of the Hamiltonian matrix
scaling is that the key localization approximations are normally elements. The other main difference is in the fact that finite
done by ‘thresholding’, i.e. by neglecting the matrix elements support offers matrix elements that are exactly zero when their
with values below a given threshold, while other approaches support regions do not overlap, while Gaussians tails do not
impose localization regions for the relevant functions from strictly vanish, and matrix elements have to be neglected, as
the outset. Another important difference with respect to for the ‘thresholding’ described above. This difference is
tendencies in physics is that these methods obtain non- only important if enforcing tight localization: it is a valid
local exact exchange routinely, with an effort similar to the approximation in both cases, but in the finite-support case, the
computation of the Hartree term, given that they calculate both calculations remain exactly within the defined Hilbert space,
within the atomic-orbital representation by computing large while in the other case, the operators deviate from it, with
amounts of bi-electronic integrals over Gaussian functions. possible numerical instabilities for large thresholds.
They have thus linear-scaling access to Hartree–Fock and The SIESTA method has been evolving since its inception,
hybrid-functional Hamiltonians at the cost of possibly larger and has been applied by an ever expanding user community
scaling pre-factors for comparable basis sizes (quantitative to a large variety of problems in many different fields. Note
comparisons are beyond the scope of this work). that we are distinguishing here between the SIESTA method
Among physicists, the early proposal of finite-elements and its program implementation: the former is the set of
calculations (the so-called blip functions, distributed in a algorithms and ideas as published (see [20] for a detailed
Cartesian way around atoms) together with grid integrals [12] account), and thus what defines the calculations, while there
has been implemented in the CONQUEST code [13]. The are different implementations by other groups in addition to
scalability of this program for parallel computing is the SIESTA program developed by the SIESTA team. In this paper
outstanding. It offers an unbiased basis set that can be we report on the present status of the method and its main
improved systematically just by introducing more, smaller implementation, other developments that interface with it, and
blips homogeneously in space. The same advantage is its applicability.
offered by the much more recent development of the ONETEP
scheme [14]. It works on a real-space discretization in a similar
2. Recent developments in SIESTA
fashion as a finite-difference method would, but computes the
kinetic energy and other matrix elements with fast Fourier 2.1. Multigrid solver for electrostatics
transforms on different ‘boxes’, i.e. regions of space ample
enough to offer a good approximation to the matrix element The SIESTA method as described before [20] and used nowadays
at hand. A strictly finite-difference approach to linear-scaling relies on periodic boundary conditions (PBC) in the three
DFT has also been developed by Fattebert and co-workers [15], directions of space, and so, clusters, molecules or point defects
and there are also other developments happening around (0D), chains, tubes or line defects (1D), and surfaces, interfaces
finite-element methods [16]. The capacity for an unbiased or plane defects (slabs, 2D) are treated with 3D supercells. The
convergence of the basis set, and, not least, for having a clear only algorithm that imposes these boundary conditions is the
and unique procedure to ensure such convergence comes at a fast Fourier transform (FFT) used for the Hartree term, which
price, which is again the higher computational cost reflected in scales as N log N . To replace FFTs, a multigrid solver [21] for
a higher pre-factor to the linear scaling. The pre-factor scales the Poisson equation on the grid has been implemented [22].
with the number of basis functions per atom, with powers that It allows the SIESTA method to become strictly linear-scaling,
depend on the particular method and implementation. Wavelets and offers substantial flexibility with the boundary conditions,
provide basis sets with localization in both real and reciprocal allowing for truly 0D, 1D, 2D and 3D calculations. Dirichlet
space, and are also used in this context [17]. boundary conditions are used on single clusters or molecules

2
J. Phys.: Condens. Matter 20 (2008) 064208 E Artacho et al

1.41
Dipole (Debye)

1.40

1.39

1.38
FFT
1.37 MG (Dirichlet BC)

0
Energy (meV)

-2

-4

0 10 20 30
Box edge length (A)
Figure 2. Top panel: eggbox effect (energy variation with position)
Figure 1. Convergence with box size of energy and dipole moment for a Mo atom on a 70 Ryd grid. Lower panels: eggbox amplitude
for a water molecule in a cubic box, using FFT or multigrid (MG) as versus grid cut-off for I, B, and Os, with (dashed line) and without
the Poisson solver. filtering (continuous line).
(This figure is in colour only in the electronic version)

than finite-difference methods since some key matrix-element

(0D), whereby the cluster is not repeated, but the electrostatic calculations do not use the grid, a filtering of the neutral-atom
potential at the box boundaries is calculated according to the potential, the partial-core density, and the basis functions has
lower moments of the charge density of the system. The been proposed [24]. It preserves strict confinement in real
multigrid method is compatible with 3D PBC as well as hybrid space while maximizing the confinement in reciprocal space
BC (open in some directions, periodic in others), allowing for (a goal shared with wavelet methods). Figure 2 shows the size
of the eggbox versus grid cut-off for different atoms, with and
unrepeated tubes and slabs.
without filtering. The effect is dramatic. Although the filtered
For a molecule or cluster, the absence of spurious
basis functions are slightly different from the original ones,
interactions between replicas gives rapid convergence with box
they offer better energy conservation in ab initio molecular
size, substantially reducing the computational cost (figure 1).
dynamics (AIMD), more efficient relaxations, and much
In addition, charged clusters or molecules can be treated
smoother energy surfaces for finite-difference derivatives.
exactly, without compensating backgrounds. The present
multigrid implementation in its PBC form is around five times
slower than FFTs for small box sizes and taking the same 2.3. Parallelization
box for both, although the difference tends to disappear and The parallelization of the SIESTA code means parallelizing
even revert with size. The smaller sizes needed by multigrid (i) linear-scaling solvers, (ii) 3D grid operations, and
compensates for the difference in many cases. However, the (iii) diagonalization (maintained as an option for those
computational effort devoted to this is not an issue since either situations not amenable to linear scaling, such as metals
of them represents a minute amount (∼1%) of the global effort. or spectral studies). The first two parts are parallelized
using a real-space domain decomposition, which allows good
2.2. Fourier filtering onto the 3D grid scalability. Benchmarks with a linear-scaling solver [25] for
crystalline Si and a minimal basis went up to 524 288 atoms in
Grid methods are known to suffer from space rippling 256 nodes of an SGI Altix machine in 2001 [26].
problems: the homogeneity of space is lost by the The diagonalization does not scale as well. It is
discretization, and spurious forces appear. This is most parallelized using ScaLAPACK [27] based on a block cyclic
apparent in the oscillation of the energy when a single atom distribution of orbitals over nodes [26]. A substantial
is moved across the simulation box (‘eggbox’ effect, figure 2). improvement was obtained by Hein in the British HPCx
The problem diminishes for finer grids, but convergence to supercomputer centre, by using a 2D distribution of orbitals
acceptable levels can be slow. A significant improvement is instead of the original 1D scheme [28]. The great
obtained by grid-cell sampling [20], by averaging over shifted efficiency of the 2D implementation of ScaLAPACK more
positions of the grid points. than compensates for the extra manipulation it requires. A
In strictly grid based (finite-difference) methods, a realistic calculation of a pilin protein of 944 atoms, with
very good solution is provided by Fourier filtering the diagonalization and a good double-ζ polarized basis, shows
pseudopotentials to the plane-wave cut-off corresponding to speed-up factors 2, 1.82, 1.67 and 1.53, when going from
the grid [23]. Although SIESTA has a smaller rippling 8 to 16, 16 to 32, 32 to 64, and 64 to 128 processors,

3
J. Phys.: Condens. Matter 20 (2008) 064208 E Artacho et al

respectively [29], on the Cambridge HPCF supercomputer, of this method: TRANSIESTA, as a utility of SIESTA, TRANSIESTAC,
where an average time step (with 10 SCF steps) takes 25 min in distributed commercially by other authors, SMEAGOL [39], and
32 processors. In the absence of high-quality communications TRANSAMPA [40]. Recent developments include the possibility
the practical concurrency limit is 8–16 processors. A revision of inelastic scattering of electrons by phonons [41], as well as
of the parallelization strategies is being carried out at the the calculation of electromigration effects [42]. It should be
Barcelona Supercomputer Centre [30], in which, in addition remembered, however, that the Kohn–Sham spectrum, with its
to an optimization of domain decompositions, instances of known limitations, is behind the electron transport properties
global communications are being minimized and substituted obtained. Otherwise it overcomes many of the limitations of
by the minimum set of point-to-point communication calls, as previous approaches to transport.
obtained from graph-theory analysis.
3.3. Time-dependent DFT
3. Developments around SIESTA
Improved electronic spectra within the SIESTA context are
3.1. Phonons, polarization, effective charges and infrared obtained within the framework of time-dependent DFT (TD-
spectra DFT) in the frequency domain [43]. It takes advantage of
the localization of the response functions in real space, as
Normal modes of vibration are calculated by obtaining the discussed in [44]. TD-DFT was also implemented in the time
dynamical matrix with a finite-difference derivation of the domain [45], and was used initially to obtain the non-linear
forces [31], automatized in the VIBRA utility. Phonons at any optical response of C60 . It has been recently used to study the
point in the BZ can be computed using suitable supercells. This non-adiabatic response of insulators to ions moving through
is instead of the more popular way of calculating phonons, them at high velocity [46], the first ab initio calculation of
which uses density-functional perturbation theory (linear electronic stopping power in insulators. Time-evolving TD-
response) [32]. The latter represents a more elegant method DFT has recently been merged with the new treatment of finite
which can treat any phonon independently from the others. electric fields in periodic systems [47] to study the response of
The finite-difference approach scales better with system size, dielectric media to time-dependent electric fields.
however, allowing for linear-scaling calculations of the phonon
modes, including the electronic computation [33]. Linear-
response DFT has also been implemented on SIESTA for 4. Basis sets
the calculation of phonons [34]. There is no fundamental
The SIESTA method requires the use of finite-support atomic-
difference in the accuracy achievable by both methods, which
like basis functions, i.e., functions composed of a spherical
is limited by the underlying DFT theory, but finite-difference
harmonic and a radial function that becomes zero beyond a
calculations should be done with enough precision to get
radius rc . The user has absolute freedom in everything else:
sensible second derivatives: the eggbox should be minimized
where to centre them (on atomic nuclei or not), how many
as much as possible, and the structure used as reference should
angular momentum channels around any given centre, how
be at the energy minimum within a much tighter tolerance than
many radial functions for a given channel, what rc and what
that used for a structural study.
radial shape to use for each basis orbital. Wisely chosen basis
The infrared activity of the modes is obtained from the
Born effective charges (BEC; derivatives of the dielectric sizes and shapes optimize the efficiency–accuracy dichotomy.
polarization with respect to atom displacements). They are also There is a vast literature on this in the QC community, but
obtained from finite differences by calculating the polarization given the finite-support constraint and the numerical flexibility
at every displacement used in the force calculations for in the radial shape, the SIESTA community builds their own
the dynamical matrix. The BECs allow the calculation of bases following QC strategies. The SIESTA code incorporates
the splitting of the LO and TO phonon bands [35]. The different ways to introduce basis functions, and also offers
calculation of the polarization [36] uses by default the Berry- preset algorithms and criteria to define reasonable basis sets
phase formalism, but it becomes very demanding for large automatically [48].
systems. An efficient alternative uses the centres of the More accurate finite-support bases can be obtained
localized solution functions obtained by the linear-scaling variationally [49, 50]. Although a systematic prescription
solver of Ordejón et al [25], using the fact that they are has not been found, some practical rules have emerged from
Wannier-like functions [37]. experience: (i) the rc values for the first-ζ orbitals should not
be shorter than 5 bohr, and there is normally no point in their
being longer than 7 bohr. The energy-shift criterion [48] gives
3.2. Ballistic transport
too short values for orbitals of anions and for internal orbitals
Keldysh’s method of calculating ballistic currents based on (e.g. in partially filled 3d shells) but unnecessary large for the
non-equilibrium Green’s functions has been implemented on valence orbitals of alkalis. A value of 6.5 bohr is sensible in
top of the SIESTA method [38], exploiting the suitability of local general. (ii) Smooth is preferable to hard confinement [49],
bases for transport calculations. It allows for finite currents with an onset between 80% and 90% of rc , and V0 between
through junctions or interfaces, beyond the linear regime, by 50 and 100 Ryd. Their variations are not critical, except for
obtaining self-consistency in the presence of the bias voltage orbitals that are (close to) unbound in the free atom. (iii) The
and current. There are at least four different implementations quality of the basis is quite sensitive to the cut-off (matching)

4
J. Phys.: Condens. Matter 20 (2008) 064208 E Artacho et al

time step (10 SCF steps) taking 5 h in 128 processors [29].

8
Alternatively, hybrid QM/MM (quantum mechanics/molecular
mechanics) approaches have also been developed for SIESTA,
7 allowing for smaller quantum active sites embedded within a
r (a.u.)

much larger classical environment [61].

6
Complex oxides frequently appear in materials and earth
5 sciences, including transition-metal and rare-earth ions with
highly correlated electrons. Even the systems/properties
4
that can be reliably computed within present GGAs impose
computational demands significantly higher than for other
Na Mg Al Si P S Cl
materials (a relaxation of a slab of manganite became
Element
very heavy for only 224 atoms [62]). Many highly
Figure 3. Matching radii for the second-ζ of the 3s (filled symbols) correlated materials are not well described by GGAs. Other
and 3p (open symbols) orbitals across period III of the periodic table. functionals and corrections have been thus implemented,
Circles correspond to bases obtained variationally for Si(diamond), namely, LDA + U [63] and the pseudo-form of the self-
AlP (zinc-blende), MgS(NaCl), and NaCl. Squares refer to standard
bases with a 0.1 eV energy shift [48]. Percentile deviations of their interaction correction [64]. Exact exchange for Hartree–
respective lattice parameters are 0.0 (1.0), 0.7 (1.3), 0.2 (1.2), and Fock and hybrid functionals is being tested, although the
0.3 (0.0), for the variational (standard) bases, with respect to computational effort implied is substantially larger.
plane-wave results. Many studies demand a chemical analysis of their results.
Crystal orbital overlap and Hamiltonian populations [65]
have been recently implemented [66], as well as Voronoi
radii rm for second and higher ζ s. Figure 3 compares optimal deformation density charges, and Mulliken and Hirshfeld
values for rm with the ones obtained by default by SIESTA across populations [67]. Finally, of increasing importance for
the third period of the periodic table. The values towards the any method is the interoperability between different codes.
right are reasonably reproduced by the standard procedures, if SIESTA participates in various efforts, especially in data
slightly lower, but deviate strongly towards the left. rm values standardization. In addition to CML [68] (chemical markup
stay in the neighbourhood of 5 bohr. language) input/output [69], pseudopotential formats have
The preceding lines refer to varying parameters as in [49], been standardized [70] for the joint use of SIESTA and
which is a highly non-linear way of changing orbital shapes. AbINIT [71].
Ozaki [51] uses a linear procedure, by obtaining the basis
orbitals as linear combinations of a set of eigenstates of the 6. Conclusions and outlook
soft-confined atomic problem. An alternative is the expansion
of the desired orbitals in terms of the eigenstates of the filtering We have shown how SIESTA is developing, and how many new
projectors of [24]. This will be presented elsewhere. techniques are appearing around it exploiting its technology.
However, calculations using full linear scaling are still a small
5. Applicability fraction of the total, in spite of ten years of history. The main
reason for this is the wall encountered behind the N 3 barrier:
SIESTA is being used by around a thousand groups worldwide, the configurational complexity of the nuclear degrees of
on a great variety of systems [52, 53]. We give here just a freedom. It currently limits the systems that can be addressed
few illustrative examples. SIESTA has been used in nanoscience to sizes that rarely benefit from linear scaling. Regardless
since its inception [2], with many calculations of tubes [54], of foreseeable progress addressing nuclear complexity, much
wires, and clusters, as low connectivity increases SIESTA’s more efficient DFT methods are needed for complex systems.
efficiency. A very recent study [55] shows results for tapered We thus face the challenge of substantially reducing the linear-
Si wires with 4044 atoms, obtained with up to only four scaling pre-factor, while keeping standards of accuracy and
processors. robustness. Notwithstanding the possibility of new ideas, we
Although they demand high accuracy in addition to think there is scope for substantial progress in the exploration
their inherent size and complexity, biomolecules were treated of the many ideas already proposed, and of variants and
from the early days [56]. Hydrogen bonds are successfully synergies thereof.
reproduced and the dynamics of liquid water has been shown
to be amenable within the method [57]. For dispersion Acknowledgments
interactions a non-local functional [58] has been implemented
into SIESTA, by means of an algorithm based on a factorization We are grateful to O Paz, J Hein, J M Cela, S Riikonen,
of the integration kernel, and FFTs, well suited to the grid A Postnikov, T O White and many others for their many
structure in the method. Details on the implementation and contributions. JDG thanks the Government of Western
performance will be presented elsewhere [59]. Presently, short Australia for support under the Premier’s Research Fellowship
AIMD runs or relaxations of proteins of a few thousand atoms scheme, as well as that of iVEC and APAC. EA thanks
can be done. The pilin of Geobacter sulfurreducens [60] the Miller Institute for its support and the Department of
is being calculated in solution, with 4545 atoms, a typical Chemistry of the University of California at Berkeley for

5
J. Phys.: Condens. Matter 20 (2008) 064208 E Artacho et al

its hospitality. This work was partially supported by grant [37] Marzari N and Vanderbilt D 1997 Phys. Rev. B 56 12847
FIS2006-12117-C04 from Spain’s MEC. [38] Brandbyge M, Mozos J L, Ordejón P, Taylor J and
Stokbro K 2002 Phys. Rev. B 65 165401
[39] Rocha A R, Garcı́a-Suarez V, Bailey S W, Lambert C J,
References Ferrer J and Sanvito S 2005 Nat. Mater. 4 335
[40] Novaes F D, da Silva A J R and Fazzio A 2006 Braz. J. Phys.
[1] Goedecker S 1999 Rev. Mod. Phys. 71 1085 36 799
Ordejón P 2000 Phys. Status Solidi b 217 335 [41] Frederiksen T, Brandbyge M, Lorente N and Jauho A 2004
[2] Ordejón P, Artacho E and Soler J M 1996 Phys. Rev. B Phys. Rev. Lett. 93 256601
33 10441 [42] Brandbyge M, Stokbro K, Taylor J, Mozos J L and
[3] Sankey O F and Niklewski D J 1989 Phys. Rev. B 40 3979 Ordejón P 2003 Phys. Rev. B 67 193104
[4] Sánchez-Portal D, Ordejón P, Artacho E and Soler J M 1997 [43] Blase X and Ordejón P 2004 Phys. Rev. B 69 085111
Int. J. Quantum Chem. 65 453 [44] Blase X, Rubio A, Louie S G and Cohen M L 1995 Phys. Rev.
[5] Pearson M, Smargiassi E and Madden P A 1993 J. Phys.: B 52 R2225
Condens. Matter 5 3221 [45] Tsolakidis A, Sánchez-Portal D and Martin R M 2002 Phys.
[6] White C A and Head-Gordon M 1994 J. Chem. Phys. 101 6593 Rev. B 66 235416
[7] Shao Y et al 2006 Phys. Chem. Chem. Phys. 8 3172 [46] Pruneda J M, Sánchez-Portal D, Arnau A, Juaristi J I and
[8] Weber V and Challacombe M 2006 J. Chem. Phys.
Artacho E 2007 Phys. Rev. Lett. at press Preprint
125 104110 and references therein
arXiv:0706.1803
[9] http://www.t12.lanl.gov/home/mchalla/MondoSCF.html
[47] Souza I, Íñiguez J and Vanderbilt D 2004 Phys. Rev. B
[10] Strain M C, Scuseria G E and Frisch M J G 1996 Science
69 085106
271 51
[11] http://www.gaussian.com [48] Artacho E, Sánchez-Portal D, Ordejón P, Garcı́a A and
[12] Hernandez E and Gillan M J 1995 Phys. Rev. B 51 10157 Soler J M 1999 Phys. Status Solidi b 215 809
[13] Bowler D R, Choudhury R, Gillan M J and Miyazaki T 2006 [49] Junquera J, Paz O, Sánchez-Portal D and Artacho E 2001 Phys.
Phys. Status Solidi b 243 989 Rev. B 64 235111
[14] Skylaris C-K, Haynes P D, Mostofi A A and Payne M C 2005 [50] Anglada E, Soler J M, Junquera J and Artacho E 2002 Phys.
J. Chem. Phys. 122 084119 Rev. B 66 205101
[15] Fattebert J L and Gygi F 2006 Phys. Rev. B 73 115124 [51] Ozaki T 2003 Phys. Rev. B 67 155108
[16] Tsuchida E 2007 J. Phys. Soc. Japan 76 034708 [52] Sánchez-Portal D, Ordejón P and Canadell E 2004 Struct.
[17] Neelov A I and Goedecker S 2006 J. Comput. Phys. 217 312 Bond. 113 103
[18] VandeVondele J, Krack M, Mohamed F, Parrinello M, [53] http://www.uam.es/siesta
Chassaing T and Hutter J 2005 Comput. Phys. Commun. [54] Mazzoni M S C, Chacham H, Ordejón P, Sánchez-Portal D,
167 103 Soler J M and Artacho E 1999 Phys. Rev. B 60 R2208
[19] Lippert G, Hutter J and Parrinello M 1997 Mol. Phys. 92 477 [55] Wu Z, Neaton J B and Grossman J C 2007 submitted
[20] Soler J M, Artacho E, Gale J D, Garcı́a A, Junquera J, [56] Artacho E, Machado M, Sánchez-Portal D, Ordejón P and
Ordejón P and Sánchez-Portal D 2002 J. Phys.: Condens. Soler J M 2003 Mol. Phys. 101 1587
Matter 14 2745 [57] Fernandez-Serra M V and Artacho E 2004 J. Chem. Phys.
[21] Briggs W L, Henson V E and McCormick S F 2000 A 121 11136
Multigrid Tutorial (Norwood, MA: SIAM) [58] Dion M, Rydberg H, Schroder E, Langreth D C and
[22] Diéguez O and Ordejón P 2007 unpublished Lundqvist B I 2004 Phys. Rev. Lett. 92 246401
[23] Briggs E L, Sullivan D J and Bernholc J 1996 Phys. Rev. B Dion M, Rydberg H, Schroder E, Langreth D C and
54 14362 Lundqvist B I 2005 Phys. Rev. Lett. 95 109902
[24] Anglada E and Soler J M 2006 Phys. Rev. B 73 115122 [59] Soler J M et al 2007 in preparation
[25] Ordejón P, Drabold D A, Grumbach M P and Martin R M 1993 [60] Reguera G, McCarthy K D, Mehta T, Nicoll J S, Tuominen M T
Phys. Rev. B 48 14646 and Lovley D R 2005 Nature 435 1098
[26] Gale J D 2001 unpublished
[61] Crespo A, Scherlis D A, Marti M A, Ordejón P, Roitberg A E
[27] http://www.netlib.org/scalapack
and Estrin D A 2003 J. Phys. Chem. B 107 13728
[28] Hein J 2004 Technical Report of HPCx HPCxTR0410
[62] Pruneda J M, Ferrari V, Rurali R, Littlewood P B, Spaldin N A
[29] Feliciano G T, Reguera G and Artacho E 2007 in preparation
and Artacho E 2007 Preprint 0705.0335
[30] See http://www.bsc.es
[31] Ackland G J, Warren M C and Clark S J 1997 J. Phys.: [63] Riikonen S and Sánchez-Portal D 2007 PhD Thesis
Condens. Matter 9 7861 Universidad del Pais Vasco, Spain
[32] Baroni S, Gianozzi P and Testa A 1987 Phys. Rev. Lett. [64] Pemmaraju C D, Archer T, Sánchez-Portal D and
58 1861 Sanvito S 2007 Phys. Rev. B 75 045101
[33] Ordejón P, Drabold D A, Martin R M and Itoh S 1995 Phys. [65] Dronskowski R and Blochl P E 1993 J. Phys. Chem. 97 8617
Rev. Lett. 75 1324 [66] Llunell M, Garcı́a A, Alemany P and Canadell E 2007
[34] Pruneda J M, Estreicher S K, Junquera J, Ferrer J and unpublished
Ordejón P 2002 Phys. Rev. B 65 075210 [67] Guerra C F, Handgraaf J W, Baerends E J and Bickelhaupt F M
[35] Fernandez-Torre D, Escribano R, Archer T D, Pruneda J M and 2003 J. Comput. Chem. 25 189
Artacho E 2004 J. Phys. Chem. A 108 10535 and references [68] Guha R et al 2006 J. Chem. Inf. Mod. 46 991
therein [69] http://www.eminerals.org/tools/xml.html
[36] Sánchez-Portal D, Souza I and Martin R M 2000 Fundamental [70] Junquera J, Verstraete M, Gonze X and Garcı́a A 2006
Physics of Ferroelectrics 2000 (AIP Conf. Proc. vol 535) ed unpublished
R E Cohen pp 111–20 and references therein [71] Gonze X et al 2002 Comput. Mater. Sci. 25 478