This is the accepted version of the article: García A., Papior N., Akhtar A., Artacho E., Blum V., Bosoni E., Brandimarte P., Brandbyge M., Cerdá J.I., Corsetti F., Cuadrado R., Dikan V., Ferrer J., Gale J., García-Fernández P., García-Suárez V.M., García S., Huhs G., Illera S., Korytár R., Koval P.. Siesta: Recent developments and applications. The Journal of chemical physics, (2020). 152. : 204108 - . 10.1063/5.0005077. Available at: https://dx.doi.org/10.1063/5.0005077 S IESTA Siesta: recent developments and applications Alberto García,1, a) Nick Papior,2, b) Arsalan Akhtar,3, c) Emilio Artacho,4, 5, 6, 7, d) Volker Blum,8, 9, e) Emanuele Bosoni,1, f) Pedro Brandimarte,5, g) Mads Brandbyge,10, h) J. I. Cerdá,11, i) Fabiano Corsetti,4, j) Ramón Cuadrado,3, k) Vladimir Dikan,1, l) Jaime Ferrer,12, 13, m) Julian Gale,14, n) Pablo García-Fernández,15, o) V. M. García-Suárez,12, 13, p) Sandra García,3, q) Georg Huhs,16, r) Sergio Illera,3, s) Richard Korytár,17, t) Peter Koval,18, u) Irina Lebedeva,4, v) Lin Lin,19, 20, w) Pablo López-Tarifa,21, x) Sara G. Mayo,22, y) Stephan Mohr,16, z) Pablo Ordejón,3, aa) Andrei Postnikov,23, bb) Yann Pouillon,15, cc) Miguel Pruneda,3, dd) Roberto Robles,21, ee) Daniel Sánchez-Portal,21, 5, ff) Jose M. Soler,22, 24, gg) Rafi Ullah,4, 25, hh) Victor Wen-zhe Yu,8, ii) and Javier Junquera15, jj) 1) Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Bellaterra E-08193, Spain 2) DTU Computing Center, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark 3) Catalan Institute of Nanoscience and Nanotechnology - ICN2, CSIC and BIST, Campus UAB, 08193 Bellaterra, Spain 4) CIC Nanogune BRTA, Tolosa Hiribidea 76, 20018 San Sebastián, Spain 5) Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastian, Spain 6) Ikerbasque, Basque Foundation for Science, 48011 Bilbao, Spain 7) Theory of Condensed Matter, Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom 8) Department of Mechanical Engineering and Materials Science, Duke University, Durham, NC 27708, USA 9) Department of Chemistry, Duke University, Durham, NC 27708, USA 10) DTU Physics, Center for Nanostructured Graphene (CNG), Technical University of Denmark, Kgs. Lyngby, DK-2800, Denmark 11) Instituto de Ciencia de Materiales de Madrid ICMM-CSIC, Cantoblanco, 28049 Madrid, Spain 12) Department of Physics, University of Oviedo, Oviedo, 33007, Spain 13) Nanomaterials and Nanotechnology Research Center, CSIC - Universidad de Oviedo, Oviedo, 33007, Spain 14) Curtin Institute for Computation, Institute for Geoscience Research (TIGeR), School of Molecular and Life Sciences, Curtin University, PO Box U1987, Perth, WA 6845, Australia 15) Departamento de Ciencias de la Tierra y Física de la Materia Condensada, Universidad de Cantabria, Cantabria Campus Internacional, Avenida de los Castros s/n, 39005 Santander, Spain 16) Barcelona Supercomputing Center, c/ Jordi Girona, 29, 08034 Barcelona, Spain 17) Department of Condensed Matter Physics, Faculty of Mathematics and Physics, Charles University, Ke Karlovu 5, 121 16 Praha 2, Czech Republic 18) Simune Atomistics S.L., Tolosa Hiribidea, 76, 20018, Donostia-San Sebastian, Spain 19) Department of Mathematics, University of California, Berkeley, CA 94720, USA 20) Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA 21) Centro de Física de Materiales, Centro Mixto CSIC-UPV/EHU, Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastian, Spain 22) Departamento de Física de la Materia Condensada, Universidad Autónoma de Madrid, 28049 Madrid, Spain 23) LCP-A2MC, Université de Lorraine, 1 Bd Arago, F-57078 Metz, France 24) Instituto de Física de la Materia Condensada (IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain 25) Departamento de Física de Materiales, UPV/EHU, Paseo Manuel de Lardizabal 3, 20018 Donostia-San Sebastián, Spain (Dated: April 20, 2020. Accepted by Jour. of Chem. Phys. After publication it can be found at https://doi.org/10.1063/5.0005077) A review of the present status, recent enhancements, and applicability of the S IESTA program is presented. Since its debut in the mid-nineties, S IESTA’s flexibility, efficiency and free distribution has given advanced materials simulation capabilities to many groups worldwide. The core methodological scheme of S IESTA combines finite-support pseudo- atomic orbitals as basis sets, norm-conserving pseudopotentials, and a real-space grid for the representation of charge density and potentials and the computation of their associated matrix elements. Here we describe the more recent imple- mentations on top of that core scheme, which include: full spin-orbit interaction, non-repeated and multiple-contact bal- listic electron transport, DFT+U and hybrid functionals, time-dependent DFT, novel reduced-scaling solvers, density- functional perturbation theory, efficient Van der Waals non-local density functionals, and enhanced molecular-dynamics options. In addition, a substantial effort has been made in enhancing interoperability and interfacing with other codes and utilities, such as WANNIER 90 and the second-principles modelling it can be used for, an AiiDA plugin for workflow automatization, interface to Lua for steering S IESTA runs, and various postprocessing utilities. S IESTA has also been S IESTA 2 engaged in the Electronic Structure Library effort from its inception, which has allowed the sharing of various low level libraries, as well as data standards and support for them, in particular the PSML definition and library for transferable pseudopotentials, and the interface to the ELSI library of solvers. Code sharing is made easier by the new open-source licensing model of the program. This review also presents examples of application of the capabilities of the code, as well as a view of on-going and future developments. I. INTRODUCTION. As we shall see, the improvements touch many areas. We can underline the implementation of new core electronic- The possibility of treating large systems with first- structure features (DFT+U, spin-orbit interaction, hybrid principles electronic-structure methods has opened up new re- functionals), modes of operation (improved time-dependent search avenues in many disciplines. The S IESTA method and density functional theory (TD-DFT), density functional per- its implementation have been key in this development, offer- turbation theory (DFPT), and analysis methods and proce- ing an efficient and flexible simulation paradigm based on the dures to access new properties. A major effort has been spent use of strictly localized basis sets. This approach enables the in enhancing the interoperability of the code at various levels implementation of reduced scaling algorithms, and its accu- (sharing of pseudopotentials, a new wannierization interface racy and cost can be tuned in a wide range, from quick ex- opening the way to sophisticated post-processing, and an in- ploratory calculations to highly accurate simulations matching terface to multiscale methods). Very significant performance the quality of other approaches, such as plane-wave methods. enhancements have been made, notably to the T RAN S IESTA The S IESTA method has been described in detail in Ref. 1, module through improved algorithms, and to the core elec- with an update in Ref. 2. In this paper we shall describe its tronic structure problem through the development of inter- present status, highlighting its strengths and documenting the faces to new solvers. These advances have put S IESTA in a steps that have recently been taken to improve its capabili- prominent place in the high-performance electronic-structure ties, performance, ease of use, and visibility in the electronic- simulation scene, a role reinforced by its participation in im- structure community. portant international initiatives and by its new open-source li- censing model. The manuscript is organized as follows. We provide an overview of the underlying methodology and the capabilities a) Electronic mail: albertog@icmab.es b) Electronic of S IESTA in section II, which serves to place the code in the mail: nicpa@dtu.dk c) Electronic mail: arsalan.akhtar@icn2.cat wider ecosystem of electronic-structure materials simulation. d) Electronic mail: ea245@cam.ac.uk Section III presents the recent developments in and around e) Electronic mail: volker.blum@duke.edu the code, which are covered in sub-sections. To showcase f) Electronic mail: ebosoni@icmab.es S IESTA’s utility in the context of electronic-structure calcula- g) Electronic mail: pedro_brandimarte001@ehu.eus tions, we present briefly some relevant applications and survey h) Electronic mail: mabr@dtu.dk i) Electronic mail: jcerda@icmm.csic.es a few areas in which S IESTA is being profitably used in sec- j) Electronic mail: fabiano.corsetti@gmail.com tion IV. Plans for the future evolution of S IESTA are outlined k) Electronic mail: ramon.cuadrado@gmail.com in section V. l) Electronic mail: vdikan@icmab.es m) Electronic mail: ferrer@uniovi.es n) Electronic mail: J.Gale@curtin.edu.au o) Electronic mail: garciapa@unican.es II. KEY CONCEPTS OF SIESTA p) Electronic mail: vm.garciasuarez@gmail.com q) Electronic mail: sandragil@gmail.com A. Theory background and context r) Electronic mail: ghuhs@physik.hu-berlin.de s) Electronic mail: sergiollera22@gmail.com t) Electronic mail: korytar@karlov.mff.cuni.cz S IESTA appeared as a consequence of the push for linear- u) Electronic mail: koval.peter@gmail.com scaling electronic structure methods of the mid nineties, v) Electronic mail: i.lebedeva@nanogune.eu which has been reviewed, for example, in Refs. 3 and 4. w) Electronic mail: linlin@math.berkeley.edu x) Electronic mail: pablolopeztarifa@gmail.com S IESTA was the first linear-scaling self-consistent implemen- y) Electronic mail: sara.garciamayo@uam.es tation of density functional theory (DFT).5,6 z) Electronic mail: stephan.mohr@bsc.es The S IESTA method relies on atomic-like functions of fi- aa) Electronic mail: pablo.ordejon@icn2.cat nite support as basis sets7,8 – of arbitrary number, angular bb) Electronic mail: andrei.postnikov@univ-lorraine.fr momentum, radial shape, and centers – combined with a dis- cc) Electronic mail: yann.pouillon@unican.es dd) Electronic mail: miguel.pruneda@icn2.cat cretization of space for the computation of the Kohn-Sham ee) Electronic mail: roberto.robles@ehu.eus Hamiltonian terms that involve more than two centers. The ff) Electronic mail: daniel.sanchez@ehu.eus electron-ion interaction is represented by norm-conserving gg) Electronic mail: jose.soler@uam.es pseudopotentials. These key ingredients, through the opti- hh) Electronic mail: ullah1@llnl.gov mized handling of sparse matrices, are used to compute the ii) Electronic mail: wenzhe.yu@duke.edu self-consistent Hamiltonian and overlap matrices with a com- jj) Electronic mail: javier.junquera@unican.es putational expense that scales linearly with system size. The S IESTA 3 method is completed with a choice of solvers for that Hamilto- energies, forces, molecular-dynamics simulations, band struc- nian, from optimized (but cube-scaling) diagonalization meth- tures, densities of states, etc., and shares with those codes the ods, to reduced-scaling solvers of different flavors. basic current limitations of DFT (notably the description of The orbitals in the S IESTA basis set are made of the product strongly-correlated systems). of a real spherical harmonic and a radial function, which is nu- What makes S IESTA different from most other codes, and is merically tabulated in a grid. The shape of the radial part is in at the root of its key strengths, is the atomic-like, and strictly principle totally arbitrary, but the experience accumulated has localized, character of its basis set. The use of a “good first proven that the numerical solution of the Schrodinger equa- approximation” to the full problem implies, first, that a much tion for a (confined) isolated atom with the corresponding smaller number of basis functions is needed. Second, the pseudopotential is a very good choice in terms of accuracy finite-support of the orbitals leads to sparsity and the possibil- versus computational cost. Fuller descriptions of the mecha- ity to use reduced-scaling methods. Thus high performance nisms to generate and optimize these pseudo-atomic orbitals emerges almost by default. (PAOs) are given in Refs. 8–10. Take first the basis cardinality: the number of basis orbital The auxiliary real-space grid is an essential ingredient of per atom in a typical S IESTA calculation is of the order of the method, as it allows the efficient representation of charge 10-20. This is to be compared with a few hundred in the typ- densities and potentials, as well as the computation of the ma- ical plane-wave (PW) calculation. Furthermore, for systems trix elements of the Hamiltonian that cannot be handled as whose description needs a vacuum region (e.g., slabs for sur- two-center integrals. This grid can be seen as the reciprocal face calculations, 2D monolayers, etc), empty space is essen- space of a set of plane waves, and its fineness is most conve- tially “free” for S IESTA, whereas PW codes still need a basis niently parametrized by an energy cutoff (the “density” cutoff set determined by the total size of the simulation cell. S IESTA of plane-wave methods). There are limits to the softness of the is then quite capable of dealing with systems composed of functions that can be described with such a grid, so core elec- dozens to hundreds of atoms on modest hardware, even when trons are not considered (although semi-core electrons usually using cubic-scaling diagonalization solvers, which are the de- are), and their effect is incorporated into pseudopotentials. fault as they are universally applicable. The real-space grid is also used to solve the Poisson equation Electronic-structure solvers with a more favorable size- involved in the computation of the electrostatic potential from scaling can be applied to suitable systems. For example, one the charge density, through the use of a fast-fourier-transform of S IESTA’s earlier calculations, in 1996, was a linear-scaling method. This means that S IESTA uses periodic boundary con- run for a strand of DNA with 650 atoms, performed on a desk- ditions. Non periodic systems, such as molecules, tubes, or top workstation of the era.6 Reduced size-scaling is also a slabs, are treated using appropriate supercells. feature of the PEXSI solver described in section III G 1 be- S IESTA is now a mature code with more than 20 years of low, and of the NTPoly solver mentioned in section III G 2. existence. In this period, the most important algorithms be- In addition to time-to-solution efficiency, these solvers have hind our implementation have been already fully described a smaller memory footprint than diagonalization, as the rele- and documented in a series of papers. Readers interested in vant matrices are kept in sparse form rather than converted to the details of how the basic elements defining the method are a dense format. combined, as well as other relevant implementation details Crucially, S IESTA’s baseline efficiency can be scaled up to that make the method practical, can find them in the main ever-larger systems by parallelization. Both distributed (MPI) S IESTA reference1 , and in the update with the new capabil- and shared-memory (OpenMP) parallelization options are im- ities of the code2 . plemented in the code. As some of the examples in section IV We note that the term S IESTA is regularly used to describe show, non-trivial calculations with thousands of atoms are both the method (as outlined in the earliest papers5,6 ) and its used in applications in different contexts, from molecular bi- implementation in a computer program. The S IESTA method ology to electronic transport. is at the basis of later independent implementations, such Work on the performance aspects of the code is continu- as OpenMX,11 and QuantumATK.12 Other subsequent codes ous, mostly on the solvers, which usually take most of the built on the method, revising some of the fundamental ingre- computer time due to the very high efficiency of the Hamilto- dients. This is the case of FHI-aims,13 which uses a more nian setup module in S IESTA. This task is facilitated (see sec- sophisticated real-space grid (atom-centered), thus extending tion III O) by leveraging external libraries and developments the core scheme to all-electron calculations. generated by a number of international initiatives in which In this paper we describe new additions to the S IESTA code, S IESTA participates. The code can still run efficiently in mod- based on independent methodological advances, either pre- est hardware, while also being able to exploit massive levels existent or specifically developed for S IESTA, as specified and of parallelism in large supercomputers (see Fig. 1). cited in each section below. It is worth noting also that the atomic character of the basis set enables the use of a very intuitive suite of analysis tools, sice most of the concepts relating to chemical bonding use the B. Overview of Siesta capabilities language of atomic orbitals. Hence S IESTA has a natural ad- vantage in this area. Partial densities of states and atomic and As a general purpose implementation, S IESTA can provide crystal populations (COOP/COHP) are routinely used to gain the standard functionality available in mainstream DFT codes: insights into the stability and other properties of materials. for ScaLAPACK, Siesta-PEXSI needs only 144 cores for C-BN0.00 and 64 for DNA-25. In the case of DNA even this minimal configuration is more than four times faster proach Sthan IESTAdiagonalization with 5120 processors. 4 Strong scaling the base the community, shared by some of the most popular electronic cost with structure codes. ase of the Fourth, with regard to S IESTA specific approximations, par- trong scal- ticularly the basis set, it should be stressed that S IESTA is lim- alculating ited to basis sets composed of functions that are product of a Hamilto- 180,000 orbs radial part and spherical harmonics, but it does not constrain ion of the 170,000 orbs 2D, sp=0.91% on how many, where such functions are centered, and the size aLAPACK 1D, sp=0.27% of their finite-support region. Calculations can flexibly range h requires from quick exploration to very high-quality simulations (one ration. In may recall that accuracy gold standards in electronic struc- the chem- ture are provided by quantum-chemistry methods, based on number of , reducing LCAO). The use of an atomic-orbital basis set implies however the gh perfor- limitation of non-uniformity of convergence. As opposed to allelization plane-wave methods, in which a single energy cutoff param- ntly. The FIG.1.2. Parallel FIG. Strong strong scalingscaling of C-BNof 0.00 and DNA-25 SIESTA-PEXSI andbased on the (Scala- eter monotonically determines the quality of the calculation, uning the the total pack) time for the diagonalization first SCF approach for step. a DNAThe various chain and alines for Graphene- there is no univocal procedure for the choice of an appropri- number of PEXSINitride Boron result stack, from using di↵erent prototypes numbers of large of processors (hundreds per of of thousands ate basis set. It is a well-known problem, shared by the whole pole (ppp), orbitals) quasiwhile the points on one-dimensional andeach curve belongsystems. two-dimensional to compu- “ppp” is demon- quantum-chemistry community, on which there is widely used tationsforwith stands 1, 2, 5, 10, the number 20, and of MPI 40 poles processes usedininparallel. each pole computa- C-BN sys- and tested know-how. As Fig. 2 shows, it is possible to attain me ppp are tion, and “sp” the sparsity of the Hamiltonian. (For more details, see in practice an accuracy comparable to that of well-converged ling. The PEXSI’s section III G 1)beneficial scaling with the system size, as de- plane-wave calculations. The reader is also referred to sec- ation over scribed in section II C, guarantees that for large enough tions IV B and IV C 1 for showcase examples of the accuracy paralleliza- systems Siesta-PEXSI will always be faster than diag- of the code, among many others in the literature. come from onalization. The scaling of the computational cost is For a recent example, demonstrated for DNA seeand Ref.C-BN 14. Similarly, in figure 3. an atomic basis To close this section, we stress that it has been a tradi- only a lim- tional and deliberate attitude by the S IESTA team that, al- thus does provides In all tests full parallelization over poles is for a very natural and adequate language the In used. first- principles simulation of electronic ballistic this configuration the influence of the symbolic factoriza- transport in nano- though proposing sensible starting points to users as defaults, able. This currently sized systems, tion would via the change the Green’s-function character of the method. based Keldysh Becausefor- the choice of fundamental approximations and inputs to the e time for malism in futureimplemented this influence in T RAN will SbeIESTA ,15 a partthe negligible, of time the S IESTA for program (not only basis sets, but also density functionals, and nal role. package. symbolic factorization is not taken into account for the pseudopotentials) is a responsibility of the users, who retain e, demon- analysis. full control and the flexibility to adapt the code to their specific The very high number of citations of the S IESTA papers Thetonumbers of processes for each needs. Nevertheless, tools for basis optimization are provided es treated testify the successful application of system the codesize are cho-dif- to widely sen to be an efficient trade-o↵ of reducing the time with the program, new curated databases of pseudopotentials rther, but ferent systems. With regard to specific capabilities andtothe on. [LL: solution while keeping the cost, which increases with the are coming online, and new ways to ameliorate the correlation levels of accuracy achievable, we can distinguish several lev- the scaling number of processes due to inefficiencies, as small as pos- problem are being implemented. Some of these developments els. First, S IESTA implements DFT, one of the most versatile time by a sible. Following this guideline it turns out, that for C- are described in the following sections. materials BN one can simulation use moreframeworks. processors with DFTSiesta-PEXSI has its shortcomings, than both ex- notably with ScaLAPACK. This also means, that the advantagesys- in regard to the description of strongly-correlated hroughout tems, but these areinbeing of Siesta-PEXSI termsaddressed (see sections of solution-time is even onlarger DFT+U III. RECENT DEVELOPMENTS IN SIESTA of Siesta- and thanhybrid functionals the benefit of cost.below). For verySecond, sparseS IESTA problems, useslike pseu- PEXSI ap- dopotentials the largest DNA examples, the amount of processors thatThe to represent the electron-ion interaction. pseudopotential can be used is approach similar for is firmly rooted in abut both methods, sound physi- Siesta- A. New distribution model and development infrastructure allows an the scalig cal approximation PEXSI is about two (thatorders bonding of effects magnitude depend mostly faster. on the More this. For details electrons); valence are listed inhowever, table II.it is at a disadvantage when core- A few years ago, in 2016, a decision was made to change ce this ex- The analysis electrons effects areshows, besides(but important Siesta-PEXSI’s see section IIIfavorable N 6 below). the licensing model for S IESTA: traditionally it had always Third, S IESTA employs periodic boundary conditions (PBC) been free of charge to academics, but non-academic use re- for the solution of the Poisson problem, sharing with plane- quired a special license and redistribution was not permitted. wave codes the need to resort to repeated supercells for the Now S IESTA is formally an open-source program, distributed study of low-dimensional systems, and to special techniques according to the terms of the GPL license18 . At the same for the treatment of charged systems. It is important to note, time, the development infrastructure was made more trans- however, that, unlike plane-wave codes, S IESTA is only bound parent and scalable, using first the Launchpad platform19 and to PBC because of the present treatment of the Hartree term now the Gitlab service20 . The net effect of the changes has of the single-particle Hamiltonian. This limitation is lifted by been a more fertile and dynamic development, with more con- the incorporation of alternative Poisson solvers, as described tributors who can have direct access to the various branches in Sec. III O, which allow for open boundary conditions, as of development, and a better experience for users, who can for isolated nano-systems, and hybrid open/periodic bound- download code and raise issues in an integrated platform. ary conditions in different dimensions, as for isolated wires These changes have been substantial for the core develop- and slabs. It should be remembered that the three approxi- ers, and the transitory period is still being felt. The main code mations mentioned in this paragraph are very widely used in base is gradually absorbing new developments, both those that S IESTA 5 norm-conserving pseudopotentials, easing in most cases the 0.0 task of pseudopotential quality control. 250 (meV) -0.4 PW C. DFT+U for correlated systems Eb - Eb Eb (meV) -0.8 240 The LDA+U method, initially developed by Anisimov and 70 75 80 85 90 coworkers28 with the objective to improve the treatment of the Number of basis functions electron-electron interaction for localized electrons within the bare LDA description, has been implemented in S IESTA. The 230 PW idea behind the LDA+U consists in describing the “strongly correlated” electronic states of a system (typically, localized 30 50 70 90 d or f orbitals) using the Hubbard model, whereas the rest of Number of basis functions valence electrons are treated at the level of “standard” approx- imate DFT functionals.29 In the current version of S IESTA the FIG. 2. Basis set convergence for the binding energy (Eb ) of a water implementation is based on the simplified rotationally invari- dimer. Details can be found in Ref. 16. The horizontal dotted line ant functional proposed by Dudarev and coworkers.30 Here, represents the converged plane-wave (PW) calculation (1300 eV cut- the corrections are made invariant under rotation of the atomic off) for the same system (dimer geometry and box), pseudopotentials orbitals used to define the occupation number of the correlated and density functional, using the ABINIT code.17 Inset: deviation of subspace, at the cost of retaining only the lowest order Slater Eb versus the PW reference. The deviation for the last point is of 10 integrals in the factorization of the integrals of the Coulomb µeV. kernel of the electron-electron interaction, and neglecting the higher order ones (i.e. taking the exchange interaction J as 0). The expression of the corrective term as a functional of the were planned long in advance, and new ones made possible by occupation number nIσ `m of the localized correlated orbital `m the greater openness and fluidity of the development model. with spin σ within the atom I is given by Most of the new features described below are already part of public releases, but a few are undergoing the last stages of U I`   testing before release. The work-flow is also moving from EU = ∑ nIσ 1 − nIσ  , (1) ∑ `m `m long-lived releases, hard to maintain with bug-fixes, to more Iσ ` 2 m frequent releases that will be maintained for a shorter time. where only one interaction parameter U I` is needed to specify the interaction per atom and `-shell. In the practical S IESTA implementation, the populations on the correlated orbitals are B. New pseudopotential format for interoperability computed using non-overlapping (i. e. orthogonal) localized projectors. They can be generated using either (i) the same PSML (for PSeudopotential Markup Language)21,22 is a file algorithm used to produce the first-ζ orbitals of the basis set, format for norm-conserving pseudopotential data which is de- but with a larger energy shift, or (ii) cutting the exact solution signed to encapsulate as much as possible the abstract con- of the pseudoatom with a Fermi function. cepts involved in the domain, and to provide appropriate meta- The results of the LDA+U method are sensitively depen- data and provenance information. This extra level of formal- dent on the numerical value of the effective on-site electronic ization aims at removing the interoperability problems associ- interaction, the Hubbard U. Although in principle the value ated to bespoke pseudopotential formats, which usually were of U can be computed from first principles using linear re- designed to serve the needs of specific generators and client sponse methods,31 a common practice is to tune it semiempir- codes, and thus contain implicit assumptions about the mean- ically, seeking agreement of certain properties (for instance ing of the data or lack information not considered relevant. band gaps or magnetic moments) with available experimen- PSML files can be produced by the ONCVPSP23 and tal measurements. Then, the fitted U is used in subsequent ATOM24 pseudopotential generator programs, and are a calculations to predict other properties. download-format option in the Pseudo-Dojo database of cu- The LDA+U corrects localized states, for which the self- rated pseudopotentials25,26 . interaction correction is expected to be stronger, and is an ef- The software library libPSML21,22 can be used by elec- fective method to improve the description of the (underesti- tronic structure codes to transparently extract the information mated) band gap of insulators, as shown in Fig. 3 for the case in a PSML file and incorporate it into their own data struc- of NiO. Once the Hubbard correction is switched on, the op- tures, or to create converters for other formats. It is cur- tical band gap increases up to 3.08 eV (from the bare GGA- rently used by S IESTA and A BINIT,17,27 making possible a PBE value of 1.08 eV), very close to the experimental value full pseudopotential interoperability and facilitating compar- for the onset of optical absorption in NiO32 (3.10 eV). The isons of calculation results. magnetic moment on the Ni atom is also properly described, The use of this new format opens the door to benefit from with a value of 1.67 µB which lies well within the experi- the availability of a periodic table of reliable and accurate mental range of values (between 1.64 µB 33 and 1.9 µB 34 ), and S IESTA 6 improves on the result of 1.39 µB obtained with a bare GGA- E. Hybrid functionals PBE functional. The screened hybrid functional HSE0642–44 has been im- 5 plemented in S IESTA building on the work of Ref. 45. This functional is the result of adding nonlocal Hartree-Fock type exact exchange (HFX) into semilocal density functionals. The Coulomb potential that appears in the exchange interaction is 0 screened, so it has a shorter range than 1/r. Here, to reduce the big prefactor involved in the computation of the HFX po- Energy (eV) tential matrix elements, we fit the NAO of the basis set with Gaussian-type orbitals, specially suited to computing the four center electron repulsion integrals (ERIs) in a straightforward -5 and efficient analytical way. An example of this fitting for the 2s and 2p atomic orbitals basis set of the oxygen is shown in Fig. 4. The LIBINT package46 is required to calculate (a) primitive ERIs, where recursive schemes of the Obara-Saika47 -10 Γ L K T Γ X method and the Head-Gordon and Pople’s variation48 thereof 5 are implemented. ERIs are calculated in the first SCF cycle and then stored in disk. Only the ERIs with non-negligible contributions are calculated, keeping the HFX Hamiltonian also sparse. 0 This HSE06 functional has been used to compute the band Energy (eV) structure of bulk Si [diamond structure; Fig. 5(a)] and BaTiO3 [cubic structure; Fig. 5(b)] with a double-zeta polarized basis set at the equilibrium lattice constant of the Perdew-Burke- -5 Ernzerhof functional49 within the Generalized Gradient Ap- proximation (5.499 Å for Si and 4.033 Å for BaTiO3 ). In both cases, the gap is opened with respect to the value ob- (b) tained with the semilocal functional. In bulk Si the band gap -10 is indirect: the top of the valence band is located at Γ and the Γ L K T Γ X bottom of the conduction band at a point along the Γ → X high-symmetry line. It increases from 0.64 eV within GGA FIG. 3. Band structure of NiO in the undistorted rock-salt type struc- to 1.00 eV with the hybrid functional, in good agreement with ture with rhombohedral symmetry introduced by a type-II antiferro- the experimental value of 1.17 eV50 . For the case of the per- magnetic order. The experimental lattice spacing is used. The bands obtained within GGA-Perdew-Burke-Ernzerhof functional (panel a), ovskite oxide BaTiO3 , the band gap is also indirect, from R to and with a Hubbard U correction of 4.6 eV applied on the d-orbitals Γ, and its value increases from 1.87 eV with GGA to 3.28 eV of Ni (panel b), as in Ref. 31, are shown. The zero of the energy is with the HSE06 functional, almost matching the experimental set at the top of the valence band. value of 3.2 eV estimated by Wemple in the cubic phase51 . F. Spin-orbit coupling The capability to include the spin–orbit (SO) interaction in D. Van der Waals functionals S IESTA and in the analysis tools is seen as a strategic asset for the project in view of the recent interest in topological in- An efficient calculation of van der Waals (vdW) sulators and quasi–two–dimensional systems with important functionals35,36 was developed and first implemented in spin–orbit effects, like some of the transition metal dichalco- S IESTA using a polynomial expansion in the local variables genides. Also, it brings the possibility to obtain the magnetic (q1 , q2 ) of the nonlocal interaction kernel Φ(q1 , q2 , r12 ) and a crystalline anisotropy (MCA) (change in the total energy of Fourier expansion in the relative position r12 37 . As a result, the system upon changing the spin quantization axis). the scaling of the vdW computation decreases from O(N 2 ) to In a standard collinear-spin DFT calculation, the total KS O(N log N) and it becomes marginal within the overall cost. Hamiltonian is represented by two independent spin–blocks, This scheme was later extended16 to a more complex kernel38 σ σ [σ =↑, ↓]. However, when the SO coupling is included, Ĥµν of the form Φ(n1 , |∇n1 |, n2 , |∇n2 |, r12 ), and it has been applied off–diagonal spin blocks arise (i.e., there are non–zero cou- to a large variety of systems, like carbon nanotubes37 , hydro- plings between the two spin components). Therefore, and gen adsorption39,40 , or liquid water41 . similar to the non-collinear spin case, the Hamiltonian be- S IESTA 7 10 1.4 1.2 1.0 5 0.8 Rnl (r) rl Energy (eV) 0.6 0 0.4 0.2 -5 0.0 0.5 1.0 1.5 2.0 2.5 r (Bohr) -10 5.0 (a) 4.0 -15 L Γ X W K Γ 8 Rnl (r) 3.0 rl 2.0 6 1.0 4 0.0 0.5 1.0 1.5 2.0 2.5 r (Bohr) Energy (eV) 2 FIG. 4. Gaussian fits of the radial part of oxygen 2s (a) and 2p (b) or- bitals using 6 Gaussian functions. The orbitals to fit are represented 0 by blue dots and the corresponding Gaussian expansions by green continuous lines. Dashed vertical lines represent the standard devi- ations of individual Gaussians and a red continuous line marks their -2 upper limit. The orbitals are set to zero in the yellow area, marking their cutoff radii. -4 (b) comes a full 2×2 matrix in spin space -6 Γ X M R Γ M X R ↑↑ ↑↓ ! KS Ĥµν Ĥµν Ĥµν = ↓↑ ↓↓ (2) FIG. 5. Band structure of (a) bulk Si in the diamond structure, and Ĥµν Ĥµν (b) bulk BaTiO3 in the cubic structure obtained with the Perdew- Burke-Ernzerhof functional (red lines) and with the HSE06 hybrid where µν subindexes refer to the S IESTA basis orbitals. The functional (black lines). The zero of energies have been set to the fully relativistic Hamiltonian Ĥ KS is expressed as a sum of the valence band maximum. kinetic energy T̂ , the scalar-relativistic pseudo-potential part in the form of Kleinman–Bylander projectors V̂ KB , the spin– orbit V̂ SO term and the Hartree V̂ H and exchange–correlation In the current implementation the SO term is included V̂ XC potentials: non-perturbatively, so that the fully relativistic Hamiltonian is solved self-consistently after extending the Kohn–Sham Ĥ KS = T̂ + V̂ KB + V̂ SO + V̂ H + V̂ XC (3) wave–functions to full spinors. Two different approaches have been implemented in S IESTA to account for the SO term, V̂ SO : The first three terms of the right hand side do not depend on - on–site approximation: the charge density, ρ(r), and therefore do not change in the Based on the work of Fernández-Seivane et al. 53,54 , self–consistent cycle, while V̂ SO and V̂ XC are the only spin– only the intra–atomic SO contributions within each l– dependent terms that couple both spin components. shell of each atom are considered. In this approach the In order to compute the MCAs, different orientations of SO terms are obtained from analytical simple expres- the spin quantization axis need to be considered. This may sions for the angular integrals while the radial integrals be done by rotating either V̂ SO (as done by Cuadrado and are computed numerically. Cerdá 52 ) or the density matrix, which is the approach cur- rently followed by S IESTA for compatibility with the non– - off–site approach: collinear case. Here, V̂ SO is built following the Hemstreet S IESTA 8 formalism52,55 whereby a fully-relativistic pseudo- • The Fermi Operator Expansion method (FOE)63 uses potential (FR-PP) operator is constructed in a fully the formal relationship between Hamiltonian and separable form, i.e., non–local in the radial part density-matrix, ρ̂ = fFD (Ĥ − µ), where fFD is the as well as in the angular variables, in order to Fermi-Dirac function. A simple polynomial expansion substantially reduce the computational cost. The of fFD can then be used to obtain ρ̂ without diagonal- necessary l j Kleinman–Bylander projectors may be ization. This method is implemented in the CheSS li- either constructed by S IESTA itself from relativistic brary64 , developed within the BigDFT project65 . semilocal PPs, or directly read from appropriately generated PSML files, as provided by the Pseudo– • The PEXSI method66,67 uses a pole expansion of fFD to Dojo project25,26 . Moreover, we note that the FR-PP get ρ̂ in the form: formalism (as well as the original one implemented in Ref. 52) uses the correct normalization constants ! P ρ Cl±1/2 , in contrast with what was erroneously stated in ωl ρ̂ = Im ∑ (4) Ref. 56. l=1 H − (zl + µ)S Although we consider the off–site approach more accurate, ρ as it includes inter-shell and inter-atomic SO couplings, both where ωl and zl are the weights and poles for the cor- approximations yield very similar results in most of the tested responding expansion of the Fermi-Dirac function. The systems, with relevant qualitative differences only found in a number of poles needed is significantly smaller than for few specific cases. Furthermore, the construction of the Vµν SO the polynomial version of the FOE, as its dependence matrix is very fast under both schemes and involves a tiny on the spectrum size is only logarithmic. fraction of the entire self-consistent calculation. It would appear that having to invert matrices would still render this approach cubic-scaling, but in fact only selected elements of ρ̂ have to be actually com- G. New electronic-structure solvers puted. This “pole expansion and selected inversion” method offers a reduced complexity (at most O(N 2 ) For most problems, S IESTA spends the largest fraction of for dense systems, and O(N) for quasi-one-dimensional cpu-time in the solver stage (solution of the generalized eigen- systems), and trivial parallelization over poles, so it is value problem HΦ = εSΦ). The stage devoted to the calcu- well-suited for very large problems on large machines. lation of the hamiltonian H and overlap S is typically much For example,68 computed the electronic structure of lighter weight, as those matrices are intrinsically sparse due large (up to 11,700 atoms) graphene nanoflakes using to the use of a finite-support basis set. Accordingly, S IESTA’s S IESTA-PEXSI. performance is almost completely linked to the use of appro- priate external solver libraries. • The electronic structure problem can also be cast as a Over the past few years we have expanded the choices avail- minimization problem (of an extended functional) with- able to users and refined the relevant interfaces. Initially, we out orthogonalization. When additional localization added support for new individual solvers as detailed below, constraints are put in place, the original linear-scaling but recently we have consolidated some of the most impor- method in S IESTA results. Without the extra localiza- tant functionality under a new common interface to the ELSI tion constraints, the cubic-scaling Orbital Minimization library of solvers57,58 . Method (OMM)69 can be competitive with respect to diagonalization, as data can be reused across scf-cycle steps. 1. Solvers with a native interface Diagonalization (solution of the generalized eigenproblem 2. The ELSI interface appropriate for non-orthogonal orbitals) is the default method for obtaining the density-matrix in S IESTA. A number of stan- We have considerably extended the range of solver choices dard routines are contained in the S CA LAPACK library59 , and the performance enhancement possibilities of the code but more efficient alternatives are possible. In particular, the with the integration of the open-source ELSI library (https: ELPA library60–62 uses an extra intermediate step in the tridi- //elsi-interchange.org), that provides a unified soft- agonal conversion of the matrices to obtain better scalability ware interface that connects electronic structure codes to var- and significant speedups over S CA LAPACK. An interface to ious high-performance solver libraries to solve or circumvent ELPA is offered in S IESTA, so this solver can be used as a eigenproblems encountered in electronic structure theory57 . drop-in replacement for S CA LAPACK throughout the code. ELSI also ships with its own tested versions of the individual In addition, S IESTA has implemented interfaces to several solver libraries, but additionally, linking against already com- methods not based on diagonalization. In most cases, the piled upstream versions from each solver library is supported use of a finite-support basis set, leading to the appearance of as much as possible. sparse matrices, is a significant factor to achieve good perfor- The ELPA, OMM, and PEXSI solvers, which had their mance: own ad-hoc interfaces as described in the previous section, S IESTA 9 are now available through ELSI, which also supports other over trial points for the determination of the chemical- conventional dense eigensolvers (EigenExa70 , MAGMA71 ), potential. sparse iterative eigensolvers (SLEPc72 ), and linear scaling density matrix purification methods (NTPoly73 ). As sketched • Mixed-precision support: The ELPA solver can be in- in Fig. 6, an electronic structure code interfacing to ELSI au- voked in single-precision mode, which can speed up the tomatically has access to all the eigensolvers and density ma- initial steps of the electronic self-consistent-field (scf) trix solvers supported in ELSI. In addition, the ELSI inter- cycle. face is able to convert arbitrarily distributed dense and sparse • Accelerator offloading: The ELPA library offers GPU matrices to the specification expected by the solvers, taking support in some kernels62 , and there is scope for ex- this burden away from the electronic structure code. A com- tending it to more kernels. ELSI also offers an interface prehensive review of the capabilities in the latest version of to the accelerator-enabled MAGMA library. Finally, the ELSI, including parallel solution of problems found in spin- PEXSI developers are working on adding GPU support polarized systems (two spin channels) and periodic systems to the solver. (multiple k-points), scalable matrix I/O, density matrix ex- trapolation, iterative eigensolvers in a reverse communication 1024 interface (RCI) framework, has recently been completed58 . Scalapack Siesta-ELPA ELSI-ELPA multi-k-point 512 Ideal CPU time (s) 256 128 64 32 64 128 256 512 1024 MPI ranks FIG. 7. Performance improvement from the use of the extra level of parallelization over k-points in S IESTA using the ELSI interface with the ELPA solver, compared to the previous diagonalization scheme (using both the standard S CA LAPACK solver and the existing ELPA interface in S IESTA). The system is bulk Si with H impurities, with 1040 atoms, 13328 orbitals, and a sampling of 8 k-points. The multi- k scheme is able to stay closer to ideal scalability for larger numbers of MPI processes. FIG. 6. Interaction of the ELSI interface with electronic structure codes. ELSI serves as a bridge between electronic structure codes and solver libraries. An electronic structure code has access to vari- ous eigensolvers and density matrix solvers via the ELSI API. When- H. Time dependent DFT ever necessary, ELSI handles the conversion between different units, conventions, matrix formats, and programming languages. Time-dependent density-functional theory (TD-DFT) was first implemented into S IESTA in its real-time propagating With the common interface in place, any additions and en- form. It was first described in Ref. 74, and then briefly in hancements to the supported solvers can be used in S IESTA Ref. 2. It was based on the Crank-Nicolson algorithm, by with almost no code changes. This is particularly relevant for which, the effect of the evolution operator for an infinitesimal performance enhancements. For example: time step • Further levels of parallelization: A feature common in   Û(t0 + ∆t,t0 ) = exp −iĤ(t)∆t (5) principle to all solvers is that the S IESTA-ELSI inter- face can exploit the full parallelization over k-points on the wave-function coefficients matrix at a given time t0 , and spins mentioned above. This means that these cal- c(t0 ) is approximated by culations can use two extra levels of parallelization in ∆t −1     the solver step beyond the standard one of paralleliza- ∆t c(t0 + ∆t) = S + iH(t0 + ∆t) S − iH(t0 ) c(t0 ) tion over orbitals (see Fig. 7). 2 2 (6) • The new version of the PEXSI solver integrated in where ∆t represents the finite time-step resulting from time ELSI can achieve the same level of precision with discretization, and S and H represent the overlap and Hamil- fewer poles, and offers an extra level of parallelization tonian matrices, respectively, in the representation given by S IESTA 10 a non-orthogonal basis set, as used by S IESTA. That expres- steps per atomic motion step, if the nuclei are still significantly sion is obtained from equating the first-order evolution of the slower than electrons, for instance. The implementation of the coefficients forward, from t0 to t0 +∆t/2, to the backward evo- Crank-Nicolson part is the same for both the fixed and moving lution from t0 + ∆t to the same intermediate step. basis. It can be further simplified to The square root and inverse square root are calculated by first computing its eigenvalues and eigenvectors, ∆t −1     ∆t c(t0 + ∆t) = S + iH(t0 ) S − iH(t0 ) c(t0 ) (7) S = U E U †, (12) 2 2 for a smooth-enough variation of the Hamiltonian itself and where E is a diagonal matrix with the eigenvalues of S. And a small enough ∆t, thereby avoiding the self-consistency im- U is a square matrix with the eigenvectors of S as its columns. plied in propagation using Eq. 6. In Section III H 3 below, Then, recent developments on efficient treatments of Eq. 6 beyond Eq. 7 are presented. Here we describe the parallelization and S1/2 = U E 1/2U † , and S−1/2 = U E −1/2U † related features in the TD-DFT implementation now found in where E 1/2 and E −1/2 are obtained by replacing diagonal el- standard S IESTA releases. ements of E with their square root and inverse square root The implemented propagation is based on Eq. (7) (with the (in the latter case neglecting those eigenvalues below certain improvement possibilities described in Section III H 3), but threshold value), respectively. proper consideration must be taken of the fact that not only The two-stage algorithm has been implemented in S IESTA the coefficients change in time, but also the basis set and the in parallel, allowing for k-point sampling and for collinear Hilbert space spanned by it when the atoms move. An analy- spin. The initial occupied states to be propagated are read sis of the geometrical implications of this fact is presented in from a file. S IESTA is prepared to run a conventional DFT Ref. 75. The time-dependent Kohn-Sham equation calculation of whatever relevant initial state, and write a wave- function continuation file that acts as initialization of the ul- H|ψi = i∂t |ψi (8) terior S IESTA run in real-time TD-DFT mode. As it stands, becomes S IESTA evolves states defined as fully occupied; partial occu- pations are not currently supported. Hc = i S (∂t + D)c (9) where H, S, and c are the Hamiltonian, overlap and coeffi- 1. Parallelization cients matrices, respectively, as before, and the D matrix is the connection in the manifold given by the evolving Hilbert The two-step procedure described above requires matrix- space75 , Dµν = hφµ |∂t φν i, for φµ and φν basis functions. matrix and matrix-scalar multiplication, matrix addition, and A way of taking such evolution into account in the matrix inversion, plus the diagonalization of the overlap ma- discretized implementation was proposed by Tomfohr and trix for the Löwdin step. Since only the occupied states are Sankey76 , and is based on a Löwdin orthonormalization. The propagated, the c matrix is rectangular N × N , that is, num- scheme consists of two steps. First the wavefunctions are ber of propagating states × number of basis functions, while propagated using both S and H at time t0 using Eq. (7), but S and H are square, N × N . The computation of the overlap to an auxiliary coefficient matrix c̃, and Hamiltonian matrices is handled by pre-existing S IESTA routines, which are already parallelized and well-optimized ∆t −1 for HPC environments69,81 .     ∆t c̃(t0 + ∆t) = S + iH(t0 ) S − iH(t0 ) c(t0 ) . (10) The parallelization of the propagation following Eqs. (10) 2 2 and (11) is done simply exploiting the MatrixSwitch library82 , Then, the propagation is followed by a change of basis opera- which allows for an abstracted manipulation of matrices, the tion (only needed if the ionic positions have changed), details of parallelization, data formats, conversions etc. being taken care of underneath. In this case, MatrixSwitch is called 1 1 c(t0 + ∆t) = S− 2 (t0 + ∆t) S 2 (t0 )c̃(t0 + ∆t). (11) to use the BLACS83 and S CA LAPACK84 libraries, meaning that this part of the code is run on dense-matrix infrastructure, This algorithm is unitary by construction, and so the preser- as already done with conventional diagonalization solvers. As vation of orthonormality is guaranteed, regardless of the size for the latter, although the H and S matrices are sparse, the c of ∆t. As discussed in detail in Ref. 75, this algorithm can be matrix is dense. shown not to be entirely consistent with the connection rep- Conversion between storage formats is an important con- resented by the D matrix defined above. Nevertheless, the sideration here. The native matrix storage format employed discrepancies due to the mentioned inconsistency have been by S IESTA is a compressed sparse column (CSC) scheme shown to be small in a series of studies using this formal- with a one-dimensional block-cyclic distribution (1D-BCD) ism77–80 , at least for low atomic velocities. The practical ben- over MPI processes. A block-cyclic distribution is needed by efit of separating the two procedures is to perform the change BLACS and S CA LAPACK package. The matrices can there- of basis only when necessary, allowing for many electronic fore be temporarily converted from sparse to dense using the S IESTA 11 same parallel distribution; this is a very efficient operation, 2. TD-DFT Beyond the released version since no MPI communication is necessary. It should be noted that a two-dimensional (2D) BCD is known to be more effi- There are many possible (and feasible) improvements on cient in terms of parallel scaling69 . The conversion from 1D to what has been described above, some of them in the pipeline. 2D does however carry a heavier cost, as MPI communication From a fundamental point of view, the Löwdin step will be is inevitable. replaced by another basis-changing step, in the direction of The parallel efficiency of our implementation is therefore what was proposed in Ref. 75. It is needed if atoms move at chiefly determined by that of the underlying S CA LAPACK velocities of around 1 atomic unit or more (1 a.u ∼ c/137, drivers. The matrix inversion in Eq. (10) is performed using being c the speed of light). In that case the diagonalization LU factorization. For the diagonalization of the overlap ma- step may be avoided (or replaced by the N ×N diagonalization trix we have implemented the option of using either a standard of the overlap matrix for the evolving states, instead of the diagonalization approach (tridiagonal reduction followed by N × N for the basis set overlap). the implicit QR algorithm) or a divide-and-conquer algorithm For fast moving atoms within Ehrenfest dynamics, there is as described in Ref.85. The latter is known to scale better with also the need to implement correction terms to the forces re- system size. lated to both the change of basis and the rotation of the time- The scaling with number of processors is very similar to the dependent Hilbert space (the intrinsic curvature of the man- scaling of a conventional DFT S IESTA run using diagonaliza- ifold). These terms are well known,86 and their geometrical tion as the solver option, since both procedures are run on meaning in terms of the relevant curvature75 will appear in routines of analogous scaling within the same dense-matrix- Ref. 87. They are being tested and should be incorporated algebra library. Figure 8 shows the relative share in the to- shortly to a visible branch in the open source repository, to tal running time of the three main procedures involved: the be later merged into the trunk, and further incorporated into Crank-Nicolson algorithm, the change of basis, and the cal- S IESTA releases. culation of the SCF Hamiltonian plus other minor processes For efficiency, iterative inversion options will be explored in S IESTAsuch as building the density matrix. This was per- replacing the present LU implementation in S CA LAPACK, formed for a system of 5000 Ge + 1 He atoms described with a and quite a few possibilities exist to incorporate more so- single-zeta polarized basis set. The Crank-Nicolson algorithm phisticated algorithms to the described operations. What has takes about 18% of the total time on 30 processors, which in- been described is robust and quite transparent, but the Ma- creases to 25% on 316 processors. Instead, the change of basis trixSwitch abstraction should allow easy implementation of procedure takes about 38% of the total time on 30 processors, other techniques. which decreases as the parallelization increases, reflecting its better scaling properties. The Löwdin step is the most ex- 3. Improved real-time propagators pensive operation on all numbers of processors, which affects TD-DFT simulations (and only those steps) involving atomic motion. Eq. (7) represents a fast approach of electronic propaga- tion in real-time TD-DFT, especially suited to study systems where the perturbation of the electronic density is relatively Crank-Nicolson small (e.g. optical linear response74 ). If one is interested Change-basis Hscf + the rest in simulating systems with heavily perturbed electronic den- 1 sities by external forces (like those exerted by intense laser fields or fast atom collisions for instance), one should choose Relative share (time) 0.8 a more elaborate propagation scheme that preserves better the time-reversibility of the propagator operator. Some of the au- 0.6 thors introduced in Ref. 79 an extrapolation algorithm to study the stopping power of prototype semiconductors. Briefly, the 0.4 method uses Eq. (7) with an extrapolated Hamiltonian: ∆t −1     ∆t 0.2 c(t0 + ∆t) = S + iHext S − iHext c(t0 ), (13) 2 2 0 where the extrapolated Hamiltonian Hext reads 50 100 150 200 250 300 1 Processes Hext = H(t0 ) + ∆H (14) 2 and FIG. 8. The relative share of the total running time for the Crank- Nicolson algorithm, the Löwdin step, and the rest of the program op- ∆H = H(t0 ) − H(t0 − ∆t). (15) erations (including the building of the SCF Hamiltonian) for a system of 5000 Ge atoms and one He projectile, using 30-316 processors. Additionally, the user is given the option to divide each prop- agation step ∆t into n sub-steps in an effort to increase the ac- curacy of the first-order expansion underlying the derivation S IESTA 12 of Eq. (13). In this case, the final equation reads code89 in the case of a double ionization of a uracil molecule −1  in the gas phase90,91 . Simulations of this type provide ac- n   j ∆t j ∆t cess to the ultrafast electronic dynamics that occurs at the atto c(t0 + ∆t) = ∏ S + iHext S − iHext c(t0 ), (16) j=1 2n 2n and femto time-scales in the ionized genetic material (DNA and RNA) as a consequence of collisions with proton or car- with bon beams92 . This particular simulation addresses the frag- mentation pattern of a doubly-ionized uracil molecule (its j 1 1 deepest Kohn-Sham orbital is empty) using the BLYP den- Hext = H(t0 ) + ( j − )∆H (17) n 2 sity functional93,94 . The technical details for the CPMD simu- Recently we introduced a third algorithm for propagation. lation used as a reference can be found in Refs. 90 and 91. Leaving aside in this description the complications associated S IESTA calculations using the same density functional and with the possible subdivision of each time-step, the new al- the integrators described above use a DZP basis set. As can gorithm is based on a two-step scheme where the electronic be seen in Fig. 9 the standard S IESTA implementation cannot wavefunction is first propagated until half of the step, ∆t/2, properly deal with such a highly excited system. The extrap- using extrapolation as in Eq. (13), olation scheme in Eq. (13) already provides a large improve- ment and gives an energy conservation similar to the CPMD ∆t −1 simulations in Refs. 90 and 91. Finally, the two-step algo-     ∆t ∆t c(t0 + ) = S + iHext S − iHext c(t0 ), (18) rithm further improves the energy conservation. For smaller 2 4 4 time steps the improvements given by the two-step scheme are then an explicit calculation of the half-step Hamiltonian, even more clear, as shown in Fig. 10. H(t0 + ∆t/2), is performed using the coefficients c(t0 + ∆t/2) obtained from Eq. (18). In a second step, the coefficients are evolved from the beginning of the step, c(t0 ), to the full step, c(t0 + ∆t), using the half-step Hamiltonian: ∆t ∆t −1     ∆t ∆t c(t0 +∆t) = S + iH(t0 + ) S − iH(t0 + ) c(t0 ) 2 2 2 2 (19) This approach, although increasing the CPU time by around ∼35% as compared to the two previous schemes, allows for better energy conservation for highly perturbed systems where the Kohn-Sham potential heavily varies in time. FIG. 10. Energy drift comparison as in Fig. 9, for time steps ∆t = 0.24 and 1.2 attoseconds. 4. Electronic stopping of atomic projectiles Let us finish the TD-DFT section with a brief mention of its successful application to the problem of simulating the ex- citation of the electrons of a condensed matter system when traversed by a high-energy atomic projectile (so-called elec- FIG. 9. Energy drift in an energy conserving TDDFT simulation tronic stopping, since the electrons slow down the projectile). of ionized-core uracil for the three propagation methods described This physical problem is very relevant to questions of interest here, as compared with the CPMD implementation, for a time step to the nuclear and aerospace industries, as well as to the treat- of ∆t = 0.24 attoseconds (0.01 atomic units). ment of cancer. In spite of its great relevance and of its having been researched since Rutherford’s experiment in 1911, the In order to provide a more quantitative comparison between understanding of electronic stopping processes has been es- the three schemes, namely, the default propagation of Eq. (7), sentially limited to either weak effects in the linear-response the extrapolation propagation of Eq. (13) and the two-step regime or beyond linear but only for target systems close to propagation of Eqs. (18) and (19), we compare their perfor- the homogeneous electron liquid (jellium). mance versus the P-TDDFT implementation88 of the CPMD An earlier version of the TD-DFT implementation in S IESTA 13 S IESTA74 , allowed the first explicit first-principles simulation c∗jα ∆iαβ ciβ and ∆iαβ = ∂ Hαβ − εi ∂ Sαβ . The   where Ai j = ∑αβ εi −ε of electronic stopping, for protons and antiprotons in LiF, a j change in the density matrix is then given by wide-band-gap insulator, which was quite successful.95 The difference of sign between protons and antiprotons produces a significant difference in the stopping power (rate of energy all  ∂ ρµν = ∑ ni ∂ c∗iµ ciν + ni c∗iµ ∂ ciν + ∂ ni c∗iµ ciν =  excitation) beyond the linear-response paradigm (the Barkas effect), and the insulating character of the target makes it in- i accessible to the jellium paradigm. The success stimulated all all   all further studies along this line77–80 using improved versions of = ∑ ni ∑ A∗i j c∗jµ ciν + c∗iµ Ai j c jν + ∑ ∂ ni c∗iµ ciν = i j i TD-DFT in S IESTA, as described here. Fig. 11 displays the electron deformation density around a proton displacing in a all unocc  A∗i j c∗jµ ciν + c∗iµ Ai j c jν  bulk Ge target.79 They were also followed by analogous sim- = ∑ ni ∑ i j ulations using plane-wave codes by an increasing number of occ all groups (for a review see Ref. 96), although the latter calcula- tions do demand considerably larger computational resources. − ∑ n j c∗iµ c jν ∑ c∗jα ∂ Sαβ ciβ + ∑ ∂ ni c∗iµ ciν (21) ij αβ i and a similar expression applies to the change in the energy- density matrix. The change in the occupation of the electronic state can be computed from the change in its eigenenergy εi = ∑αβ ciα ∆iαβ ciβ and it is relevant in metals, for states close to the Fermi level. The Fermi level can also be shifted by the perturbation, and it can be determined through conservation of the number of electrons in the system, Ne . FIG. 11. Electron deformation density isosurfaces (blue positive, red negative) for a proton displacing leftwards, at a velocity of 1 a.u. in the bulk of a Ge crystal. I. Density Functional Perturbation Theory The original implementation of Density Functional Pertur- bation Theory, as a post-processing and independent code (LINRES97 ) has been recently merged into S IESTA. It allows to compute the phonon dispersions using a supercell approach FIG. 12. (Left panel): Comparison of the performance of the new (Γ-point phonons). Both LDA and GGA functionals can be DFPT approach with the conventional finite-differences method in used (through LibXC). Calculation of the perturbed Hamilto- S IESTA. The time required to compute a whole row of the dynam- ical matrix (derivative of the forces on all atoms when one atom is nian and overlap matrix elements follows the same method- displaced in x, y, z directions) is plotted as a function of the number ology as for ground-state calculations, with similar computa- of orbitals in carbon fullerenes of different sizes. (Right panel) Per- tional costs, which are comparable to those obtained with a formance of the alternative algorithm described in the text (based on finite difference approach. Figure 12 shows a comparison be- Eq. (22), blue circles), as compared to the original implementation tween both methods for model fullerene-type systems of dif- based on Eq. (20) (black diamonds). The initialization (computing Ξ ferent sizes. and Ω) is the most time-consuming step, although it has to be per- The solution of the Sternheimer equation, and calculation formed only once, and can be used for all perturbations (each atomic of the perturbed density matrix is the most demanding step. It displacement). The new algorithm becomes more efficient for sys- requires the perturbed coefficients of the electronic wavefunc- tem sizes larger than the threshold value (green dashed line). tions to be obtained, . Obtaining ∂ ρµν is the most computationally expensive part ∂ Hαβ − εi ∂ Sαβ   ∂ ciµ = ∑ ∑ c∗jα ciβ c jµ = ∑ Ai j c jµ (20) of the code. While the computation of ∂ Hµν basically has j αβ εi − ε j j linear scaling with the system size, the matrix Ai j scales as S IESTA 14 Nb2 · M, where Nb is the number of basis functions, and M is for Ne ≥ 1 electrodes. the maximum number of neighbour orbitals for any orbital in the system. Equation (20) then requires Nb3 operations for ρ = ρeeq + ∑ ∆ee0 ≡ ρeneq , (23) 0 e 6=e each atomic perturbation, and the change in the density matrix requires Nb2 · M loops. An alternative approach that offers a i ZZ h i better computational scaling for systems with a gap has also ρeeq ≡ dkdε Gk (ε) − G†k (ε) nF,e (ε)e−ik·R , (24) 2π BZ c∗jα c jβ been tested. If we define Ξiαβ = ∑ j 1 ZZ εi −ε j , we obtain: e ∆e0 ≡ dkdε Ae0 ,k (ε)e−ik·R nF,e0 (ε) − nF,e (ε) ,   2π BZ   (25) ∂ ciµ = ∑ ∑ α µ αβ ciβ = ∑ Λiµβ ciβ Ξ i ∆i (22) where ρeeq is the equilibrium density matrix for electrode e, G β α β is the retarded Green’s function matrix, and ∆ee0 is the cor- where Λiµν = ∑η Ξiη µ ∆iην is a smooth function of εi and can rection to the equilibrium part. The spectral function Ae = be described by a Chebyshev’s expansion with a few selected GΓe G† and carries electrons from the electrode e. Finally, energy points and their corresponding weights: nF,e is the Fermi function with chemical potential denoted by the electrode e. It is evident that the Fermi functions depend (l) (l) (l) on the chemical potential and the electronic temperature in ∂ ciµ = ∑ ω̃l,i ∑ Λµβ ciβ = ∑ ω̃l,i ∑ Ξα µ ∆αβ ciβ l β l αβ the associated electrodes. By using different temperatures for each electrode one can calculate thermoelectric effects due (l) Notice that Ξαβ is perturbation-independent, and could be to different reservoirs having separate electronic temperatures self-consistently. computed only once and used for all the possible atomic dis- We note that T RAN S IESTA uses a multiple complex energy- placements, with a cost that scales as Nb2 · Nl , with Nl being contour algorithm to more accurately describe the total den- the number of Chebyshev’ polynomials. The change in the sity ρ. It does this by weighing each ρeneq using a simple electronic density is then given by scheme98 (Sec. 3.2). So far, few multi-electrode calculations +∂ ρµν ∼ ∑ c∗iµ ∂ ciν = ∑ Λνη ∑ c∗iµ ω̃l,i ciη = (l) have been performed so the importance of the multiple con- i lη i tour algorithm is currently unknown100–102 . However, for the (l) (l) (l) (l) (l) well-known 2-electrode problem it allows smoother conver- ∑ Λνη Ωη µ = ∑ Ξγν ∆γη Ωη µ gence properties15 . lη lηγ In the latest T RAN S IESTA we implement 3 different inver- sion algorithms; i) a block-tri-diagonal algorithm (BTD), ii) where only the central term requires self-consistency, and (i) MUMPS sparse algorithm and iii) a dense algorithm (LA- Ωαβ = ∑ j c∗jα ω̃i, j c jβ . Although computing the change in the PACK). The performance of these (speedup compared to density scales as O(Nb2 M), most of the computational cost T RAN S IESTA 3.1) is summarized in Fig. 13. Since the BTD is required in an initialization step to obtain Ξ and Ω that algorithm is linear scaling for constant width it can easily out- are perturbation-independent, enabling the extraction of the perform the dense algorithm by a factor 100. This perfor- whole dynamical matrix with O(Nb3 ) operations. A prelimi- mance gain is also important for the memory footprint which nary serial calculation for Cn fullerenes shows that the thresh- enables even larger systems. The BTD algorithm favors long old system size for the new algorithm to become more ef- and narrow systems, but uses less memory for all types of sys- ficient than the original implementation lies at around 650 tems. atoms. This value can be conveniently reduced by an efficient A recent addition to the T RAN S IESTA package is the use parallelization of the initialization step. of real space self-energy terms99,103 . These self-energies are semi-infinite in more than 1 direction and can thus be used as surrounding electrodes, for e.g. single defects in 2D or 3D J. TranSiesta structures or line defects. Real space self-energies are supe- rior to BZ integrated quantities since they correctly describe The transport code T RAN S IESTA, initially developed by the infinite bulk by leaving out image couplings and also re- Brandbyge and co-workers15 , enables open-boundary con- moves the need for k-point sampling. When taking into ac- dition calculations by extending periodic regions with bulk count the complete procedure for a T RAN S IESTA calculation electrodes. It is based on the non-equilibrium Green func- the real space self-energies provide an increased throughput tion formalism which allows biased calculations. T RAN - since the SCF k-point sampling and the subsequent k-point S IESTA has been completely re-written and now uses ad- sampled transport calculation are completely removed99 . vanced inversion algorithms, enables Ne ≥ 1 electrodes, al- Additionally, TB TRANS enables calculations of user de- lows thermo-electric calculations, performing real-space cal- fined tight-binding models and also interfaces to phonon culations (without k-points) and adds phonon transport calcu- transport using the Hessian matrix (program named PH- lations using the Hessian98,99 , TRANS). The phonon Green function is similar to the electron The non-equilibrium Green function formalism can be −1 Gq (ω) = (ω 2 + iη 2 )I − Dq − Σq (ω) ,  summarized in the following equations which are generalized (26) S IESTA 15 LAPACK MUMPS minimized.108 In a practical procedure to construct Wannier a) BTD b) c) functions, a set of localized functions is used to generate an 60 100 10 90 9 initial guess for the unitaryLinear transformations. These localized 50 0V 0.75 V Opteron2354 80 8functions should be roughly located on sites where Wannier 70 Opteron6136 40 7functions are expected to be centered and have appropriate Speedup Speedup Speedup 60 6 Xeon2665 30 50 5angular character. In our implementation, Xeon2660v3 we can directly use 40 4the localized atomic orbitals of the basis, or the hydrogenoid 20 30 3 10 20 2 localized functions (including hybrid orbitals), as suggested 10 1in WANNIER 90. 0 0 1,500 3,000 4,500 6,000 1,500 3,000 4,500 6,000 1The 2 3 Wannier 4 5 6 7 functions 8 9 10 provide an exact tight-binding rep- # of orbitals # of orbitals Threads resentation of the dispersion of the Bloch bands. This prop- erty will be exploited to extract in an automatic and user blind FIG. 13. Performance characterization of T RAN S IESTA using a pris- way the parameters required to run multiscale simulations as tine graphene cell (24 atoms wide). Speedup for (a) EGF and (b) described in Sec. III L. NEGF calculations of pristine graphene compared against the dense implementation. The BTD method exhibits more than 40 times the speed of the LAPACK implementation for the largest size. MUMPS gains speed after 5, 000 orbitals. where D is the Hessian and ω the phonon frequency. Finally, inelastic transport involving phonon-excitation can be treated with perturbation theory in a postprocessing step with the I N - ELASTICA package104,105 . K. Wannierization The interface between S IESTA and WANNIER 90106,107 (ver- sion 3.0.0) has been implemented, so the latter code can be called directly from S IESTA as a library, or used as a post- processing tool. WANNIER 90 is an open-source code for gen- erating maximally-localized Wannier functions108,109 and us- ing them to compute advanced electronic properties of mate- FIG. 14. Maximally localized Wannier functions (MLWFs) for rials with high efficiency and accuracy. graphene. Top panel displays the character of σ -bonded combi- The Wannier functions can be considered as a unitary trans- nations of sp2 hybrids. Bottom panel displays the π character of formation (more precisely, a Fourier transformation) of a set the bands with weight on the pz orbitals. Isosurfaces of different of Bloch functions associated with a given manifold of bands. colors correspond to two opposite values for the amplitudes of the We can view the Bloch and Wannier functions as providing real-valued MLWFs. Yellow spheres represent the position of the C atoms, while smaller blue spheres mark the center of the bonding. two different basis sets describing the same manifold of states associated with the electron band manifold in question. The Wannier functions display a number of very interesting prop- Currently, WANNIER 90 can be used as a post-processing erties.110 Among them, we can enumerate: (i) they are local- tool, or it can be directly called from S IESTA in a library ized in real space, each of them concentrated around a given mode. Within this last approach, the unitary matrices that unit cell (see Fig. 14); (ii) Wannier functions centered on dif- transform the Bloch states into Wannier functions are directly ferent cells are translational images of one another; (iii) they accessible in S IESTA, allowing a clear and straightforward in- form an orthonormal basis set; (iv) they span the same sub- terconnection between the two alternatives to span the Hilbert space of the Hilbert space as is spanned by the Bloch func- space. Besides, the use of Wannier functions opens the door to tions from which they are constructed. Because of the gauge a wide range of potential applications. Already implemented freedom in the definition of the phases of the Bloch functions, in S IESTA is the possibility of performing SCF convergence the Wannier functions are not unique. However, the location under the constraint of a rigid shift on the energy associated of their centers in the home unit cell is unique to within a lat- with a given Wannier function to be used to calculate electron- tice vector, i.e. they are gauge invariant.110 The high degree of electron interactions for multiscale simulations as detailed in arbitrariness in the definition of the phases can be exploited to Sec. III L. The interface with the self-consistent dynamical produce unitary transformation matrices between Bloch and mean field theory DMFT W DFT code111 using MLWF has Wannier functions in such a way that a localization functional been already implemented112 . Also, alternative approaches that measures the sum of the quadratic spreads of the Wan- to compute the exact Hartree-Fock exchange in extended in- nier functions in the home unit cell around their centers is sulating systems with a linear scaling computational cost us- S IESTA 16 ing MLWFs have been proposed, being another interesting re- lattice part has successfully been used in several applications. search line for the future.113 One of the main fields of research has been thermal conduc- tivity in perovskites. In particular it was employed to study the electrophononic coupling in SrTiO3 116 and PbTiO3 117 and L. Multiscale methods the proposal of a thermal switch in PbTiO3 .118 It has also been used to study the competition between various ferro- Density Functional Theory can be used as the basis for pa- electric domain structures in PbTiO3 /SrTiO3 superlattices as rameterized multiscale methods, that can be used to carry out a function of strain.119 As a result it was found that tensile simulations including tens or even hundreds of thousands of strains lead to the appearance of chiral ferroelectric vortices atoms.114 First-principles methods are used to produce de- while ferroelectric skyrmions were predicted and experimen- tailed models that are subsequently used to predict properties tally observed for more compressive strain values.119 The cal- that require large-scale simulations. The models are created culated dielectric properties of these superlattices120 are in for specific materials and their accuracy can be systemati- very good agreement with measured values and show very cally improved to converge towards DFT precision. Given the large electric susceptibility consistent with regions of nega- dependence on first-principles, we refer to these methods as tive, static electric permittivity situated at the core of the vor- second-principles DFT (SPDFT) and are run on an indepen- tices and the PbTiO3 /SrTiO3 interfaces. dent code called S CALE -U P.114 SPDFT is based on a division of the total electronic den- sity, n(~r), into a reference (n0 (~r)) and a deformation (δ n(~r)) M. Scripting and integration in external frameworks contributions, n(~r) = n0 (~r) + δ n(~r), (27) An ongoing trend in many areas of computational science is to move away from rigid and monolithic codes, favoring in- where δ n(~r) is considered as a small perturbation with respect stead a more flexible approach in which the internal function- to n0 (~r) that, in non-magnetic cases, represents the ground ality of a program is somehow exposed to the outside world. state of the system114 . This division is then used114 to expand If done in a proper and well-documented way, this can serve to the DFT energy with δ n finding that the zeroth order term, enhance the interoperability of codes with different function- E (0) , corresponds with the full DFT energy for the reference alities, playing to the relative strengths of each, and/or to im- density. The corrections to this reference energy only depend plement new functionalities by combining the available basic on δ n (and parametrically on n0 ) which, given its smallness, blocks. In S IESTA we have followed two different but comple- can be efficiently calculated leading to a fast and accurate ap- mentary routes to these ends: the development of an internal proximation of the full DFT energy. The expansion is usually scripting framework based on the Lua language, which en- taken to second-order, ables new functionality without code recompilation, and the implementation of a formal interface to the AiiDA platform. E ≈ E (0) + E (1) + E (2) + ..., (28) resulting in a stationary problem that is equivalent to Hartree- Fock with the important distinction that the interactions are 1. Lua interface screened by the exchange-correlation potential. In order to keep δ n small the application of the method is restricted to Lua121 is an easy-to-learn and fast scripting language built problems where atomic bonds are not created or destroyed, for embedding. It is very lightweight (its memory footprint i.e. to processes that display an invariant bond topology. is less than 300 kB), and provides very simple ways to inter- The E (0) term represents the exact DFT energy for the refer- face to the data structures and routines of a host program. A ence density. We represent it for a variety of geometries with Lua script, interpreted by the Lua interpreter embedded in the an accurate force-field115 that allows for fast evaluation. The program, can then control the flow of execution and the data. E (1) and E (2) terms account for the changes in the electronic Different user-level scripts can implement new functionalities, structure that are represented by geometry-dependent Wan- without recompilation of the host code. The strategy we have nier functions. Under this basis E (1) becomes a tight-binding followed in S IESTA is based on handling control to the Lua model while E (2) represents electron-electron interactions. interpreter at specific relevant points in the program flow (e.g. The interconnection between the first (S IESTA) and the sec- at the beginning of a geometry step, at the end of a scf step, ond (S CALE -U P) principles simulations is carried out through etc). Lua scripts implement handlers appropriate to the point a python script, M ODELMAKER. Taking a few cutoff dis- they want to hook into, and can request access to specific data tances M ODELMAKER is able to produce a model’s terms and structures. For example, a script intended to implement a bet- automatically carry out DFT simulations with S IESTA to de- ter scf mixing algorithm would be executed after every scf termine the force field, a Wannier Hamiltonian to represent the step, inspecting the convergence data, and changing mixing bands, electron-lattice terms to account how the bands change parameters or schemes, as appropriate. As another example, with geometry and electron-electron interactions to describe, convergence checks over mesh-cutoffs and k-point sampling for example, magnetism. can be performed automatically. While, so far, few publications with SPDFT methods in- The above mixing scenario exemplifies an important area of clude explicit treatment of electronic degrees of freedom, the usefulness of the approach: the prototyping in Lua, (followed S IESTA 17 eventually by a full implementation), of new ideas and algo- Code (1603) 'siesta-4.0-new' Str (3148) Code (2845) 'plstm-4.0' StructureData (3149) CH rithms. We have implemented a number of custom molecu- code protocol stm_code WorkCalculation (3150) structure lar dynamics modes, geometry relaxation algorithms, and ad- CREATE CREATE CREATE CREATE CREATE CREATE CREATE vanced optimization schemes, in a pure Lua library FLOS122 . KpointsData (3155) ParameterData (3153) ParameterData (3154) ParameterData (3156) ParameterData (3157) CALL Bool (3151) Int (3152) 2x2x2 (+0.0,0.0,0.0) parameters basis settings options clean_workdir max_iterations kpoints The code in the library can be re-used, or taken as starting output_structure WorkCalculation (3158) kpoints CALL CALL point for other implementations by users. These user-level output_structure output_parameters retrieved SiestaCalculation (3162) FINISHED remote_folder stm_array scripts can in turn be shared, opening the way to the devel- StructureData (3166) CH output_structure ParameterData (3165) output_parameters retrieved FolderData (3164) output_array ArrayData (3167) remote_folder RemoteData (3163) opment of new functionality with faster turnaround that the parent_calc_folder STMCalculation (3170) traditional approach that needs a careful integration into the output_parameters FINISHED remote_folder retrieved stm_array program’s code base. ParameterData (3173) RemoteData (3171) FolderData (3172) ArrayData (3174) As a specific showcase of the power of the Lua embedding, we have developed a number of variations of the nudged- elastic band method (NEB)123,124 for transition-state search. Previously proposed implementations in S IESTA involved sig- STM image of benzene molecule nificant, hard to maintain code changes, and did not make into the mainstream version. With Lua, we have been able to implement, non-intrusively, not only the standard algorithm, FIG. 15. Automatically generated graph for the execution of an Ai- but a Double Nudged Elastic Band (DNEB)125 variation, and iDA workflow for simulation of STM images. also another version which treats atomic coordinates and lat- tice variables on an equal footing (the variable-cell NEB, or VC-NEB, method126 ). different codes. The integration of Lua functionality in S IESTA has been made possible by the development of an intermediate layer, FLOOK 127 , (for “fortran-Lua-hook”), which provides wrap- N. Utilities for post-processing and supplementary features pers for access to Fortran data structures and subroutines. S IESTA offers several features beyond the core functional- 2. AiiDA plugins and workflows ity of solving the electronic structure problem and performing optional geometry relaxations and molecular dynamics runs. It is worth noting in particular that the atomic character of the The AiiDA framework128–130 provides support for high- basis set enables the use of a very intuitive suite of analysis throughput computations in materials science, keeping full tools, which take advantage of the fact that most of the con- provenance of the calculations and facilitating data handling cepts relating to chemical bonding use the language of atomic and sharing. The framework is open-source, written in orbitals. Python, and designed to support arbitrary codes via a plugin The (partial) density of states, atomic and orbital popula- interface. A plugin for S IESTA has been implemented and tions, and other useful output can be obtained directly from is distributed as the open-source package aiida-siesta131 . the program. The S IESTA distribution includes also several The plugin provides the basic operations of preparing the in- tools in the Util directory for band-structure and wavefunc- put files for a calculation using AiiDA-specific input objects, tion plotting, bonding analysis, etc. Beyond these, special tool and parsing the results and generating AiiDA output objects. packages that implement a specific feature that extends the The AiiDA data are stored in a graph database that keeps a functionality of the main program, or that provide extra op- permanent record of the inputs and outputs of the calculation, tions for visualization or post-processing in general, are avail- and is fully searchable for, e.g. data analytics purposes. able in alternate distribution points. We describe in what fol- AiiDA also provides robust support for the creation of lows the most relevant developments. workflows that incorporate all the necessary steps in the cal- culation of potentially complex properties, together with the proper heuristics and fail-safe features. The aiida-siesta package provides a base workflow and a few workflows for 1. Updates to core utilities standard materials properties, such as band structures. Fig. 15 shows the execution graph of a workflow designed to gener- A number of improvements, enhancements, and additions ate a synthetic STM image from a given structure. Work is have been made to the core utilities shipped with the S IESTA ongoing to implement more complex ones. distribution. In addition to an interface to the computational capabili- There is now a “fat-bands” feature, by which bands can ties of the S IESTA code via the plugin and workflows, the be decorated with information about the relative weight of aiida-siesta package also provides an implementation of given orbitals in each state. The wave-function-related anal- basic objects representing pseudopotential files, notably one ysis tools have been extended to the non-collinear and spin- for PSML. Families of pseudopotentials can be uploaded to orbit case. This includes the COOP/COHP bonding analy- an AiiDA database and shared via the provided mechanisms sis, band-structures, and a new tool spin-texture calculation. for data export and import, facilitating the interoperability of There have been also improvements to band-structure plotting S IESTA 18 utilities and to the visualization of charge densities, potentials, for interpolating the data from the S IESTA output grid (fixed and other magnitudes represented in a real-space grid. by the unit cell dimensions and the MeshCutoff parameter) A band unfolding utility has been added. Based on the onto an arbitrarily cut (and possibly rotated or resampled) par- Fourier decomposition of the Bloch wave functions, it allows allelepipedic box. XCrysDen provides a number of display to perform a “full unfolding” even for non-periodic systems options, including contour lines over grid planes, or isosur- (e. g. liquids) calculated with a large simulation cell. By re- faces. A special feature available in XCrySDen is plotting the folding the fully unfolded bands, from the reciprocal supercell Fermi surfaces. A special script, eig2bxsf, serves to analyze of a perturbed or defective crystal, into the reciprocal unit cell the list of k-points handled by S IESTA, expanding it onto a of the primitive crystal, one recovers the conventional “un- regular sequence, and writing the respective band energies in folded” bands132 . the necessary format. The tools concerning the lattice dynamics have been de- veloped having in mind the Γ phonons calculated for a large 2. sisl enough supercell, that is a typical case in a simulation of molecular crystals or disordered substitutional alloys. For SISL is a Python toolbox that was initially conceived to han- visualization, vib2xsf and vib2vesta place arrows at the dle and manipulate S IESTA/T RAN S IESTA output103 . It has atoms according to the vibration pattern stored in the eigen- since been extended to support other DFT codes, with the aim vectors file (.vectors), produced by the core Vibra util- of offering equivalent operations for them. ity, and can also be used to make animations (sequences of By reading the LCAO outputs from S IESTA one can post- snapshots) of selected vibration modes. Both vib2xsf and process the Hamiltonian and calculate e.g. Brillouin zone in- vib2vesta tools allow the selection of a part of the system to tegrated DOS, wave functions expanded on grids, eigenval- be exposed. ues, band velocities and many more. SISL can process nearly The phdos tool is designed for analyzing zone-center vi- all the S IESTA output files. In particular, it is also able to bration results. As the system is supposed to be large (e.g., post-process data on the real-space grid. Its command line in- a supercell chosen for a periodic crystal), the (artificially terface allows data format changes, e.g. conversion of S IESTA broadened, for convenience) discrete spectrum may serve as XV files to xyz/xsf files or S IESTA binary grid data (VH, VT, a fair approximation to the total density of modes, and if . . . ) to cube/xsf files. weighted with (squared) components of eigenvectors at dif- As it can process density matrices from S IESTA, one can ferent atoms – provide a decomposition into contributions of also use SISL to prepare an input electronic-structure for new different atoms in the total density of vibration modes. calculations, which may be helpful to reduce initial SCF steps. A more sophisticated option is the projection of different SISL also allows creation of custom tight-binding models eigenvectors according to various criteria. The typical system (both orthogonal and non-orthogonal), and since it extracts under study is a supercell in which e.g. an alloying, or some the DFT Hamiltonian matrix one can manipulate the Hamilto- kind of deformation, breaks the underlying perfect periodic- nian to retain certain band-structure features and thus perform ity. Still, some trends related to the latter can be revealed by large-scale simulations133 . This allows calculating far-field appropriate projections. The two obvious cases are the pro- currents using reduced basis-sets with very little loss of accu- jections onto (1) q-vectors of the underlying lattice and (2) racy. irreducible representations of the space group of the underly- The Atomic Simulation Environment (ASE)134 and SISL ing lattice; the corresponding formulas and some results can have a certain degree of overlap in terms of geometry han- be found in Ref. 137. The first type of projection, if done dling functionality. One can easily convert to and from ASE for a sequence of q values, helps to reveal “phonon disper- objects in SISL, thus allowing seamless interaction. sions”, obviously blurred by the broken periodicity, also mak- ing distinction between transversal and longitudinal modes – see Ref. 138 for an example of use. To make the trends more 3. Other post-processing and visualization utilities pronounced, the supercell needs to be sufficiently long in the direction concerned – see, e.g., Fig. 16. The simplest case, The body of utilities contributed by non-core developers a projection onto a single q=0 value, may also be of interest, and other S IESTA users has continued to expand. In partic- since it enhances the modes which are expected to dominate ular, we feature in this section two suites of utilities, one deal- the infrared or Raman spectra, and thus facilitates their com- ing with alternate visualization tools for some S IESTA results, parison with experiment. and another one specifically dealing with lattice dynamics. The symmetry projection may help to isolate in a possibly For structures, the xv2xsf and xv2vesta converters pro- complex spectrum those modes which are expected to domi- cess data from the S IESTA .XV file into the native formats nate according to a given selection rule, again in view of their of XCrySDen135 and VESTA136 , respectively. Each of these verification against the experiments. The group-symmetry in- two codes offers many options of graphical representation of formation needed for the projections is available e.g. from structures, adding translations, clipping fragments etc. Three- the Bilbao Crystallographic Server,139 and the technical de- dimensional spatial functions (e.g., charge density, local den- tails are explained in the documentation included in the tools. sity of states integrated throughout the chosen energy range), The vibent tool performs a straightforward calculation computed by S IESTA on a real-space grid. Tools are provided (see, e.g., Sec. II.C in Ref. 140, or Sec. 5.3 in Ref. 141) of S IESTA 19 Be32Zn64Se96 ; Zn contribution 0.05 "Ni4": molecular dynamics at 600 K projected vibration DOS for q = ( 0 0 qz) -1 qz (Bohr ) "Ni4": frozen phonon calculation 0 TO 002-reduced LO -0.05 Velocity autocorrelation function Be Ni 0 100 200 300 400 500 600 700 800 900 1000 0.00 H Mo MD steps Zn 0.01 O 0.006 0.02 Se 0.03 0.004 0.04 cos - VCF 0.05 0.002 0.06 0.07 0.08 0 0 500 1000 1500 2000 2500 3000 3500 4000 -1 0.09 Frequency (cm ) 0 500 1000 1500 2000 2500 3000 3500 4000 0.10 -1 Frequency (cm ) 0.11 0.12 0.13 FIG. 18. Vibration properties of “Ni4 ” molecular magnet, 0.14 0.15 [Mo12 O30 (µ2 -OH)10 H2 {Ni(H2 O)3 }4 ]·14H2 O. Left panel: density 0.16 of modes from frozen phonon calculation; right panel: velocity au- 0.17 0.18 tocorrelation function, its Fourier transform and hence resulting den- 0.19 sity of vibration modes. 0 100 200 300 400 -1 wave number (cm ) MD steps at 600 K) is shown in Fig. 18 in comparison with FIG. 16. Left panel: a 192-at. quasirandom supercell representative frozen phonon results, revealing similarities of the spectra ob- for the Be1/3 Zn2/3 Se solid solution; right panel: density of modes within the frequency range of Zn-Se vibrations, extracted with phdos tained. and projected onto different values of qz and different polarisations, parallel (labelled LO) and perpendicular (TO) to q. These results were partially shown in Fig. 4 of Ref. 138 and discussed in that work. 4. Optical properties of finite systems: linear response TDDFT starting from Siesta orbitals Cu2ZnSnS4: 64-at. supercell, total density of modes Vibrational properties of point defects in CZTS: doped minus perfect Arb. units, per supercell Σ (64at.) Cu (16 at.) 0 200 400 600 800 1000 The S IESTA package offers at least two ways of obtain- Zn (8 at.) Vibrational free energy (eV) Sn (8 at.) 0.1 Temperature (K) ing optical properties of finite systems. The first way uses S (32 at.) 0 real-time TD-DFT propagation by applying an external elec- -0.1 -0.2 tric field with a simple time dependence (e.g., a Heaviside step-function)74 . The second way is by computing the non- +0 ZnCu(Sn-pl.) +0 -0.3 ZnCu(Zn-pl.) 0 50 100 150 200 250 300 350 400 interacting dielectric function1,2 . Both methods are imple- +0 CuZn +0 ZnCu(Sn-pl.) -0.4 +0 ZnSn _ +0 +0 ZnCu(Zn-pl.) -0.5 |_| Cu(Sn-pl.) _ +0 |_| Cu(Zn-pl.) mented in S IESTA and can be employed without any external Difference from perfect _ +0 |_| Zn tools. However, they are limited in different aspects. The Vibrational entropy (meV/K) +0 CuZn _ +0 |_| S 0.4 ZnSn +0 _ +0 non-interacting dielectric function often underestimates the |_| Cu(Sn-pl.) 0.2 _ +0 |_| Cu(Zn-pl.) HOMO-LUMO gaps and calls for the use of the phenomeno- 0 _ +0 |_| Zn logical scissor-shift operator. Real-time propagation makes -0.2 cumbersome the analysis of the optical response properties in _ +0 |_| S -0.4 the frequency domain. Furthermore, the frequency resolution 0 50 100 150 200 250 300 350 400 0 200 400 600 800 1000 Frequency (cm ) -1 Temperature (K) scales with the duration of the real-time simulation. Thus, ac- curate spectra require long simulations. FIG. 17. Vibration properties of Cu2 ZnSnS4 (CZTS) with substi- Fortunately, there are two efficient implementations of tutional impurities, used in Ref. 140. Left panel: densities of modes linear-response TDDFT that use the Kohn-Sham orbitals from (extracted with phdos); right panel: vibration contributions to the S IESTA as a starting point and are available for the open- free energy and entropy (calculated with vibent). Adapted from source community144,145 . In both packages, the linear den- Fig. 5.3 and 5.4 of Ref. 141. sity response δ n(r, ω) is obtained directly in the frequency domain which makes straightforward the analysis of derived properties. However, there are differences between both im- temperature-dependent vibration contributions to the free en- plementations on the construction of the auxiliary basis nec- ergy and entropy – see Fig. 17 as an example. The necessary essary to expand the orbital products. These differences can input information is the vibration spectrum, originating from severely affect the computational cost of the calculation. the Vibra frozen phonon calculation on a sufficiently large The linear-response TDDFT is built on the concept of the system. induced electronic density δ n(r, ω) in response to a small per- The velcf tool calculates the velocity autocorrelation turbation of the external potential δVext (r, ω). The integral function and its Fourier transform from a (presumably suf- operator connecting δ n(r, ω) to δVext (r, ω) is the interacting ficiently long) molecular dynamics (MD) history, recorded in density response function χ(r, r 0 , ω). By virtue of the KS the .MD or .ANI file. This technique142 can be used to obtain equations, χ(r, r 0 , ω) can be connected to the non-interacting phonon frequencies, and was applied along with a S IESTA cal- density response function χ0 (r, r 0 , ω)146 with a Dyson equa- culation in Ref. 143. An example of such simulation (1000 tion S IESTA 20 2 χ(ω) = χ0 (ω) + χ0 (ω)Kχ(ω), (29) Ag561 1.6 Ag309 Clusters where the interaction kernel K(r, r 0 ) contains the bare Ag147 σopt (nm2) Ag55 Coulomb interaction and the so-called exchange and correla- 1.2 Ag13 tion kernel Kxc , which is a known operator for simple func- tionals like LDA and GGA. The non-interacting response 0.8 function χ0 (r, r 0 , ω) can be expressed as a sum over electron- hole excitations within the basis formed by the KS orbitals 0.4 Ψn (r)145–147 0 2 3 4 5 6 7 8 Ψn (r)Ψm (r)Ψm (r 0 )Ψ n (r 0 ) χ0 (r, r 0 , ω) = ∑( fn − fm ) , ω (eV) nm ω − Em + En (30) where fn are occupations of the KS orbitals and En are their FIG. 19. The absorption cross sections of silver clusters of icosa- hedral shape. One can recognize sharp surface-plasmon resonances energies. around 3–4 eV and a broad resonance at 6–7 eV. The optical polarizability R tensor α(ω) is related to the in- duced density by α(ω) = rδ n(r, ω)dr or alternatively ZZ α(ω) = rχ0 (r, r 0 , ω)δVs (r 0 , ω) drdr 0 , (31) where due to Eq. (29) and using the dipole approximation for the electron-photon coupling, the screened effective perturba- tion δVs (r 0 , ω) satisfies the linear integral equation (I − Kχ0 (ω))δVs (ω) = r. (32) The efficiency of the methods presented in References 144 and 145 comes from solving iteratively Eq. (32) for δVs (ω) instead of using standard matrix inversion to obtain χ(ω) FIG. 20. The isosurfaces of density change Re(δ n(r, ω)) of the from Eq. (29). Once δVs (ω) is known, Eqs. (30) and (31) Ag147 cluster close to the frequency of the surface-plasmon reso- allow the computation of the optical properties of the system. nance of the cluster (3.4 eV). Furthermore, it is also possible to perform different types of analysis. For example, it is easy to partition the polarizability tensor α(ω) in terms of electron-hole contributions147,148 due 5. Thermal transport by the AEMD method to existence of the sum over the electron-hole pairs in Eq. (30). Similarly, one can achieve other types of Mulliken-like anal- The approach to equilibrium molecular dynamics (AEMD) ysis145,147,149 of the optical polarizability tensor α(ω) or the method153 has been implemented to obtain the thermal con- induced density δ n(r, ω). ductivity. In the first stage of the method, the system is de- The Python implementation of linear response TDDFT in composed in two different regions, each one equilibrated to a the PySCF-NAO package as described in Ref. 145 is con- different initial temperature (canonical run with Bose, or An- venient to use and rather potent. It is capable of computing neal MD). Then, a microcanonical run (Verlet) is carried out the optical properties of compact metallic objects containing for the whole system, and the average temperature of each up to several hundreds of atoms147,150,151 . For example, we subsystem is monitored. This temperature transient regime is were able to track down the different size-dependence of the then used to extract the thermal conductivity from the exact plasmon resonance in sodium and silver clusters due to the solution of the heat transport equation.154 screening effect of silver d-orbitals in the latter case152 . In those calculations, using an optimized version that incorpo- rates some additional memory-saving features not present in the currently distributed version of PySCF-NAO, icosahedral 6. Core level shifts silver and sodium clusters containing up to 5043 atoms were studied. Core-level shifts can serve to analyze changes in the lo- In Figures 19 and 20, we show the photo-absorption cross cal and chemical environment of atoms of a given species. sections of a series of compact silver clusters147 and the real Density-functional-theory calculations have proved to be part of induced density change in the cluster Ag147 close to its quite useful in complementing the experimental information, surface-plasmon frequency (3.4 eV), respectively. which is sometimes hard to interpret. Two schemes have been S IESTA 21 implemented in S IESTA for the calculation of core-level shifts they have enabled a much simpler bookkeeping of the data within a pseudopotential approach155 . structures needed for a richer control of molecular-mechanics In the so-called initial-state approximation the electronic and scf iterations. relaxation in the presence of the core hole is neglected, and Regarding performance-oriented developments, in the re- the photo-electron’s binding energy is directly related to the cent past we have implemented a mixed MPI/OpenMP pro- eigenvalue of the core level. A pseudopotential calculation gramming model, which allows, for suitable systems, to better obviously cannot compute the latter, but differences in core balance arithmetic intensity and communications needs. The eigenvalues in different environments can be estimated by the deployment of this model is more advanced in the TranSiesta changes in the expectation value of the crystal potential using module, and significant speedups have been obtained for large the core state’s atomic wavefunctions ψnlm at different sites. systems. These can be extracted from the matrix elements Some of the above software-engineering developments have been enabled and strengthened by the participation of Z S IESTA in a number of international partnerships, notably the mm0 0 V = d 3 r (ψnlm (~r −~τ))∗V (~r)ψnlm (~r −~τ) (33) MaX (Materials at the eXascale) EU center of excellence158 and the Electronic Structure Library initiative159 . The “sepa- with a further step of averaging to remove the splittings stem- ration of concerns” described above in the context of modu- ming from the loss of spherical symmetry. larization is an example of the so-called “open-innovation” In the final-state approximation, the relaxation is explicitly paradigm, at the foundation of the ESL strategy for code taken into account, and the experimental shifts (measured via reusability, and is also a cornerstone of MaX’s efforts to the kinetic energy of an exiting electron) are correlated with achieve exascale-readiness for its flagship materials science the differences in the energy of the crystal with a “core-hole” codes (with S IESTA among them): performance-enhancement in different sites. For this, a special pseudopotential with a efforts are to be focused on relevant domain-specific modules. missing core electron has to be generated, and a full S IESTA A number of modules from S IESTA have been turned calculation is needed for each different site. into stand-alone libraries which now feature in the ESL: The implemented methodology has been used to study, for libGridXC for exchange and correlation calculations, example, the shifts induced by hydrogen bonding in organic libPSML as a handler of PSML files, xmlf90 for general molecules156 . purpose handling of XML files, etc. Conversely, S IESTA uses some of the libraries offered by the ESL, notably the ELSI library of electronic-structure solvers mentioned in O. Software-engineering advances and partnerships Sect. III G 2, whose development, including its API design and internal data organization, has been in turn influenced by contributions and feedback from the S IESTA project, among The traditional development model for scientific codes in others. There are also plans to incorporate the PSolver li- academic settings has been typically based on multiple con- brary160 for the solution of the Poisson problem, a contribu- tributions with various levels of programming competence, tion to the ESL from the BigDFT project. and with very little time to plan ahead in the face of press- We should mention that the renewed dynamism of S IESTA ing scientific demands. S IESTA has been no exception, and development and the advances made possible by the interac- has grown in features and complexity over the years. It is tion with community initiatives are both a blessing and a chal- very important to keep complexity under control, or else a lenge. It is non-trivial, for example, to handle the building project becomes un-maintainable and cannot survive. It is not process of a code that relies on a number of different external simple, however, to balance the need of incorporation of new libraries, programming models, and special features such as features, and the need to increase the computing performance the embedded Lua interpreter. Luckily, as will be discussed in a landscape of constantly evolving hardware and program- in Sect. V, these are issues that are being addressed in wider ming models. One essential route is modularization, which contexts, and S IESTA is well placed to take advantage of it. allows the separation of concerns at various levels. In the context of a code like S IESTA , this means that the scientific ideas and algorithms should be handled at a high level, call- IV. APPLICATIONS ing on lower-level modules for specific functionality (domain- specific libraries, mathematical libraries, communication pro- tocols, etc). These lower-level modules can hopefully be re- We present here a few showcase applications that illustrate used by different codes and, most importantly, can be focused the capabilities of S IESTA, in breadth, efficiency, and accu- on by highly-skilled programmers for optimization on rele- racy. vant architectures. Another important method of taming complexity involves the streamlining of the data structures of the code. This is an A. 4 terminal NEGF on germanium surface ongoing process (see Sect. V), but has already taken a very significant step by the introduction of reference-counted data Breakthrough simulations using the new multi-terminal im- structures. They build on a well-known and not particularly plementation on T RAN S IESTA were fundamental to elucidate advanced technique of memory-handling157 , but in S IESTA the electronic transport mechanism on a novel and complex S IESTA 22 B. Novel topological phases in ferroelectric materials In material systems with several interacting degrees of free- dom (such as spin, charge and lattice distortions), the com- plex interplay between these factors can give rise to exotic phases. A prototypical example are the superlattices of al- ternating lead titanate and strontium titanate layers. Simula- tions on such PbTiO3 /SrTiO3 heterostructures, consisting on n unit cells of PbTiO3 and n unit cells of SrTiO3 stacked along the [001] direction were carried out with S IESTA. As a func- tion of the periodicity, the superlattices undergo a phase tran- sition from a monodomain configuration (small periodicity, n . 3 − 4) with a normal component of the polarization that is FIG. 21. First-principles transport simulations for the two-probe ex- preserved throughout the structure, to a multidomain config- periments. a) Representation of the four-terminal setup. The elec- uration (large periodicity, n & 3 − 4) with alternating up and trode regions are highlighted by blue boxes, two of them located at down domains.161 In order to further reduce the electrostatic each Ge(001)-c(4 × 2) slab terminations (leads left and right) and the energy costs, the local dipoles within the PbTiO3 layer con- other two at each Au model tip (leads tip1 and tip2). The 50 Ge atoms closest to each tip were allowed to fully relax, adapted from tinuously rotate forming a sequence of clock-wise/counter- Ref. 101. clockwise array of vortices along the [100] direction. The theoretical predictions, done with S IESTA162 after the relax- ation of supercells of up to 1000 atoms, were experimentally confirmed five years later by atomic-scale mapping of the po- lar atomic displacements by scanning transmission electron microscopy163 (Fig. 22) Moreover, the appearance of an axial component of the polarization pointing in the direction of the vortices make the systems chiral and optically active, as lately confirmed by circular dichroism experiments164 . experiment.101 For the first time a two-probe scanning tun- neling microscopy/spectroscopy (STM/STS) with probes op- erating in tunneling conditions over the same atomic-scale C. 1D and 2D systems system was used to extract detailed information of in-plane electronic transport. The addressed system was the recon- S IESTA is particularly well suited to study low dimen- structed (001) surface of germanium, where electrons injected sional nanostructures, such as 1D and 2D systems where a from one STM tip at a position determined with atomic pre- large vacuum region is needed within the simulation cell. cision were collected at the same Ge dimer row at a distance When, in addition, a large number of atoms is required to as short as 30 nm. The experiment was theoretically mod- study particular physical effects is where S IESTA could ex- eled by a system composed of a twelve-layer Ge(001)-c(4×2) cel with respect to other methods. There is extensive liter- slab contacted by Au tips oriented along the (100) direction ature on simulations of graphene and other exfoliated mate- (Fig. 21). On this self-consistent 4-terminal treatment, two rials, where the properties of point defects, edges or grain Ge electrodes were connected at each slab termination and boundaries are of much relevance. To list a few examples, other two at the Au model tips. The whole system was de- the magnetic properties of impurities,165,166 and edges167 , fined by 4924 atoms (36442 atomic orbitals), in a super-cell but also electronic properties, including transport charac- 3 of dimensions ∼ 32 × 160 × 80 Å , and where 5 different tip- teristics, in grain boundaries168,169 , ribbons170 , nanoporous to-sample distances were considered. Besides the large di- graphene171 , large graphene flakes68,172 , or the effect of sub- mensions of the system, another important challenge of such strates173 . Other materials, such as mono- and multi-layered simulation was the level alignment between the metallic and dichalcogenides174,175 or phosphorene176,177 , are also being semiconducting leads and the scattering region, for which a widely studied, including optical properties in nanoflakes with method had to be devised. A remarkable agreement was found up to a few thousand atoms.178 . between the calculated transmission function and the experi- mental transconductance spectra, allowing the identification and assignment of the observed resonances to transport chan- 1. CDWs nels existing along the surface Ge dimer rows. Moreover, the simulations elucidated the transport directionality of the in- A number of recent studies on charge density waves (CDW) jected hot electrons, revealing a transition from 2D to quasi- in low dimensional materials illustrates the impressive accu- 1D coherent transport regime as a function of the carrier’s en- racy that can be obtained with S IESTA for systems with very ergy. This work shows that complex experiment setups com- subtle electronic structures.179 For example, in 2H-NbSe2 bined with advanced calculations can provide new insights S IESTA calculations were able to predict the existence of six into transport properties at the nanoscale. different atomic structures within a narrow energy range of a Panorama S IESTA 23 EL EC TR ÓN I C A experiments), providing, for the first time, a quantitative evi- Condensadores dence of the weak electron-phonon coupling scenario for the Peierls transition. con capacidad negativa Un estudio encuentra el origen microscópico de la capacidad D. negativa, Siesta in biology: pilin proteins as conductors una exótica propiedad electrónica que aparece en ciertos materiales. El hallazgo augura el diseño de nuevos transistores más eficientes S IESTA’s efficiency and the clear band gaps of Pablo García-Fernández y Javier Junquera biomolecules in general have made molecular biology a very suitable field for S IESTA since the beginning,182 and have stimulated targeted developments of the code for the field, such as QM/MM.183,184 An interesting illustration of its suitability in an all-quantum biological problem is the study of the electrostatics around the pilin protein in aqueous solution.185 The pilin considered here is the main protein in the pili (external filaments) of the geobacter sulfurreducens bacterium, which have been shown to be able to transmit electronic current, allowing the microbe to feed by remote redox reactions on ferrous mineral particles in the soil. As a nanowire designed by natural evolution, understanding the mechanism for charge transport is of obvious interest. Peculiar to this protein is the fact that its main alpha helix, the main feature of this elongated protein, is singly oriented, that is, there is no back alpha helix (as in a common hair- pin configuration) that would counter the polarization of the single alpha helix: In an alpha helix all peptide-bond dipoles point in the same direction along the axis of the helix, which, FIG. 22. Top panel: local polarization profile of polydomain struc- in solid-state parlance, represents a polarization, with clear tures in (PbTiO3 )n /(SrTiO3 )n with n=6 obtained from an atomic re- electrostatic implications. Indeed, a DFT calculation of the laxation with S IESTA. The PbTiO3 and SrTiO3 are depicted as molecule in vacuum shows a well defined electrostatic poten- grey and white regions respectively. Clockwise and counterclock- tial ramp along the protein, which tends to close the effective wise vortices within the PbTiO3 are clearly visible. Red dashed band gap. The question is then, how does an aqueous environ- square in the SrTiO3 layers mark the position where antivortices are formed. Reprinted with permission from Aguado-Puente and Jun- ment affect this depolarizing field. quera 162 Phys. REMOLINOS Rev. B 85,Mapa DE POLARIZACIÓN: microscópico 184105 (2012). de laBottom polarización (flechasex- panel: Longde titanato amarillas) en capas molecular de plomomechanics (MM) simulations were per- y titanato deobservation perimental estroncio. La inusual reacción de estos materiales of vortex–antivortex ante unin structures campo externo genera a cross- formed for the zonas donde protein in a suitable solution of NaCl at a con- la capacidad CORTESÍA DE RAmAmOORThy RAmESh, UNIVERSIDAD DE CALIFORNIA EN BERKELEy eléctrica es negativa, un fenómeno considerado imposible hasta hace pocos años. sectional high-resolution scanning transmission electron microscopy centration of 0.1 M. The protein’s residues had charge states image with an overlay of the polar displacement vectors for a corresponding to pH = 7, and the MM field was validated T odos nos los aparatos (SrTiO rodean se basan en tex–antivortex electrónicos 3 )10 /(PbTiO el uso pairs que deen, Walter 3 )10 superlattice, tran- en isdepresent in los each showingBrattain Laboratorios thaty an William Bell, fueCourtesy PbTiO3 layer. of vor- según arrayShockley sistores: diminutos dispositivos de meca- en 1956 con el premio Nobel de física y cada dos años. la cual with SelIESTA número de transistoresin calculations reconocidoof un circuito se duplica aproximadamente 3 3 en vacuum (944-atom dynamic re- laxation in a 104.4 Å box). The wet system contained 4580 R. Ramesh, adapted from Ref. 163. nismo aparentemente simple. Un transis- supuso el inicio de la revolución tecnoló- Enatoms, and the la actualidad, lasstatistical dimensionesaverage de of the electrostatic potential tor puede entenderse como un pequeño gica de la que disfrutamos hoy. estosaround componentes electrónicos(see the molecule se sitúan Fig. 23) was obtained from a sample interruptor con dos posiciones: «encendi- Desde sus comienzos, la industria entre los 15 y los 50 nanómetros, una esca- of full S IESTA calculations of statistically independent snap- do» (permite la circulación de corriente) microelectrónica se ha obsesionado con la en la que los efectos cuánticos comien- few meV, all of them compatible with the experimental 3×3 y «apagado» (no la permite). Una com- integrar un número cada vez mayor de zan ashots, taken Como ser relevantes. everyconsecuencia, 50 ps during the last 0.5 ns of the simula- CDW binación modulation. adecuada Careful de ellos permite im- analysis transistoresofentheoretical and carrera los circuitos. Esta exper- un procesador tion. moderno puede llegar a imental plementar STM cualquier imageslógica. operación for different bias potentials por la miniaturización allowed se resume en lato albergar miles de millones de transistores. Su descubrimiento, identify two efectuado años conocida en losstructures of these thatleycan de Moore, coexistenunciada in theensame1965 Sin embargo, ese aumento en el nú- por John Bar- por Moore, cofundador de Intel, y mero de componentes por circuito tiene cuarenta del siglo image.180pasado In a different work,181Gordon the temperature dependency of the electronic Lindhard response function in blue bronze K0.3 MoO 12 INVESTIGACIÓN 3 was agosto Y CIENCIA, studied. 2019 This system has a rather complex monoclinic structure, with twenty formula units per unit cell where MoO6 octahedra form chains along one direction (b- axis). The Lindhard function shows well decoupled sharp responses that correspond to intra- and interband Fermi sur- face nesting. By fitting these peaks one can obtain the coher- FIG. 23. Colour coded electrostatic potential on a plane cutting along ence length of the fluctuating 1D electron-hole pair (that de- the main axis of the geobacter sulfurreducens pilin molecule in wet termines the length scale of the experimental intrachain CDW conditions. A perspective ball rendering of the atomic strucuture of correlations), and the intrachain modulation of the response the protein is superposed. For the meaning and details on this Figure (that determines the shape of the Kohn anomaly measured in see Ref. 185 (Figure courtesy of Gustavo T. Feliciano). S IESTA 24 Fig. 23 shows how the aqueous environment kills the quite variables and encapsulate them into objects or derived types homogeneous potential ramp along the protein axis that ap- associated to particular configurations and stages of the calcu- pears in vacuum and replaces it with long-wave-length slow, lations. This encapsulation will be matched by a streamlining but quite significant fluctuations. The gap remains sizeable, of the input/output operations. This work will open the door and coherent transport is not likely. However, the frontier or- to the creation of complex workflows leveraging the strengths bitals evolve in a very suggestive way for enhanced diffusive of various codes. electron transport.185 Accelerated hybrid architectures (including, for example, GPUs) are very likely going to feature prominently in the up- coming exascale machines. In the case of S IESTA, the data E. Use of Siesta in other fields indirection associated to the handling of sparse matrices lim- its the acceleration possibilities of the section of the code that Although an exhaustive summary of all the recent results builds the Hamiltonian and overlap matrices, but the solver obtained with S IESTA is out of the scope of this work, we stage is more amenable to porting, and in fact several solver would like to point the attention of the reader to a sample of libraries used by S IESTA are being enhanced to offer GPU recent reviews in various fields in which the program is fea- support, as mentioned in Sec. III G 2. tured. These cover biological sciences186 (including interac- Modularization and the use of new programming models tion between organic and inorganic materials187,188 ), geology cause an increase in the complexity of the building and de- and materials under high-pressure189 , isotopic fractionation ployment of the code. We will leverage the ESL bundle, predictions for Martian geochemistry190 , the engineering of created to facilitate the use of the modules in the ESL col- typical core structural materials used in nuclear reactors,191 lection, to streamline S IESTA’s building process, and explore or even in astrophysical and atmospheric systems192 . The containerization as an option for deployment of the code. reactivity of metallic nanoparticles for catalysis was treated The “pseudopotential barrier to entry” has been lowered by by Viñes, Gomes, and Illas 193 , and the role of S IESTA in the the availability of curated databases supporting the PSML for- computation of the kinetic and dynamics of catalytic reaction mat. Basis sets are a perennial challenge, but new tools and at surfaces (including adsorption and desorption of reactants ideas are being explored to provide users with appropriate ba- or products) was explored in Chapter 8 of Ref. 194 by Catapan sis sets: High-throughput workflows for optimization; "tiers" and coworkers. of quality/cost, but perhaps not just of a simple “periodic ta- ble” form, as offered by other codes (e.g., FHI-aims13 ), but with a possible dependence on an approximate characteriza- V. FUTURE EVOLUTION tion of the chemical environment in which a given atom finds itself. Work on enhancing S IESTA’s capabilities, performance, Complementary to the underlying basis-set optimization and robustness is continuing, driven by a good number of de- that focuses on providing an adequate variational freedom, an velopers and collaborators. A mature and flexible develop- on-the-fly contraction of the basis set, which results in a set ment platform and practices are essential to keep them pro- of lower-cardinality adapted to the description of the occu- ductive. Our recent platform changes have forced develop- pied subspace can be exploited for increased efficiency. This ers to shift workflows twice in the past four years. Through is particularly relevant for FOE methods (see Sect. III G 1, in the changes we have learned a lot but also spent a signifi- which the number of polynomial terms depends on the extent cant amount of time on ensuring S IESTA’s continuous devel- of the spectrum. opment. At the current state we believe we have stabilized the The original claim to fame of S IESTA was based on its development platform on GitLab while we will add more in- linear-scaling solver. We are in the process of a re-design of tegrated development features in the coming years, e.g. con- the O(N) code with a new, more efficient backend, based on tinuous integration (CI) and source code checks. Using CI the DBCSR library for handling distributed block-sparse ma- will also enable easier code-style checks to conform to coding trices195,196 with the MatrixSwitch library82 acting as an in- standards. We hope that our open-platform initiative will keep termediary interface between it and high-level physical ideas external contributions coming into the program. and algorithms. A connection between the internal S IESTA Our basic-development plans include also refactoring, ap- formats and MatrixSwitch itself has been recently provided, parently unexciting but essential to streamline the code base using initially the cubic-scaling libOMM library197 as a test to enable further implementations. Also, we foresee a change bed, hence still using a dense coefficient matrix, as it corre- in the release model, moving away from coexisting long-lived sponds to the case without localization constraints in the solu- release branches whose maintenance takes up a lot of time, tion of the electronic-structure problem. The implementation and offering instead more frequent and short-maintenance re- of a sparse coefficient matrix will make it possible to perform leases. efficient O(N) calculations. The computational effort can be We plan to exploit the idea of modularization, continuing further reduced through the analysis of sparsity of the Hamil- the abstraction of relevant reusable pieces, but also dealing tonian and overlap matrices and their re-organization in the with a higher-level, exposing the core electronic-structure ca- block-compressed sparse form. pabilities of S IESTA to other programs. It will be necessary to Other developments in the pipeline are linear-response cal- redesign some of the internal data structures to remove global culations for arbitrary distortions, electronic transport calcu- S IESTA 25 lations with spin-orbit coupling, thermal transport with the 7 O. F. Sankey and D. J. Niklewski, “Ab initio multicenter tight-binding Green-Kubo formalism, as described in Ref. 198, a redesign model for molecular-dynamics simulations and other applications in cova- of the molecular dynamics subsystem, and the development lent systems,” Phys. Rev. B 40, 3979–3995 (1989). 8 E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler, of workflows for the generation of data for S CALE -U P. “Linear-scaling ab-initio calculations for large and complex systems,” Phys. Status Solidi (b) 215, 809–817 (1999). 9 J. Junquera, O. Paz, D. Sánchez-Portal, and E. Artacho, “Numerical ACKNOWLEDGMENTS atomic orbitals for linear-scaling calculations,” Phys. Rev. B 64, 235111 (2001). 10 E. Anglada, J. M. Soler, J. Junquera, and E. Artacho, “Systematic genera- S IESTA development has been historically supported by tion of finite-range atomic basis sets for linear-scaling calculations,” Phys. different Spanish National Plan projects: MEC-DGES- Rev. B 66, 205101 (2002). 11 http://www.openmx-square.org/. PB95-0202, MCyT-BFM2000-1312, MEC-BFM2003-03372, 12 https://www.synopsys.com/silicon/quantumatk.html. FIS2006-12117, FIS2009-12721, FIS2012-37549, FIS2015- 13 V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, 64886-P, and RTC-2016-5681-7, the latter one together with and M. Scheffler, “Ab initio molecular simulations with numeric atom- Simune Atomistics Ltd. Currently, we thank financial support centered orbitals,” Comput. Phys. Commun. 180, 2175 – 2196 (2009). 14 A. Carreras, S. Conejeros, A. Camón, A. García, N. Casañ-Pastor, P. Ale- from the Spanish Ministry of Science, Innovation and Univer- many, and E. Canadell, “Charge delocalization, oxidation states, and sil- sities through the grant No. PGC2018-096955-B. ver mobility in the mixed silver–copper oxide agcuo2,” Inorg. Chem. 58, We acknowledge the Severo Ochoa Centers of Ex- 7026–7035 (2019). 15 M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, cellence Program under Grants No. SEV-2015-0496 (ICMAB), and SEV-2017-0706 (ICN2), the GenCat Grant No. “Density-functional method for nonequilibrium electron transport,” Phys. Rev. B 65, 165401 (2002). 2017SGR1506, and the European Union MaX Center of Ex- 16 F. Corsetti, E. Artacho, J. M. Soler, S. S. Alexandre, and M.-V. Fernández- cellence (EU-H2020 Grant No. 824143). Serra, “Room temperature compressibility and diffusivity of liquid water P.G.-F. acknowledges support from Ramón y Cajal Grant from first principles,” J. Chem. Phys. 139, 194502 (2013). 17 X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, No. RyC-2013-12515. J.I.C acknowledges RTI2018-097895- B-C41. P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, T. Deutsch, L. Genovese, P. Ghosez, M. Giantomassi, S. Goedecker, D. Hamann, R.C. acknowledges to the European Union’s Horizon P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, 2020 research and innovation program under the Marie M. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, D. San- Skłodoswka–Curie grant agreement no. 665919. galli, R. Shaltaf, M. Torrent, M. Verstraete, G. Zerah, and J. Zwanziger, D.S.P, P.K, and P.B acknowledge MAT2016-78293-C6, “Abinit: First-principles approach to material and nanosystem properties,” Comput. Phys. Commun. 180, 2582 – 2615 (2009). FET-Open No. 863098, and UPV-EHU Grant IT1246-19. 18 See: https://www.gnu.org/licenses/gpl-3.0.html. V. Yu was supported by a MolSSI fellowship (U.S. NSF 19 See: https://launchpad.net/siesta. 20 See: https://gitlab.com/siesta-project. award 1547580), and the ELSI development (V.B.,V.Yu) by 21 A. García, M. J. Verstraete, Y. Pouillon, and J. Junquera, “The psml format NSF award 1450280. We also acknowledge Honghui Shang and library for norm-conserving pseudopotential data curation and inter- and Xinming Qin for giving us access to the H ONPAS code, operability,” Comput. Phys. Commun. 227, 51 – 71 (2018). where a preliminary version of the hybrid functionals support 22 See: https://siesta-project.github.io/psml-docs, accessed described here was implemented. November 2019. 23 D. R. Hamann, “Optimized norm-conserving vanderbilt pseudopotentials,” We are indebted to other contributors to the S IESTA project, whose names can be seen in the file in the Phys. Rev. B 88, 085117 (2013). 24 ATOM code for the generation of norm-conserving pseudopotentials. The Docs/Contributors.txt file of the S IESTA distribution, version maintained by the S IESTA project can be accessed at http: and we thank those, too many to list, contributing fixes, com- //icmab.es/siesta/Pseudopotentials/index.html. An alterna- ments, clarifications, and documentation for the code. tive version is available at http://bohr.inesc-mn.pt/~jlm/pseudo. The data that support the findings of this study are available html. (Accessed July 2017). 25 M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. from the corresponding author upon reasonable request. Hamann, X. Gonze, and G. M. Rignanese, “The PSEUDODOJO: Train- 1 J. M. Soler, ing and grading a 85 element optimized norm-conserving pseudopotential E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and table,” Comput. Phys. Commun. 226, 39–54 (2018). D. Sánchez-Portal, “The SIESTA method for ab initio order-n materials 26 See: http://www.pseudo-dojo.org. simulation,” J. Phys.: Condens. Matter 14, 2745–2779 (2002). 27 X. Gonze, F. Jollet, F. Abreu Araujo, D. Adams, B. Amadon, T. Applen- 2 E. Artacho, E. Anglada, O. Diéguez, J. D. Gale, A. García, J. Junquera, court, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, E. Bousquet, R. M. Martin, P. Ordejón, J. M. Pruneda, D. Sánchez-Portal, and J. M. F. Bruneval, D. Caliste, M. Côté, F. Dahm, F. Da Pieve, M. Delaveau, Soler, “The SIESTA method; developments and applicability,” J. Phys.: M. Di Gennaro, B. Dorado, C. Espejo, G. Geneste, L. Genovese, Condens. Matter 20, 064208 (2008). A. Gerossier, M. Giantomassi, Y. Gillet, D. Hamann, L. He, G. Jo- 3 G. Galli, “Linear scaling methods for electronic structure calculations and mard, J. Laflamme Janssen, S. Le Roux, A. Levitt, A. Lherbier, F. Liu, quantum molecular dynamics simulations,” Curr. Opin. Solid State Mater. I. Lukačević, A. Martin, C. Martins, M. Oliveira, S. Poncé, Y. Pouil- Sci. 1, 864 – 874 (1996). lon, T. Rangel, G.-M. Rignanese, A. Romero, B. Rousseau, O. Rubel, 4 S. Goedecker, “Linear scaling electronic structure methods,” Rev. Mod. A. Shukri, M. Stankovski, M. Torrent, M. Van Setten, B. Van Tro- Phys. 71, 1085–1123 (1999). eye, M. Verstraete, D. Waroquiers, J. Wiktor, B. Xu, A. Zhou, and 5 P. Ordejón, E. Artacho, and J. M. Soler, “Self-consistent order-n density- J. Zwanziger, “Recent developments in the ABINIT software package,” functional calculations for very large systems,” Phys. Rev. B 53, R10441– Comput. Phys. Commun. 205, 106–131 (2016). R10444 (1996). 28 V. I. Anisimov, J. Zaanen, and O. K. Andersen, “Band theory and mott in- 6 D. Sánchez-Portal, P. Ordejón, E. Artacho, and J. M. Soler, “Density- sulators: Hubbard u instead of stoner i,” Phys. Rev. B 44, 943–954 (1991). functional method for very large systems with lcao basis sets,” Int. J. Quan- 29 B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, “Hubbard- tum Chem. 65, 453–461 (1997). S IESTA 26 corrected dft energy functionals: The lda+u description of correlated sys- 55 L. A. Hemstreet, C. Y. Fong, and J. S. Nelson, “First-principles calcula- tems,” Int. J. Quantum Chem. 114, 14–49 (2014). tions of spin-orbit splittings in solids using nonlocal separable pseudopo- 30 S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. tentials,” Phys. Rev. B 47, 4238 (1993). Sutton, “Electron-energy-loss spectra and the structural stability of nickel 56 F. Zirkelbach, P.-Y. Prodhomme, P. Han, R. Cherian, and G. Bester, oxide: An lsda+u study,” Phys. Rev. B 57, 1505–1509 (1998). “Large-scale atomic effective pseudopotential program including an effi- 31 M. Cococcioni and S. de Gironcoli, “Linear response approach to the cal- cient spin-orbit coupling treatment in real space,” Phys. Rev. B 91, 075119 culation of the effective interaction parameters in the LDA + U method,” (2015). Phys. Rev. B 71, 035105 (2005). 57 V. W. zhe Yu, F. Corsetti, A. García, W. P. Huhn, M. Jacquelin, W. Jia, 32 R. J. Powell and W. E. Spicer, “Optical properties of nio and coo,” Phys. B. Lange, L. Lin, J. Lu, W. Mi, A. Seifitokaldani, Álvaro Vázquez- Rev. B 2, 2182–2193 (1970). Mayagoitia, C. Yang, H. Yang, and V. Blum, “Elsi: A unified software 33 H. A. Alperin, J. Phys. Soc. Jpn. Suppl. B 17, 12 (1962). interface for kohn-sham electronic structure solvers,” Comput. Phys. Com- 34 A. K. Cheetham and D. A. O. Hope, “Magnetic ordering and exchange mun. 222, 267 – 285 (2018). effects in the antiferromagnetic solid solutions mnx ni1−x O,” Phys. Rev. B 58 V. W. zhe Yu, C. Campos, W. Dawson, A. García, V. Havu, B. Hourahine, 27, 6964–6967 (1983). W. P. Huhn, M. Jacquelin, W. Jia, M. Keçeli, R. Laasner, Y. Li, L. Lin, 35 M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, J. Lu, J. Moussa, J. E. Roman, Álvaro Vázquez-Mayagoitia, C. Yang, and “Van der waals density functional for general geometries,” Phys. Rev. Lett. V. Blum, “Elsi – an open infrastructure for electronic structure solvers,” 92, 246401 (2004). (2019), arXiv:1912.13403 [physics.comp-ph]. 36 K. Berland and P. Hyldgaard, “Exchange functional that tests the robust- 59 J. Choi, J. Demmel, I. Dhillon, J. Dongarra, S. Ostrouchov, A. Petitet, ness of the plasmon description of the van der waals density functional,” K. Stanley, D. Walker, and R. Whaley, “Scalapack: a portable linear alge- Phys. Rev. B 89, 035412 (2014). bra library for distributed memory computers — design issues and perfor- 37 G. Román-Pérez and J. M. Soler, “Efficient implementation of a van der mance,” Comput. Phys. Commun. 97, 1 – 15 (1996). waals density functional: Application to double-wall carbon nanotubes,” 60 T. Auckenthaler, V. Blum, H.-J. Bungartz, T. Huckle, R. Johanni, Phys. Rev. Lett. 103, 096102 (2009). L. Krämer, B. Lang, H. Lederer, and P. Willems, “Parallel solution of 38 O. A. Vydrov and T. Van Voorhis, “Nonlocal van der waals density func- partial symmetric eigenvalue problems from electronic structure calcula- tional: The simpler the better,” J. Chem. Phys. 133, 244103 (2010). tions,” Parallel Comput. 37, 783 – 794 (2011). 39 L. Kong, G. Román-Pérez, J. M. Soler, and D. C. Langreth, “Energetics 61 A. Marek, V. Blum, R. Johanni, V. Havu, B. Lang, T. Auckenthaler, A. Hei- and dynamics of h2 adsorbed in a nanoporous material at low tempera- necke, H.-J. Bungartz, and H. Lederer, “The ELPA library: scalable paral- ture,” Phys. Rev. Lett. 103, 096103 (2009). lel eigenvalue solutions for electronic structure theory and computational 40 H. González-Herrero, J. M. Gómez-Rodríguez, P. Mallet, M. Moaied, J. J. science,” J. Phys.: Condens. Matter 26, 213201 (2014). Palacios, C. Salgado, M. M. Ugeda, J.-Y. Veuillen, F. Yndurain, and I. Bri- 62 P. Ku̇s, A. Marek, S. Koecher, H.-H. Kowalski, C. Carbogno, C. Scheurer, huega, “Atomic-scale control of graphene magnetism by using hydrogen K. Reuter, M. Scheffler, and H. Lederer, “Optimizations of the eigen- atoms,” Science 352, 437–441 (2016). solvers in the elpa library,” Parallel Comput. 85, 167 – 177 (2019). 41 J. Wang, G. Román-Pérez, J. M. Soler, E. Artacho, and M.-V. Fernández- 63 S. Goedecker, “Integral representation of the fermi distribution and its ap- Serra, “Density, structure, and dynamics of water: The effect of van der plications in electronic-structure calculations,” Phys. Rev. B 48, 17573– waals interactions,” J. Chem. Phys. 134, 024516 (2011). 17575 (1993). 42 J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Hybrid functionals based on 64 S. Mohr, W. Dawson, M. Wagner, D. Caliste, T. Nakajima, and L. Gen- a screened coulomb potential,” J. Chem. Phys. 118, 8207–8215 (2003). ovese, “Efficient Computation of Sparse Matrix Functions for Large-Scale 43 J. Heyd, G. E. Scuseria, and M. Ernzerhof, “Erratum: “hybrid function- Electronic Structure Calculations: The CheSS Library,” J. als based on a screened coulomb potential”,” J. Chem. Phys. 124, 219906 Chem. Theory Comput. 13, 4684–4698 (2017). (2006). 65 L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Wil- 44 A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, “Influ- land, D. Caliste, O. Zilberberg, M. Rayson, A. Bergman, and R. Schnei- ence of the exchange screening parameter on the performance of screened der, “Daubechies wavelets as a basis set for density functional pseudopo- hybrid functionals,” J. Chem. Phys. 125, 224106 (2006). tential calculations,” J. Chem. Phys. 129, 014109 (2008). 45 H. Shang, Z. Li, and J. Yang, “Implementation of screened hybrid density 66 L. Lin, M. Chen, C. Yang, and L. He, “Accelerating atomic orbital-based functional for periodic systems with numerical atomic orbitals: Basis func- electronic structure calculation via pole expansion and selected inversion,” tion fitting and integral screening,” J. Chem. Phys. 135, 034110 (2011). J. Phys.: Condens. Matter 25, 295501 (2013). 46 (), see: https://sourceforge.net/projects/libint/, provided by 67 L. Lin, A. García, G. Huhs, and C. Yang, “SIESTA-PEXSI: massively par- E. Valeev and J. T. Fermann. allel method for efficient and accurateab initiomaterials simulation without 47 S. Obara and A. Saika, “Efficient recursive computation of molecular in- matrix diagonalization,” J. Phys.: Condens. Matter 26, 305503 (2014). tegrals over cartesian gaussian functions,” J. Chem. Phys. 84, 3963–3974 68 W. Hu, L. Lin, C. Yang, and J. Yang, “Electronic structure and aromatic- (1986). ity of large-scale hexagonal graphene nanoflakes,” J. Chem. Phys. 141, 48 M. Head-Gordon and J. A. Pople, “A method for two-electron gaussian 214704 (2014). integral and integral derivative evaluation using recurrence relations,” J. 69 F. Corsetti, “The orbital minimization method for electronic structure cal- Chem. Phys. 89, 5777–5786 (1988). culations with finite-range atomic basis sets,” Comput. Phys. Commun. 49 J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approxi- 185, 873 – 883 (2014). mation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996). 70 T. Imamura, S. Yamada, and M. Machida, “Development of a high- 50 C. Kittel, Introduction to Solid State Physics (John Wiley & Sons, New performance eigensolver on a peta-scale next-generation supercomputer York, 1986). system,” Progress in Nuclear Science and Technology 2, 643–650 (2011). 51 S. H. Wemple, “Polarization fluctuations and the optical-absorption edge 71 J. Dongarra, M. Gates, A. Haidar, J. Kurzak, P. Luszczek, S. Tomov, and in batio3 ,” Phys. Rev. B 2, 2679–2689 (1970). I. Yamazaki, “Accelerating numerical dense linear algebra calculations 52 R. Cuadrado and J. I. Cerdá, “Fully relativistic pseudopotential formal- with GPUs,” Numerical Computations with GPUs , 1–26 (2014). ism under an atomic orbital basis: spin–orbit splittings and magnetic 72 V. Hernandez, J. E. Roman, and V. Vidal, “SLEPc: A scalable and flexible anisotropies,” J. Phys.: Condens. Matter 24, 086005 (2012). toolkit for the solution of eigenvalue problems,” ACM Transactions on 53 L. Fernández-Seivane, M. A. Oliveira, S. Sanvito, and J. Ferrer, “On- Mathematical Software 31, 351–362 (2005). site approximation for spin–orbit coupling in linear combination of atomic 73 W. Dawson and T. Nakajima, “Massively parallel sparse matrix function orbitals density functional methods,” J. Phys.: Condens. Matter 18, 7999– calculations with ntpoly,” Comput. Phys. Commun. 225, 154 – 165 (2018). 8013 (2006). 74 A. Tsolakidis, D. Sánchez-Portal, and R. M. Martin, “Calculation of the 54 F. Fernández-Seivane, M. A. Oliveira, S. Sanvito, and J. Ferrer, “Erratum: optical response of atomic clusters using time-dependent density func- On-site approximation for spin–orbit coupling in lcao density functional tional theory and local orbitals,” Phys. Rev. B 66, 235416 (2002). methods,” J. Phys.: Condens. Matter 19, 489001 (2007). 75 E. Artacho and D. D. O’Regan, “Quantum mechanics in an evolving S IESTA 27 hilbert space,” Phys. Rev. B 95, 115155 (2017). 100 K. W. Jacobsen, J. T. Falkenberg, N. Papior, P. Bøggild, A.-P. Jauho, 76 J. K. Tomfohr and O. F. Sankey, “Time-dependent simulation of conduc- and M. Brandbyge, “All-graphene edge contacts: Electrical resistance of tion through a molecule,” Phys. Status Solidi (b) 226, 115–123 (2001). graphene t-junctions,” Carbon 101, 101 – 106 (2016). 77 A. A. Correa, J. Kohanoff, E. Artacho, D. Sánchez-Portal, and A. Caro, 101 M. Kolmer, P. Brandimarte, J. Lis, R. Zuzak, S. Godlewski, H. Kawai, “Nonadiabatic forces in ion-solid interactions: The initial stages of radia- A. Garcia-Lekue, N. Lorente, T. Frederiksen, C. Joachim, D. Sanchez- tion damage,” Phys. Rev. Lett. 108, 213201 (2012). Portal, and M. Szymonski, “Electronic transport in planar atomic-scale 78 M. A. Zeb, J. Kohanoff, D. Sánchez-Portal, A. Arnau, J. I. Juaristi, and structures measured by two-probe scanning tunneling spectroscopy,” Nat. E. Artacho, “Electronic stopping power in gold: The role of d electrons Commun. 10, 1573 (2019). and the H/He anomaly,” Phys. Rev. Lett. 108, 225504 (2012). 102 P. Brandimarte, M. Engelund, N. Papior, A. Garcia-Lekue, T. Frederiksen, 79 R. Ullah, F. Corsetti, D. Sánchez-Portal, and E. Artacho, “Electronic stop- and D. Sánchez-Portal, “A tunable electronic beam splitter realized with ping power in a narrow band gap semiconductor from first principles,” crossed graphene nanoribbons,” J. Chem. Phys. 146, 092318 (2017). Phys. Rev. B 91, 125203 (2015). 103 N. Papior, “sisl: v0.9.7,” (2019). 80 J. Halliday and E. Artacho, “Anisotropy of electronic stopping power in 104 https://github.com/tfrederiksen/inelastica/. graphite,” Phys. Rev. B 100, 104112 (2019). 105 T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho, “Inelas- 81 F. Corsetti, “Performance analysis of electronic structure codes on hpc sys- tic transport theory from first principles: Methodology and application to tems: A case study of siesta,” PLOS ONE 9, 1–8 (2014). nanoscale devices,” Phys. Rev. B 75, 205413 (2007). 82 F. Corsetti, https://gitlab.e-cam2020.eu/esl/omm. 106 See: http://www.wannier.org. 83 E. Anderson, A. Benzoni, J. Dongarra, S. Moulton, S. Ostrouchov, 107 G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Giber- B. Tourancheau, and R. van de Geijn, “Basic linear algebra comrnuni- tini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. cation subprograms,” in The Sixth Distributed Memory Computing Con- Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. No- ference, 1991. Proceedings (1991) pp. 287–290. hara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, 84 L. S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel, I. Dhillon, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker, Mostofi, and J. R. Yates, “Wannier90 as a community code: new features and R. C. Whaley, ScaLAPACK Users’ Guide (Society for Industrial and and applications,” J. Phys.: Condens. Matter 32, 165902 (2020). Applied Mathematics, Philadelphia, PA, 1997). 108 N. Marzari and D. Vanderbilt, “Maximally localized generalized wannier 85 F. Tisseur and J. Dongarra, “A parallel divide and conquer algorithm for functions for composite energy bands,” Phys. Rev. B 56, 12847–12865 the symmetric eigenvalue problem on distributed memory architectures,” (1997). SIAM J. Sci. Comput. 20, 2223–2236 (1999). 109 N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, “Max- 86 T. N. Todorov, “Time-dependent tight binding,” J. Phys.: Condens. Matter imally localized wannier functions: Theory and applications,” Rev. Mod. 13, 10125–10148 (2001). Phys. 84, 1419–1475 (2012). 87 J. Halliday and E. Artacho, (in preparation) (2020). 110 D. Vanderbilt, Berry phases in electronic structure theory (Cambridge 88 I. Tavernelli, U. F. Röhrig, and U. Rothlisberger, “Molecular dynamics in University Press, 2018). electronically excited states using time-dependent density functional the- 111 https://github.com/romerogroup/DMFTwDFT. ory,” Mol. Phys. 103, 963–981 (2005). 112 V. Singh, U. Herath, B. Wah, X. Liao, A. H. Romero, and H. Park, 89 “Cpmd v3.13,” Copyright IBM Corp (1990-2008), copyright MPI fuer “Dmftwdft: An open-source code combining dynamical mean field Festkoerperforschung Stuttgart (1997-2001). theory with various density functional theory packages,” (2020), 90 P. López-Tarifa, M.-A. Hervé du Penhoat, R. Vuilleumier, M.-P. Gaigeot, arXiv:2002.00068v1. I. Tavernelli, A. Le Padellec, J.-P. Champeaux, M. Alcamí, P. Moretto- 113 X. Wu, A. Selloni, and R. Car, “Order-n implementation of exact exchange Capelle, F. Martín, and M.-F. Politis, “Ultrafast nonadiabatic fragmenta- in extended insulating systems,” Phys. Rev. B 79, 085102 (2009). tion dynamics of doubly charged uracil in a gas phase,” Phys. Rev. Lett. 114 P. García-Fernández, J. C. Wojdeł, J. Íñiguez, and J. Junquera, “Second- 107, 023202 (2011). principles method for materials simulations including electron and lattice 91 P. López-Tarifa, M.-A. Penhoat, R. Vuilleumier, M.-P. Gaigeot, U. Roth- degrees of freedom,” Phys. Rev. B 93, 195137 (2016). lisberger, I. Tavernelli, A. Le Padellec, J.-P. Champeaux, M. Alcami, 115 J. C. Wojdeł, P. Hermet, M. P. Ljungberg, P. Ghosez, and J. Íñiguez, “First- P. moretto capelle, F. Martín, and M. Politis, “Time-dependent density principles model potentials for lattice-dynamical studies: general method- functional theory molecular dynamics simulation of doubly charged uracil ology and example of application to ferroic perovskite oxides,” J. Phys.: in gas phase,” Cent. Eur. J. Phys. 12, 97–102 (2014). Condens. Matter 25, 305401 (2013). 92 M.-P. Gaigeot, P. Lopez-Tarifa, F. Martin, M. Alcami, R. Vuilleumier, 116 P. Torres, J. A. Seijas-Bellido, C. Escorihuela-Sayalero, J. Íñiguez, and I. Tavernelli, M.-A. Penhoat, and M. Politis, “Theoretical investigation R. Rurali, “Theoretical investigation of lattice thermal conductivity and of the ultrafast dissociation of ionised biomolecules immersed in water: electrophononic effects in srtio3 ,” Phys. Rev. Mater. 3, 044404 (2019). Direct and indirect effects,” Mutat. Res. 704, 45–53 (2010). 117 J. A. Seijas-Bellido, J. Íñiguez, and R. Rurali, “Anisotropy-driven ther- 93 A. D. Becke, “Density-functional exchange-energy approximation with mal conductivity switching and thermal hysteresis in a ferroelectric,” Appl. correct asymptotic behavior,” Phys. Rev. A 38, 3098–3100 (1988). Phys. Lett. 115, 192903 (2019). 94 C. Lee, W. Yang, and R. G. Parr, “Development of the colle-salvetti 118 J. A. Seijas-Bellido, H. Aramberri, J. Íñiguez, and R. Rurali, “Electric correlation-energy formula into a functional of the electron density,” Phys. control of the heat flux through electrophononic effects,” Phys. Rev. B 97, Rev. B 37, 785–789 (1988). 184306 (2018). 95 J. M. Pruneda, D. Sánchez-Portal, A. Arnau, J. I. Juaristi, and E. Artacho, 119 S. Das, Y. L. Tang, Z. Hong, M. A. P. Gonçalves, M. R. McCarter, “Electronic stopping power in lif from first principles,” Phys. Rev. Lett. C. Klewe, K. X. Nguyen, F. Gómez-Ortiz, P. Shafer, E. Arenholz, V. A. 99, 235501 (2007). Stoica, S.-L. Hsu, B. Wang, C. Ophus, J. F. Liu, C. T. Nelson, S. Saremi, 96 A. A. Correa, “Calculating electronic stopping power in materials from B. Prasad, A. B. Mei, D. G. Schlom, J. Íñiguez, P. García-Fernández, D. A. first principles,” Comp. Mater. Sci. 150, 291 – 303 (2018). Muller, L. Q. Chen, J. Junquera, L. W. Martin, and R. Ramesh, “Observa- 97 J. M. Pruneda, S. K. Estreicher, J. Junquera, J. Ferrer, and P. Ordejón, “Ab tion of room-temperature polar skyrmions,” Nature 568, 368–372 (2019). initio local vibrational modes of light impurities in silicon,” Phys. Rev. B 120 A. K. Yadav, K. X. Nguyen, Z. Hong, P. García-Fernández, P. Aguado- 65, 075210 (2002). Puente, C. T. Nelson, S. Das, B. Prasad, D. Kwon, S. Cheema, A. I. Khan, 98 N. Papior, N. Lorente, T. Frederiksen, A. García, and M. Brandbyge, “Im- C. Hu, J. Íñiguez, J. Junquera, L.-Q. Chen, D. A. Muller, R. Ramesh, and provements on non-equilibrium and transport green function techniques: S. Salahuddin, “Spatially resolved steady-state negative capacitance,” Na- The next-generation transiesta,” Comput. Phys. Commun. 212, 8 – 24 ture 565, 468–471 (2019). (2017). 121 R.Ierusalimschy, Programming in Lua, Fourth Edition (Feisty Duck Digi- 99 N. Papior, G. Calogero, S. Leitherer, and M. Brandbyge, “Removing all tal Book Distribution, 2016). periodic boundary conditions: Efficient nonequilibrium Green’s function 122 https://siesta-project.github.io/flos/ldoc/index.html (). calculations,” Phys. Rev. B 100, 195417 (2019). 123 S. Smidstrup, A. Pedersen, K. Stokbro, and H. Jónsson, “Improved initial S IESTA 28 guess for minimum energy path calculations,” J. Chem. Phys. 140, 214106 response of silver clusters and their hollow shells from linear-response (2014). TDDFT,” J. Phys.: Condens. Matter 28, 214001 (2016). 124 D. Sheppard, R. Terrell, and G. Henkelman, “Optimization methods for 148 Y. Li and C. A. Ullrich, “The particle–hole map: A computational tool to finding minimum energy paths,” J. Chem. Phys. 128, 134106 (2008). visualize electronic excitations,” J. Chem. Theory Comput. 11, 5838–5852 125 S. A. Trygubenko and D. J. Wales, “A doubly nudged elastic band method (2015). for finding transition states,” J. Chem. Phys. 120, 2082–2094 (2004). 149 R. Dronskowski and P. E. Blöchl, “Crystal orbital hamilton populations 126 G. R. Qian, X. Dong, X. F. Zhou, Y. Tian, A. R. Oganov, and H. T. Wang, (cohp): energy-resolved visualization of chemical bonding in solids based “Variable cell nudged elastic band method for studying solid-solid struc- on density-functional calculations,” The Journal of Physical Chemistry 97, tural phase transitions,” Comput. Phys. Commun. 184, 2111–2118 (2013). 8617–8624 (1993). 127 https://esl.cecam.org/Flook (). 150 M. Barbry, P. Koval, F. Marchesin, R. Esteban, A. G. Borisov, J. Aizpurua, 128 G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, “Ai- and D. Sánchez-Portal, “Atomistic Near-Field Nanoplasmonics: Reach- ida: automated interactive infrastructure and database for computational ing Atomic-Scale Resolution in Nanooptics,” Nano Lett. 15, 3410–3419 science,” Comput. Mater. Sci. 111, 218 – 230 (2016). (2015). 129 (), see: http://www.aiida.net/. 151 F. Marchesin, P. Koval, M. Barbry, J. Aizpurua, and D. Sánchez-Portal, 130 https://aiida.readthedocs.io/ (). “Optical response of metallic nanojunctions driven by single atom mo- 131 https://aiida-siesta-plugin.readthedocs.io/. tion,” ACS Photonics 3, 269–277 (2016). 132 S. G. Mayo, F. Yndurain, and J. M. Soler, “Band unfolding made simple,” 152 M. Barbry, Plasmons in Nanoparticles: Atomistic Ab Initio Theory for J. Phys.: Condens. Matter 32, 205902 (2020). Large Systems, Ph.D. thesis, University of Basque Country, Donostia-San 133 G. Calogero, N. Papior, M. Koleini, M. Larsen, and M. Brandbyge, Sebastián, Spain (2018), http://cfm.ehu.es/view/files/MArc_ “Multi-scale approach to first-principles electron transport beyond 100 barbry_2-1.pdf (accessed June 14 2019). nm,” Nanoscale 11, 6153–6164 (2019). 153 E. Lampin, P. L. Palla, P.-A. Francioso, and F. Cleri, “Thermal conduc- 134 A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, tivity from approach-to-equilibrium molecular dynamics,” J. Appl. Phys. M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Her- 114, 033525 (2013). mes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kols- 154 S. Illera, M. Pruneda, L. Colombo, and P. Ordejón, “Thermal and transport bjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, properties of pristine single-layer hexagonal boron nitride: A first princi- T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, ples investigation,” Phys. Rev. Mater. 1, 044006 (2017). M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, 155 S. García-Gil, A. García, and P. Ordejón, “Calculation of core level shifts and K. W. Jacobsen, “The atomic simulation environment—a python li- within dft using pseudopotentials and localized basis sets,” Eur. Phys. J. B brary for working with atoms,” J. Phys.: Condens. Matter 29, 273002 85, 239 (2012). (2017). 156 S. Garcia-Gil, A. Arnau, and A. Garcia-Lekue, “Exploring large o 1s and 135 A. Kokalj, “XCrySDen, a crystalline and molecular structure visualisation n 1s core level shifts due to intermolecular hydrogen bond formation in program,” http://www.xcrysden.org. organic molecules,” Surf. Sci. 613, 102 – 107 (2013). 136 K. Momma, “VESTA, a 3D visualization program for structural 157 See: https://en.wikipedia.org/wiki/Reference_counting. models, volumetric data, and crystal morphologies,” http://jp- 158 See: https://www.max-centre.eu/. minerals.org/vesta/en/. 159 See: https://esl.cecam.org/. 137 A. Postnikov and N. B. Mortazavi Amiri, “Calculated vibration spectrum 160 See: https://l_sim.gitlab.io/psolver/. of monoclinic Cu2 SnSe3 in comparison with kesterite-type Cu2 ZnSnSe4 ,” 161 P. Zubko, N. Jecklin, A. Torres-Pardo, P. Aguado-Puente, A. Gloter, Phys. Status Solidi A 210, 1332–1335 (2013). C. Lichtensteiger, J. Junquera, O. Stéphan, and J.-M. Triscone, “Elec- 138 M. N. Rao, D. Lamago, A. Ivanov, M. d’Astuto, A. V. Postnikov, trostatic coupling and local structural distortions at interfaces in ferroelec- R. Hajj Hussein, T. Basak, S. L. Chaplot, F. Firszt, W. Paszkowicz, S. K. tric/paraelectric superlattices,” Nano Lett. 12, 2846–2851 (2012). Deb, and O. Pagés, “Lattice dynamics of the model percolation-type 162 P. Aguado-Puente and J. Junquera, “Structural and energetic properties of (Zn,Be)Se alloy: Inelastic neutron scattering, ab initio study, and shell- domains in pbtio3 /srtio3 superlattices from first principles,” Phys. Rev. B model calculations,” Phys. Rev. B 89, 155201 (2014). 85, 184105 (2012). 139 “Bilbao Crystallographic Server → Raman and Hyper-Raman scattering 163 A. K. Yadav, C. T. Nelson, S. L. Hsu, Z. Hong, J. D. Clarkson, C. M. → Specral Active Modes,” https://www.cryst.ehu.es/rep/sam.html. Schlepütz, A. R. Damodaran, P. Shafer, E. Arenholz, L. R. Dedon, 140 V. Kosyak, N. B. Mortazavi Amiri, A. Postnikov, and M. A. Scarpulla, D. Chen, A. Vishwanath, A. M. Minor, L. Q. Chen, J. F. Scott, L. W. Mar- “Model of native point defect equilibrium in Cu2 ZnSnS4 and application tin, and R. Ramesh, “Observation of polar vortices in oxide superlattices,” to one-zone annealing,” J. Appl. Phys. 114, 124501 (2013). Nature 530, 198–201 (2016). 141 N. B. Mortazavi Amiri, Relation entre motifs structuraux et dynamique de 164 P. Shafer, P. García-Fernández, P. Aguado-Puente, A. R. Damodaran, A. K. réseau dans les cristaux mixtes Cu-Zn-Sn-Se: étude premiers principes, Yadav, C. T. Nelson, S.-L. Hsu, J. C. Wojdeł, J. Íñiguez, L. W. Martin, Ph.D. thesis, Université de Lorraine (2013). E. Arenholz, J. Junquera, and R. Ramesh, “Emergent chirality in the elec- 142 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford tric polarization texture of titanate superlattices,” Proc. Natl. Acad. Sci. University Press, Oxford, 1987). U.S.A. 115, 915–920 (2018). 143 A. V. Postnikov, “Vibrations in solids and small particles from first- 165 D. W. Boukhvalov, M. I. Katsnelson, and A. I. Lichtenstein, “Hydrogen principles calculations,” in Comput. Mater. Sci., NATO Science Series III on graphene: Electronic structure, total energy, structural distortions and (Computer and System Science), Vol. 187, edited by C. R. A. Catlow and magnetism from first-principles calculations,” Phys. Rev. B 77, 035427 E. Kotomin (IOS Press, 2003) pp. 153–166, proceedings of the NATO (2008). Advanced Study Institute “Computational Materials Science”, Il Ciocco, 166 O. V. Yazyev and L. Helm, “Defect-induced magnetism in graphene,” September 9–22, 2001. Phys. Rev. B 75, 125408 (2007). 144 O. Coulaud, P. Bordat, P. Fayon, V. Le Bris, I. Baraille, and R. Brown, 167 M. Slota, A. Keerthi, W. K. Myers, E. Tretyakov, M. Baumgarten, A. Ar- “Extensions of the Siesta DFT Code for Simulation of Molecules,” Re- davan, H. Sadeghi, C. J. Lambert, A. Narita, K. Müllen, and L. Bogani, search Report RR-8221 (INRIA, 2013). “Magnetic edge states and coherent manipulation of graphene nanorib- 145 P. Koval, M. Barbry, and D. Sánchez-Portal, “Pyscf-nao: An efficient and bons,” Nature 557, 691–695 (2018). flexible implementation of linear response time-dependent density func- 168 O. V. Yazyev and S. G. Louie, “Electronic transport in polycrystalline tional theory with numerical atomic orbitals,” Comput. Phys. Commun. graphene,” Nat. Mater. 9, 806–809 (2010). 236, 188 – 204 (2019). 169 O. V. Yazyev and S. G. Louie, “Topological defects in graphene: Disloca- 146 M. Petersilka, U. Gossmann, and E. Gross, “Excitation Energies from tions and grain boundaries,” Phys. Rev. B 81, 195420 (2010). Time-Dependent Density-Functional Theory,” Phys. Rev. Lett. 76, 1212– 170 W. Y. Kim and K. S. Kim, “Prediction of very large values of magnetore- 1215 (1996). sistance in a graphene nanoribbon device,” Nat. Nanotechnol. 3, 408–412 147 P. Koval, F. Marchesin, D. Foerster, and D. Sánchez-Portal, “Optical (2008). S IESTA 29 171 C. Moreno, M. Vilas-Varela, B. Kretz, A. Garcia-Lekue, M. V. Costache, Chem. B 107, 13728–13736 (2003). M. Paradinas, M. Panighel, G. Ceballos, S. O. Valenzuela, D. Peña, 184 C. F. Sanz-Navarro, R. Grima, A. García, E. A. Bea, A. Soba, J. M. and A. Mugarza, “Bottom-up synthesis of multifunctional nanoporous Cela, and P. Ordejón, “An efficient implementation of a qm-mm method graphene,” Science 360, 199–203 (2018). in siesta,” Theor. Chem. Acc. 128, 825–833 (2011). 172 W. Hu, Y. Huang, X. Qin, L. Lin, E. Kan, X. Li, C. Yang, and 185 G. T. Feliciano, A. J. R. da Silva, G. Reguera, and E. Artacho, “Molecular J. Yang, “Room-temperature magnetism and tunable energy gaps in edge- and electronic structure of the peptide subunit of geobacter sulfurreducens passivated zigzag graphene quantum dots,” NPJ 2D Mater. Appl. 3, 17 conductive pili from first principles,” J. Phys. Chem. A 116, 8023–8030 (2019). (2012). 173 S. Kim, J. Ihm, H. J. Choi, and Y.-W. Son, “Origin of anomalous electronic 186 D. J. Cole and N. D. M. Hine, “Applications of large-scale density func- structures of epitaxial graphene on silicon carbide,” Phys. Rev. Lett. 100, tional theory in biology,” J. Phys.: Condens. Matter 28, 393001 (2016). 176802 (2008). 187 M. Darvish Ganji, “Amino acids interacting with defected carbon nan- 174 H. S. S. Ramakrishna Matte, A. Gomathi, A. K. Manna, D. J. Late, otubes: ab initio calculations,” Journal of Pharmaceutical & Health Sci- R. Datta, S. K. Pati, and C. N. R. Rao, “Mos2 and ws2 analogues of ences 4, 157–166 (2016). graphene,” Angew. Chem. Int. Ed. 49, 4059–4062 (2010). 188 W. Li, K. Kotsis, and S. Manzhos, “Comparative density functional theory 175 I. Popov, G. Seifert, and D. Tománek, “Designing electrical contacts to and density functional tight binding study of arginine and arginine-rich cell mos2 monolayers: A computational study,” Phys. Rev. Lett. 108, 156802 penetrating peptide tat adsorption on anatase tio2,” Phys. Chem. Chem. (2012). Phys. 18, 19902–19917 (2016). 176 H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, 189 A. Hermann, “Chemical bonding at high pressure,” in Reviews in Com- “Phosphorene: An unexplored 2d semiconductor with a high hole mobil- putational Chemistry (John Wiley & Sons, Ltd, 2017) Chap. 1, pp. 1–41, ity,” ACS Nano 8, 4033–4041 (2014). https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119356059.ch1. 177 J. Guan, Z. Zhu, and D. Tománek, “Phase coexistence and metal-insulator 190 T. Liu, E. Artacho, F. Gázquez, G. Walters, and D. Hodell, “Prediction of transition in few-layer phosphorene: A computational study,” Phys. Rev. equilibrium isotopic fractionation of the gypsum/bassanite/water system Lett. 113, 046804 (2014). using first-principles calculations,” Geochim. Cosmochim. Ac. 244, 1 – 11 178 W. Hu, L. Lin, C. Yang, J. Dai, and J. Yang, “Edge-modified phosphorene (2019). nanoflake heterojunctions as highly efficient solar cells,” Nano Lett. 16, 191 E. Mayoral, A. Rey, J. Klapp, A. Gómez, and M. Mayoral, “Ab initio 1675–1682 (2016). dft calculations for materials in nuclear research,” in High Performance 179 I. B. Guster, A bird’s-eye view of charge and spin density wave from first Computing, edited by C. J. Barrios Hernández, I. Gitler, and J. Klapp principles calculations, Ph.D. thesis, Universitat Autonoma de Barcelona (Springer International Publishing, Cham, 2017) pp. 329–339. (2019). 192 R. Escribano and G. M. Muñoz Caro, “Introduction to spectroscopy and 180 B. Guster, C. Rubio-Verdú, R. Robles, J. Zaldívar, P. Dreher, M. Pruneda, astronomical observations,” in Laboratory Astrophysics, edited by G. M. J. A. Silva-Guillén, D.-J. Choi, J. I. Pascual, M. M. Ugeda, P. Ordejón, Muñoz Caro and R. Escribano (Springer International Publishing, Cham, and E. Canadell, “Coexistence of elastic modulations in the charge density 2018) pp. 27–47. wave state of 2h-nbse2,” Nano Lett. 19, 3027–3032 (2019). 193 F. Viñes, J. R. B. Gomes, and F. Illas, “Understanding the reactivity of 181 B. Guster, M. Pruneda, P. Ordejón, E. Canadell, and J.-P. Pouget, “Evi- metallic nanoparticles: beyond the extended surface model for catalysis,” dence for the weak coupling scenario of the peierls transition in the blue Chem. Soc. Rev. 43, 4922–4939 (2014). bronze,” Phys. Rev. Mater. 3, 055001 (2019). 194 F. Tao, W. Schneider, and P. Kamat, Heterogeneous Catalysis at 182 P. J. de Pablo, F. Moreno-Herrero, J. Colchero, J. Gómez Herrero, P. Her- Nanoscale for Energy Applications (Wiley, 2015). rero, A. M. Baró, P. Ordejón, J. M. Soler, and E. Artacho, “Absence of 195 https://www.cp2k.org/dbcsr. dc-conductivity in λ -dna,” Phys. Rev. Lett. 85, 4992–4995 (2000). 196 I. Sivkov, P. Seewald, A. Lazzaro, and J. Hutter, “DBCSR: A blocked 183 A. Crespo, D. A. Scherlis, M. A. Martí, P. Ordejón, A. E. Roitberg, and sparse tensor algebra library,” (2019), arXiv:1910.13555. D. A. Estrin, “A dft-based qm-mm approach designed for the treatment 197 http://esl.cecam.org/libOMM (). of large molecular systems: application to chorismate mutase,” J. Phys. 198 A. Marcolongo, P. Umari, and S. Baroni, “Microscopic theory and quan- tum simulation of atomic heat transport,” Nat. Phys. 12, 80–84 (2016).