Siesta: Recent developments and applications

SIESTA 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 SIESTA program is presented. Since its debut in the mid-nineties, SIESTA’s flexibility, efficiency and free distribution has given advanced materials simulation capabilities to many groups worldwide. The core methodological scheme of SIESTA combines finite-support pseudoatomic 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 implementations on top of that core scheme, which include: full spin-orbit interaction, non-repeated and multiple-contact ballistic electron transport, DFT+Uand hybrid functionals, time-dependent DFT, novel reduced-scaling solvers, densityfunctional 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 WANNIER90 and the second-principles modelling it can be used for, an AiiDA plugin for workflow automatization, interface to Lua for steering SIESTA runs, and various postprocessing utilities. SIESTA has also been arXiv:2006.01270v1 [physics.comp-ph] 1 Jun 2020 “This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in García, A et al., J. Chem. Phys. 152, 204108 (2020) and may be found at: https://doi.org/10.1063/5.0005077" SIESTA 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. The possibility of treating large systems with firstprinciples electronic-structure methods has opened up new research avenues in many disciplines. The SIESTA method and its implementation have been key in this development, offering an efficient and flexible simulation paradigm based on the use of strictly localized basis sets. This approach enables the implementation of reduced scaling algorithms, and its accuracy and cost can be tuned in a wide range, from quick exploratory calculations to highly accurate simulations matching the quality of other approaches, such as plane-wave methods. The SIESTA method has been described in detail in Ref. 1, with an update in Ref. 2. In this paper we shall describe its present status, highlighting its strengths and documenting the steps that have recently been taken to improve its capabilities, performance, ease of use, and visibility in the electronicstructure community. a)Electronic mail: [email protected] b)Electronic mail: [email protected] c)Electronic mail: [email protected] d)Electronic mail: [email protected] e)Electronic mail: volker[email protected] f)Electronic mail: [email protected] g)Electronic mail: [email protected] h)Electronic mail: [email protected] i)Electronic mail: [email protected] j)Electronic mail: [email protected] k)Electronic mail: [email protected] l)Electronic mail: [email protected] m)Electronic mail: [email protected] n)Electronic mail: [email protected] o)Electronic mail: [email protected] p)Electronic mail: [email protected] q)Electronic mail: [email protected] r)Electronic mail: [email protected] s)Electronic mail: [email protected] t)Electronic mail: [email protected] u)Electronic mail: ko[email protected] v)Electronic mail: i.lebede[email protected] w)Electronic mail: [email protected]y.edu x)Electronic mail: [email protected] y)Electronic mail: [email protected] z)Electronic mail: [email protected] aa)Electronic mail: [email protected] bb)Electronic mail: andrei.postniko[email protected] cc)Electronic mail: [email protected] dd)Electronic mail: [email protected] ee)Electronic mail: [email protected] ff)Electronic mail: [email protected] gg)Electronic mail: [email protected] hh)Electronic mail: [email protected] ii)Electronic mail: [email protected] jj)Electronic mail: javier[email protected] As we shall see, the improvements touch many areas. We can underline the implementation of new core electronicstructure features (DFT+U, spin-orbit interaction, hybrid functionals), modes of operation (improved time-dependent density functional theory (TD-DFT), density functional perturbation theory (DFPT), and analysis methods and procedures to access new properties. A major effort has been spent in enhancing the interoperability of the code at various levels (sharing of pseudopotentials, a new wannierization interface opening the way to sophisticated post-processing, and an interface to multiscale methods). Very significant performance enhancements have been made, notably to the TRANSIESTA module through improved algorithms, and to the core electronic structure problem through the development of interfaces to new solvers. These advances have put SIESTA in a prominent place in the high-performance electronic-structure simulation scene, a role reinforced by its participation in important international initiatives and by its new open-source licensing model. The manuscript is organized as follows. We provide an overview of the underlying methodology and the capabilities of SIESTA in section II, which serves to place the code in the wider ecosystem of electronic-structure materials simulation. Section III presents the recent developments in and around the code, which are covered in sub-sections. To showcase SIESTA’s utility in the context of electronic-structure calculations, we present briefly some relevant applications and survey a few areas in which SIESTA is being profitably used in section IV. Plans for the future evolution of SIESTA are outlined in section V. II. KEY CONCEPTS OF SIESTA A. Theory background and context SIESTA appeared as a consequence of the push for linearscaling electronic structure methods of the mid nineties, which has been reviewed, for example, in Refs. 3and 4. SIESTA was the first linear-scaling self-consistent implementation of density functional theory (DFT).5,6 The SIESTA method relies on atomic-like functions of finite support as basis sets7,8– of arbitrary number, angular momentum, radial shape, and centers – combined with a discretization of space for the computation of the Kohn-Sham Hamiltonian terms that involve more than two centers. The electron-ion interaction is represented by norm-conserving pseudopotentials. These key ingredients, through the optimized handling of sparse matrices, are used to compute the self-consistent Hamiltonian and overlap matrices with a computational expense that scales linearly with system size. The SIESTA 3 method is completed with a choice of solvers for that Hamiltonian, from optimized (but cube-scaling) diagonalization methods, to reduced-scaling solvers of different flavors. The orbitals in the SIESTA basis set are made of the product of a real spherical harmonic and a radial function, which is numerically tabulated in a grid. The shape of the radial part is in principle totally arbitrary, but the experience accumulated has proven that the numerical solution of the Schrodinger equation for a (confined) isolated atom with the corresponding pseudopotential is a very good choice in terms of accuracy versus computational cost. Fuller descriptions of the mechanisms to generate and optimize these pseudo-atomic orbitals (PAOs) are given in Refs. 8–10. The auxiliary real-space grid is an essential ingredient of the method, as it allows the efficient representation of charge densities and potentials, as well as the computation of the matrix elements of the Hamiltonian that cannot be handled as two-center integrals. This grid can be seen as the reciprocal space of a set of plane waves, and its fineness is most conveniently parametrized by an energy cutoff (the “density” cutoff of plane-wave methods). There are limits to the softness of the functions that can be described with such a grid, so core electrons are not considered (although semi-core electrons usually are), and their effect is incorporated into pseudopotentials. The real-space grid is also used to solve the Poisson equation involved in the computation of the electrostatic potential from the charge density, through the use of a fast-fourier-transform method. This means that SIESTA uses periodic boundary conditions. Non periodic systems, such as molecules, tubes, or slabs, are treated using appropriate supercells. SIESTA is now a mature code with more than 20 years of existence. In this period, the most important algorithms behind our implementation have been already fully described and documented in a series of papers. Readers interested in the details of how the basic elements defining the method are combined, as well as other relevant implementation details that make the method practical, can find them in the main SIESTA reference1, and in the update with the new capabilities of the code2. We note that the term SIESTA is regularly used to describe both the method (as outlined in the earliest papers5,6) and its implementation in a computer program. The SIESTA method is at the basis of later independent implementations, such as OpenMX,11 and QuantumATK.12 Other subsequent codes built on the method, revising some of the fundamental ingredients. This is the case of FHI-aims,13 which uses a more sophisticated real-space grid (atom-centered), thus extending the core scheme to all-electron calculations. In this paper we describe new additions to the SIESTA code, based on independent methodological advances, either preexistent or specifically developed for SIESTA, as specified and cited in each section below. B. Overview of Siesta capabilities As a general purpose implementation, SIESTA can provide the standard functionality available in mainstream DFT codes: energies, forces, molecular-dynamics simulations, band structures, densities of states, etc., and shares with those codes the basic current limitations of DFT (notably the description of strongly-correlated systems). What makes SIESTA different from most other codes, and is at the root of its key strengths, is the atomic-like, and strictly localized, character of its basis set. The use of a “good first approximation” to the full problem implies, first, that a much smaller number of basis functions is needed. Second, the finite-support of the orbitals leads to sparsity and the possibility to use reduced-scaling methods. Thus high performance emerges almost by default. Take first the basis cardinality: the number of basis orbital per atom in a typical SIESTA calculation is of the order of 10-20. This is to be compared with a few hundred in the typical plane-wave (PW) calculation. Furthermore, for systems whose description needs a vacuum region (e.g., slabs for surface calculations, 2D monolayers, etc), empty space is essentially “free” for SIESTA, whereas PW codes still need a basis set determined by the total size of the simulation cell. SIESTA is then quite capable of dealing with systems composed of dozens to hundreds of atoms on modest hardware, even when using cubic-scaling diagonalization solvers, which are the default as they are universally applicable. Electronic-structure solvers with a more favorable sizescaling can be applied to suitable systems. For example, one of SIESTA’s earlier calculations, in 1996, was a linear-scaling run for a strand of DNA with 650 atoms, performed on a desktop workstation of the era.6Reduced size-scaling is also a feature of the PEXSI solver described in section IIIG1 below, and of the NTPoly solver mentioned in section IIIG2. In addition to time-to-solution efficiency, these solvers have a smaller memory footprint than diagonalization, as the relevant matrices are kept in sparse form rather than converted to a dense format. Crucially, SIESTA’s baseline efficiency can be scaled up to ever-larger systems by parallelization. Both distributed (MPI) and shared-memory (OpenMP) parallelization options are implemented in the code. As some of the examples in section IV show, non-trivial calculations with thousands of atoms are used in applications in different contexts, from molecular biology to electronic transport. Work on the performance aspects of the code is continuous, mostly on the solvers, which usually take most of the computer time due to the very high efficiency of the Hamiltonian setup module in SIESTA. This task is facilitated (see section IIIO) by leveraging external libraries and developments generated by a number of international initiatives in which SIESTA participates. The code can still run efficiently in modest hardware, while also being able to exploit massive levels of parallelism in large supercomputers (see Fig. 1). 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 language of atomic orbitals. Hence SIESTA has a natural advantage in this area. Partial densities of states and atomic and crystal populations (COOP/COHP) are routinely used to gain insights into the stability and other properties of materials. SIESTA 4 9 B. Accuracy of the SIESTA-PEXSI approach For one insulating and one metallic system (The size of such system does not need to be very large), show the accuracy of the PEXSI approach compared to the result from diagonalization with increasing number of poles, after the SCF iteration. C. Efficiency of the SIESTA-PEXSI approach The sets of DNA and C-BN examples are the base for examining the growth of the computational cost with system size (weak scaling) as well as the increase of the solution time with the number of processors (strong scaling). The analysis is based on the time for the calculating the first SCF step, including the setup of the Hamiltonian and, in case of diagonalization, computation of the density matrix based on the results of the ScaLAPACK eigenvalue solver. PEXSI uses 40 poles, which requires for all systems two inertia counts and one µiteration. In subsequent SCF iterations information about the chemical potential can be used for lowering the number of inertia counts or even completely omitting it, reducing the time per iteration even further. Siesta-PEXSI is particularly suitable for high performance computing, since the two levels of parallelization allow using a large number processors efficiently. The total number of processes can be varied by tuning the number of processes per pole (ppp) and the number of poles treated in parallel. The e↵ect of both is demonstrated in figure 2 for the largest DNA and C-BN systems examined. Configurations using the same ppp are connected with lines and show very good scaling. The first point on each line represents no parallelization over poles, while the last point corresponds to full parallelization. The inefficiencies in this regime mainly come from the symbolic factorization. This part can use only a limited number of processors smaller than ppp and thus does not scale at all, a↵ecting the performance notable. This is only a technical issue, related to the libraries currently used, and will be resolved in future. Then the time for symbolic factorization will play a only a marginal role. Increasing the number of processors per pole, demonstrated by points with the same number of poles treated in parallel, allows reducing the time even further, but scales less efficient than the pole-parallelization. [LL: Why is this the case? For C-BN it seems that the scaling from 144ppp to 400ppp reduces the wall clock time by a factor of 2, which is reasonably good.] Due to the similar numbers of orbitals of both examples, diagonalization times are alike, but throughout the tests much higher than the sulution times of SiestaPEXSI. In the case of C-BN0.00 the Siesta-PEXSI approach is one order of magnitude faster and allows an efficient use of more than 10000 cores, while the scalig of diagonalization is limited to about half of this. For DNA-25 less processors per pole are used since this example features sparser matrices. On the other hand this sparsity makes the solver work two orders of magnitude faster than diagonalization. Another consequence of PEXSI dealing only with sparse matrices is the smaller demand of memory. While on Edison the memory of at least 1000 cores is needed 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 than diagonalization with 5120 processors. FIG. 2. Strong scaling of C-BN0.00 and DNA-25 based on the total time for the first SCF step. The various lines for PEXSI result from using di↵erent numbers of processors per pole (ppp), while the points on each curve belong to computations with 1, 2, 5, 10, 20, and 40 poles in parallel. PEXSI’s beneficial scaling with the system size, as described in section II C, guarantees that for large enough systems Siesta-PEXSI will always be faster than diagonalization. The scaling of the computational cost is demonstrated for DNA and C-BN in figure 3. In all tests full parallelization over poles is used. In this configuration the influence of the symbolic factorization would change the character of the method. Because in future this influence will be negligible, the time for symbolic factorization is not taken into account for the analysis. The numbers of processes for each system size are chosen to be an efficient trade-o↵of reducing the time to solution while keeping the cost, which increases with the number of processes due to inefficiencies, as small as possible. Following this guideline it turns out, that for CBN one can use more processors with Siesta-PEXSI than with ScaLAPACK. This also means, that the advantage of Siesta-PEXSI in terms of solution-time is even larger than the benefit of cost. For very sparse problems, like the largest DNA examples, the amount of processors that can be used is similar for both methods, but SiestaPEXSI is about two orders of magnitude faster. More details are listed in table II. The analysis shows, besides Siesta-PEXSI’s favorable Strong scaling 170,000 orbs 180,000 orbs 1D, sp=0.27% 2D, sp=0.91% FIG. 1. Parallel strong scaling of SIESTA-PEXSI and the (Scalapack) diagonalization approach for a DNA chain and a GrapheneBoron Nitride stack, prototypes of large (hundreds of thousands of orbitals) quasi one-dimensional and two-dimensional systems. “ppp” stands for the number of MPI processes used in each pole computation, and “sp” the sparsity of the Hamiltonian. (For more details, see section IIIG1) For a recent example, see Ref. 14. Similarly, an atomic basis provides a very natural and adequate language for the firstprinciples simulation of electronic ballistic transport in nanosized systems, via the Green’s-function based Keldysh formalism implemented in TRANSIESTA,15 a part of the SIESTA package. The very high number of citations of the SIESTA papers testify to the successful application of the code to widely different systems. With regard to specific capabilities and the levels of accuracy achievable, we can distinguish several levels. First, SIESTA implements DFT, one of the most versatile materials simulation frameworks. DFT has its shortcomings, notably in regard to the description of strongly-correlated systems, but these are being addressed (see sections on DFT+U and hybrid functionals below). Second, SIESTA uses pseudopotentials to represent the electron-ion interaction. The pseudopotential approach is firmly rooted in a sound physical approximation (that bonding effects depend mostly on the valence electrons); however, it is at a disadvantage when coreelectrons effects are important (but see section III N6 below). Third, SIESTA employs periodic boundary conditions (PBC) for the solution of the Poisson problem, sharing with planewave codes the need to resort to repeated supercells for the study of low-dimensional systems, and to special techniques for the treatment of charged systems. It is important to note, however, that, unlike plane-wave codes, SIESTA is only bound to PBC because of the present treatment of the Hartree term of the single-particle Hamiltonian. This limitation is lifted by the incorporation of alternative Poisson solvers, as described in Sec. IIIO, which allow for open boundary conditions, as for isolated nano-systems, and hybrid open/periodic boundary conditions in different dimensions, as for isolated wires and slabs. It should be remembered that the three approximations mentioned in this paragraph are very widely used in the community, shared by some of the most popular electronic structure codes. Fourth, with regard to SIESTA specific approximations, particularly the basis set, it should be stressed that SIESTA is limited to basis sets composed of functions that are product of a radial part and spherical harmonics, but it does not constrain on how many, where such functions are centered, and the size of their finite-support region. Calculations can flexibly range from quick exploration to very high-quality simulations (one may recall that accuracy gold standards in electronic structure are provided by quantum-chemistry methods, based on LCAO). The use of an atomic-orbital basis set implies however the limitation of non-uniformity of convergence. As opposed to plane-wave methods, in which a single energy cutoff parameter monotonically determines the quality of the calculation, there is no univocal procedure for the choice of an appropriate basis set. It is a well-known problem, shared by the whole quantum-chemistry community, on which there is widely used and tested know-how. As Fig. 2shows, it is possible to attain in practice an accuracy comparable to that of well-converged plane-wave calculations. The reader is also referred to sections IVB and IV C1 for showcase examples of the accuracy of the code, among many others in the literature. To close this section, we stress that it has been a traditional and deliberate attitude by the SIESTA team that, although proposing sensible starting points to users as defaults, the choice of fundamental approximations and inputs to the program (not only basis sets, but also density functionals, and pseudopotentials) is a responsibility of the users, who retain full control and the flexibility to adapt the code to their specific needs. Nevertheless, tools for basis optimization are provided with the program, new curated databases of pseudopotentials are coming online, and new ways to ameliorate the correlation problem are being implemented. Some of these developments are described in the following sections. III. RECENT DEVELOPMENTS IN SIESTA A. New distribution model and development infrastructure A few years ago, in 2016, a decision was made to change the licensing model for SIESTA: traditionally it had always been free of charge to academics, but non-academic use required a special license and redistribution was not permitted. Now SIESTA is formally an open-source program, distributed according to the terms of the GPL license18. At the same time, the development infrastructure was made more transparent and scalable, using first the Launchpad platform19 and now the Gitlab service20. The net effect of the changes has been a more fertile and dynamic development, with more contributors who can have direct access to the various branches of development, and a better experience for users, who can download code and raise issues in an integrated platform. These changes have been substantial for the core developers, and the transitory period is still being felt. The main code base is gradually absorbing new developments, both those that SIESTA 5 230 240 250 30 50 70 90 PW Eb (meV) Number of basis functions -0.8 -0.4 0.0 70 75 80 85 90 Eb - Eb PW (meV) Number of basis functions FIG. 2. Basis set convergence for the binding energy (Eb) of a water dimer. Details can be found in Ref. 16. The horizontal dotted line represents the converged plane-wave (PW) calculation (1300 eV cutoff) for the same system (dimer geometry and box), pseudopotentials and density functional, using the ABINIT code.17 Inset: deviation of Ebversus the PW reference. The deviation for the last point is of 10 µeV. were planned long in advance, and new ones made possible by the greater openness and fluidity of the development model. Most of the new features described below are already part of public releases, but a few are undergoing the last stages of testing before release. The work-flow is also moving from long-lived releases, hard to maintain with bug-fixes, to more frequent releases that will be maintained for a shorter time. B. New pseudopotential format for interoperability PSML (for PSeudopotential Markup Language)21,22 is a file format for norm-conserving pseudopotential data which is designed to encapsulate as much as possible the abstract concepts involved in the domain, and to provide appropriate metadata and provenance information. This extra level of formalization aims at removing the interoperability problems associated to bespoke pseudopotential formats, which usually were designed to serve the needs of specific generators and client codes, and thus contain implicit assumptions about the meaning of the data or lack information not considered relevant. PSML files can be produced by the ONCVPSP23 and ATOM24 pseudopotential generator programs, and are a download-format option in the Pseudo-Dojo database of curated pseudopotentials25,26. The software library libPSML21,22 can be used by electronic structure codes to transparently extract the information in a PSML file and incorporate it into their own data structures, or to create converters for other formats. It is currently used by SIESTA and ABINIT,17,27 making possible a full pseudopotential interoperability and facilitating comparisons of calculation results. The use of this new format opens the door to benefit from the availability of a periodic table of reliable and accurate norm-conserving pseudopotentials, easing in most cases the task of pseudopotential quality control. C. DFT+U for correlated systems The LDA+U method, initially developed by Anisimov and coworkers28 with the objective to improve the treatment of the electron-electron interaction for localized electrons within the bare LDA description, has been implemented in SIESTA. The idea behind the LDA+U consists in describing the “strongly correlated” electronic states of a system (typically, localized dor forbitals) using the Hubbard model, whereas the rest of valence electrons are treated at the level of “standard” approximate DFT functionals.29 In the current version of SIESTA the implementation is based on the simplified rotationally invariant functional proposed by Dudarev and coworkers.30 Here, the corrections are made invariant under rotation of the atomic orbitals used to define the occupation number of the correlated subspace, at the cost of retaining only the lowest order Slater integrals in the factorization of the integrals of the Coulomb 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 occupation number nIσ `mof the localized correlated orbital `m with spin σwithin the atom Iis given by EU=∑ Iσ` UI` 2∑ m nIσ `m1−nIσ `m,(1) where only one interaction parameter UI`is needed to specify the interaction per atom and `-shell. In the practical SIESTA implementation, the populations on the correlated orbitals are computed using non-overlapping (i. e. orthogonal) localized projectors. They can be generated using either (i)the same algorithm used to produce the first-ζorbitals of the basis set, but with a larger energy shift, or (ii)cutting the exact solution of the pseudoatom with a Fermi function. The results of the LDA+U method are sensitively dependent on the numerical value of the effective on-site electronic interaction, the Hubbard U. Although in principle the value of Ucan be computed from first principles using linear response methods,31 a common practice is to tune it semiempirically, seeking agreement of certain properties (for instance band gaps or magnetic moments) with available experimental measurements. Then, the fitted Uis used in subsequent calculations to predict other properties. The LDA+U corrects localized states, for which the selfinteraction correction is expected to be stronger, and is an effective method to improve the description of the (underestimated) band gap of insulators, as shown in Fig. 3for the case of NiO. Once the Hubbard correction is switched on, the optical band gap increases up to 3.08 eV (from the bare GGAPBE value of 1.08 eV), very close to the experimental value for the onset of optical absorption in NiO32 (3.10 eV). The magnetic moment on the Ni atom is also properly described, with a value of 1.67 µBwhich lies well within the experimental range of values (between 1.64 µB33 and 1.9 µB34), and SIESTA 6 improves on the result of 1.39 µBobtained with a bare GGAPBE functional. -10 -5 0 5 Energy (eV) ΓL K T ΓX (a) -10 -5 0 5 Energy (eV) ΓL K T ΓX (b) FIG. 3. Band structure of NiO in the undistorted rock-salt type structure with rhombohedral symmetry introduced by a type-II antiferromagnetic order. The experimental lattice spacing is used. The bands obtained within GGA-Perdew-Burke-Ernzerhof functional (panel a), and with a Hubbard U correction of 4.6 eV applied on the d-orbitals of Ni (panel b), as in Ref. 31, are shown. The zero of the energy is set at the top of the valence band. D. Van der Waals functionals An efficient calculation of van der Waals (vdW) functionals35,36 was developed and first implemented in SIESTA using a polynomial expansion in the local variables (q1,q2)of the nonlocal interaction kernel Φ(q1,q2,r12)and a Fourier expansion in the relative position r1237. As a result, the scaling of the vdW computation decreases from O(N2)to O(NlogN)and it becomes marginal within the overall cost. This scheme was later extended16 to a more complex kernel38 of the form Φ(n1,|∇n1|,n2,|∇n2|,r12), and it has been applied to a large variety of systems, like carbon nanotubes37, hydrogen adsorption39,40, or liquid water41. E. Hybrid functionals The screened hybrid functional HSE0642–44 has been implemented in SIESTA 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 screened, so it has a shorter range than 1/r. Here, to reduce the big prefactor involved in the computation of the HFX potential 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 and efficient analytical way. An example of this fitting for the 2sand 2patomic orbitals basis set of the oxygen is shown in Fig. 4. The LIBINT package46 is required to calculate primitive ERIs, where recursive schemes of the Obara-Saika47 method and the Head-Gordon and Pople’s variation48 thereof 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. This HSE06 functional has been used to compute the band 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-BurkeErnzerhof functional49 within the Generalized Gradient Approximation (5.499 Å for Si and 4.033 Å for BaTiO3). In both cases, the gap is opened with respect to the value obtained with the semilocal functional. In bulk Si the band gap is indirect: the top of the valence band is located at Γand the bottom of the conduction band at a point along the Γ→X high-symmetry line. It increases from 0.64 eV within GGA to 1.00 eV with the hybrid functional, in good agreement with the experimental value of 1.17 eV50. For the case of the perovskite oxide BaTiO3, the band gap is also indirect, from Rto Γ, and its value increases from 1.87 eV with GGA to 3.28 eV with the HSE06 functional, almost matching the experimental 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 SIESTA and in the analysis tools is seen as a strategic asset for the project in view of the recent interest in topological insulators and quasi–two–dimensional systems with important spin–orbit effects, like some of the transition metal dichalcogenides. Also, it brings the possibility to obtain the magnetic crystalline anisotropy (MCA) (change in the total energy of the system upon changing the spin quantization axis). In a standard collinear-spin DFT calculation, the total KS Hamiltonian is represented by two independent spin–blocks, ˆ Hσσ µν [σ=↑,↓]. However, when the SO coupling is included, off–diagonal spin blocks arise (i.e., there are non–zero couplings between the two spin components). Therefore, and similar to the non-collinear spin case, the Hamiltonian be- SIESTA 7 0.5 1.0 1.5 2.0 2.5 r (Bohr) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Rnl ( r ) rl 0.5 1.0 1.5 2.0 2.5 r (Bohr) 0.0 1.0 2.0 3.0 4.0 5.0 Rnl ( r ) rl FIG. 4. Gaussian fits of the radial part of oxygen 2s(a) and 2p(b) orbitals using 6 Gaussian functions. The orbitals to fit are represented by blue dots and the corresponding Gaussian expansions by green continuous lines. Dashed vertical lines represent the standard deviations of individual Gaussians and a red continuous line marks their upper limit. The orbitals are set to zero in the yellow area, marking their cutoff radii. comes a full 2×2 matrix in spin space ˆ HKS µν = ˆ H↑↑ µν ˆ H↑↓ µν ˆ H↓↑ µν ˆ H↓↓ µν !(2) where µν subindexes refer to the SIESTA basis orbitals. The fully relativistic Hamiltonian ˆ HKS is expressed as a sum of the kinetic energy ˆ T, the scalar-relativistic pseudo-potential part in the form of Kleinman–Bylander projectors ˆ VKB, the spin– orbit ˆ VSO term and the Hartree ˆ VHand exchange–correlation ˆ VXC potentials: ˆ HKS =ˆ T+ˆ VKB +ˆ VSO +ˆ VH+ˆ VXC (3) The first three terms of the right hand side do not depend on the charge density, ρ(r), and therefore do not change in the self–consistent cycle, while ˆ VSO and ˆ VXC are the only spin– dependent terms that couple both spin components. In order to compute the MCAs, different orientations of the spin quantization axis need to be considered. This may be done by rotating either ˆ VSO (as done by Cuadrado and Cerdá52) or the density matrix, which is the approach currently followed by SIESTA for compatibility with the non– collinear case. -15 -10 -5 0 5 10 Energy (eV) LΓXWKΓ (a) -6 -4 -2 0 2 4 6 8 Energy (eV) ΓX M R ΓM X R (b) FIG. 5. Band structure of (a) bulk Si in the diamond structure, and (b) bulk BaTiO3in the cubic structure obtained with the PerdewBurke-Ernzerhof functional (red lines) and with the HSE06 hybrid functional (black lines). The zero of energies have been set to the valence band maximum. In the current implementation the SO term is included non-perturbatively, so that the fully relativistic Hamiltonian is solved self-consistently after extending the Kohn–Sham wave–functions to full spinors. Two different approaches have been implemented in SIESTA to account for the SO term, ˆ VSO: -on–site approximation: Based on the work of Fernández-Seivane et al. 53,54, only the intra–atomic SO contributions within each l– shell of each atom are considered. In this approach the SO terms are obtained from analytical simple expressions for the angular integrals while the radial integrals are computed numerically. -off–site approach: Here, ˆ VSO is built following the Hemstreet SIESTA 8 formalism52,55 whereby a fully-relativistic pseudopotential (FR-PP) operator is constructed in a fully separable form, i.e., non–local in the radial part as well as in the angular variables, in order to substantially reduce the computational cost. The necessary l j Kleinman–Bylander projectors may be either constructed by SIESTA itself from relativistic semilocal PPs, or directly read from appropriately generated PSML files, as provided by the Pseudo– Dojo project25,26. Moreover, we note that the FR-PP formalism (as well as the original one implemented in Ref. 52) uses the correct normalization constants Cl±1/2, in contrast with what was erroneously stated in Ref. 56. Although we consider the off–site approach more accurate, as it includes inter-shell and inter-atomic SO couplings, both approximations yield very similar results in most of the tested systems, with relevant qualitative differences only found in a few specific cases. Furthermore, the construction of the VSO µν matrix is very fast under both schemes and involves a tiny fraction of the entire self-consistent calculation. G. New electronic-structure solvers For most problems, SIESTA spends the largest fraction of cpu-time in the solver stage (solution of the generalized eigenvalue problem HΦ=εSΦ). The stage devoted to the calculation of the hamiltonian H and overlap S is typically much lighter weight, as those matrices are intrinsically sparse due to the use of a finite-support basis set. Accordingly, SIESTA’s performance is almost completely linked to the use of appropriate external solver libraries. Over the past few years we have expanded the choices available to users and refined the relevant interfaces. Initially, we added support for new individual solvers as detailed below, but recently we have consolidated some of the most important functionality under a new common interface to the ELSI library of solvers57,58. 1. Solvers with a native interface Diagonalization (solution of the generalized eigenproblem appropriate for non-orthogonal orbitals) is the default method for obtaining the density-matrix in SIESTA. A number of standard routines are contained in the SCALAPACK library59, but more efficient alternatives are possible. In particular, the ELPA library60–62 uses an extra intermediate step in the tridiagonal conversion of the matrices to obtain better scalability and significant speedups over SCALAPACK. An interface to ELPA is offered in SIESTA, so this solver can be used as a drop-in replacement for SCALAPACK throughout the code. In addition, SIESTA has implemented interfaces to several methods not based on diagonalization. In most cases, the use of a finite-support basis set, leading to the appearance of sparse matrices, is a significant factor to achieve good performance: •The Fermi Operator Expansion method (FOE)63 uses the formal relationship between Hamiltonian and density-matrix, ˆ ρ=fFD(ˆ H−µ), where fFD is the Fermi-Dirac function. A simple polynomial expansion of fFD can then be used to obtain ˆ ρwithout diagonalization. This method is implemented in the CheSS library64, developed within the BigDFT project65. •The PEXSI method66,67 uses a pole expansion of fFD to get ˆ ρin the form: ˆ ρ=Im P ∑ l=1 ωρ l H−(zl+µ)S!(4) where ωρ land zlare the weights and poles for the corresponding expansion of the Fermi-Dirac function. The number of poles needed is significantly smaller than for the polynomial version of the FOE, as its dependence on the spectrum size is only logarithmic. 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 computed. This “pole expansion and selected inversion” method offers a reduced complexity (at most O(N2) for dense systems, and O(N)for quasi-one-dimensional systems), and trivial parallelization over poles, so it is well-suited for very large problems on large machines. For example,68 computed the electronic structure of large (up to 11,700 atoms) graphene nanoflakes using SIESTA-PEXSI. •The electronic structure problem can also be cast as a minimization problem (of an extended functional) without orthogonalization. When additional localization constraints are put in place, the original linear-scaling method in SIESTA results. Without the extra localization constraints, the cubic-scaling Orbital Minimization Method (OMM)69 can be competitive with respect to diagonalization, as data can be reused across scf-cycle steps. 2. The ELSI interface We have considerably extended the range of solver choices and the performance enhancement possibilities of the code with the integration of the open-source ELSI library (https: //elsi-interchange.org), that provides a unified software interface that connects electronic structure codes to various high-performance solver libraries to solve or circumvent eigenproblems encountered in electronic structure theory57. ELSI also ships with its own tested versions of the individual solver libraries, but additionally, linking against already compiled upstream versions from each solver library is supported as much as possible. The ELPA, OMM, and PEXSI solvers, which had their own ad-hoc interfaces as described in the previous section, SIESTA 9 are now available through ELSI, which also supports other conventional dense eigensolvers (EigenExa70, MAGMA71), sparse iterative eigensolvers (SLEPc72), and linear scaling density matrix purification methods (NTPoly73). As sketched in Fig. 6, an electronic structure code interfacing to ELSI automatically has access to all the eigensolvers and density matrix solvers supported in ELSI. In addition, the ELSI interface is able to convert arbitrarily distributed dense and sparse matrices to the specification expected by the solvers, taking this burden away from the electronic structure code. A comprehensive review of the capabilities in the latest version of ELSI, including parallel solution of problems found in spinpolarized systems (two spin channels) and periodic systems (multiple k-points), scalable matrix I/O, density matrix extrapolation, iterative eigensolvers in a reverse communication interface (RCI) framework, has recently been completed58. 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 various eigensolvers and density matrix solvers via the ELSI API. Whenever necessary, ELSI handles the conversion between different units, conventions, matrix formats, and programming languages. With the common interface in place, any additions and enhancements to the supported solvers can be used in SIESTA with almost no code changes. This is particularly relevant for performance enhancements. For example: •Further levels of parallelization: A feature common in principle to all solvers is that the SIESTA-ELSI interface can exploit the full parallelization over k-points and spins mentioned above. This means that these calculations can use two extra levels of parallelization in the solver step beyond the standard one of parallelization over orbitals (see Fig. 7). •The new version of the PEXSI solver integrated in ELSI can achieve the same level of precision with fewer poles, and offers an extra level of parallelization over trial points for the determination of the chemicalpotential. •Mixed-precision support: The ELPA solver can be invoked in single-precision mode, which can speed up the initial steps of the electronic self-consistent-field (scf) cycle. •Accelerator offloading: The ELPA library offers GPU support in some kernels62, and there is scope for extending it to more kernels. ELSI also offers an interface to the accelerator-enabled MAGMA library. Finally, the PEXSI developers are working on adding GPU support to the solver. 32 64 128 256 512 1024 64 128 256 512 1024 CPU time (s) MPI ranks Scalapack Siesta-ELPA ELSI-ELPA multi-k-point Ideal FIG. 7. Performance improvement from the use of the extra level of parallelization over k-points in SIESTA using the ELSI interface with the ELPA solver, compared to the previous diagonalization scheme (using both the standard SCALAPACK solver and the existing ELPA interface in SIESTA). The system is bulk Si with H impurities, with 1040 atoms, 13328 orbitals, and a sampling of 8 k-points. The multik scheme is able to stay closer to ideal scalability for larger numbers of MPI processes. H. Time dependent DFT Time-dependent density-functional theory (TD-DFT) was first implemented into SIESTA in its real-time propagating form. It was first described in Ref. 74, and then briefly in Ref. 2. It was based on the Crank-Nicolson algorithm, by which, the effect of the evolution operator for an infinitesimal time step ˆ U(t0+∆t,t0) = exp−iˆ H(t)∆t(5) on the wave-function coefficients matrix at a given time t0, c(t0)is approximated by c(t0+∆t) = S+iH(t0+∆t)∆t 2−1S−iH(t0)∆t 2c(t0) (6) where ∆trepresents the finite time-step resulting from time discretization, and Sand Hrepresent the overlap and Hamiltonian matrices, respectively, in the representation given by SIESTA 16 ing MLWFs have been proposed, being another interesting research line for the future.113 L. Multiscale methods Density Functional Theory can be used as the basis for parameterized multiscale methods, that can be used to carry out simulations including tens or even hundreds of thousands of atoms.114 First-principles methods are used to produce detailed models that are subsequently used to predict properties that require large-scale simulations. The models are created for specific materials and their accuracy can be systematically improved to converge towards DFT precision. Given the dependence on first-principles, we refer to these methods as second-principles DFT (SPDFT) and are run on an independent code called SCALE-UP.114 SPDFT is based on a division of the total electronic density, n(~r), into a reference (n0(~r)) and a deformation (δn(~r)) contributions, n(~r) = n0(~r)+δn(~r),(27) where δn(~r)is considered as a small perturbation with respect to n0(~r)that, in non-magnetic cases, represents the ground state of the system114. This division is then used114 to expand the DFT energy with δnfinding that the zeroth order term, E(0), corresponds with the full DFT energy for the reference density. The corrections to this reference energy only depend on δn(and parametrically on n0) which, given its smallness, can be efficiently calculated leading to a fast and accurate approximation of the full DFT energy. The expansion is usually taken to second-order, E≈E(0)+E(1)+E(2)+..., (28) resulting in a stationary problem that is equivalent to HartreeFock with the important distinction that the interactions are screened by the exchange-correlation potential. In order to keep δnsmall the application of the method is restricted to problems where atomic bonds are not created or destroyed, i.e. to processes that display an invariant bond topology. The E(0)term represents the exact DFT energy for the reference density. We represent it for a variety of geometries with an accurate force-field115 that allows for fast evaluation. The E(1)and E(2)terms account for the changes in the electronic structure that are represented by geometry-dependent Wannier functions. Under this basis E(1)becomes a tight-binding model while E(2)represents electron-electron interactions. The interconnection between the first (SIESTA) and the second (SCALE-UP) principles simulations is carried out through a python script, MODELMAKER. Taking a few cutoff distances MODELMAKER is able to produce a model’s terms and automatically carry out DFT simulations with SIESTA to determine the force field, a Wannier Hamiltonian to represent the bands, electron-lattice terms to account how the bands change with geometry and electron-electron interactions to describe, for example, magnetism. While, so far, few publications with SPDFT methods include explicit treatment of electronic degrees of freedom, the lattice part has successfully been used in several applications. One of the main fields of research has been thermal conductivity in perovskites. In particular it was employed to study the electrophononic coupling in SrTiO3116 and PbTiO3117 and the proposal of a thermal switch in PbTiO3.118 It has also been used to study the competition between various ferroelectric domain structures in PbTiO3/SrTiO3superlattices as a function of strain.119 As a result it was found that tensile strains lead to the appearance of chiral ferroelectric vortices while ferroelectric skyrmions were predicted and experimentally observed for more compressive strain values.119 The calculated dielectric properties of these superlattices120 are in very good agreement with measured values and show very large electric susceptibility consistent with regions of negative, static electric permittivity situated at the core of the vortices and the PbTiO3/SrTiO3interfaces. M. Scripting and integration in external frameworks An ongoing trend in many areas of computational science is to move away from rigid and monolithic codes, favoring instead a more flexible approach in which the internal functionality of a program is somehow exposed to the outside world. If done in a proper and well-documented way, this can serve to enhance the interoperability of codes with different functionalities, playing to the relative strengths of each, and/or to implement new functionalities by combining the available basic blocks. In SIESTA we have followed two different but complementary routes to these ends: the development of an internal scripting framework based on the Lua language, which enables new functionality without code recompilation, and the implementation of a formal interface to the AiiDA platform. 1. Lua interface Lua121 is an easy-to-learn and fast scripting language built for embedding. It is very lightweight (its memory footprint is less than 300kB), and provides very simple ways to interface to the data structures and routines of a host program. A Lua script, interpreted by the Lua interpreter embedded in the program, can then control the flow of execution and the data. Different user-level scripts can implement new functionalities, without recompilation of the host code. The strategy we have followed in SIESTA is based on handling control to the Lua interpreter at specific relevant points in the program flow (e.g. at the beginning of a geometry step, at the end of a scf step, etc). Lua scripts implement handlers appropriate to the point they want to hook into, and can request access to specific data structures. For example, a script intended to implement a better scf mixing algorithm would be executed after every scf step, inspecting the convergence data, and changing mixing parameters or schemes, as appropriate. As another example, convergence checks over mesh-cutoffs and k-point sampling can be performed automatically. The above mixing scenario exemplifies an important area of usefulness of the approach: the prototyping in Lua, (followed SIESTA 17 eventually by a full implementation), of new ideas and algorithms. We have implemented a number of custom molecular dynamics modes, geometry relaxation algorithms, and advanced optimization schemes, in a pure Lua library FLOS122. The code in the library can be re-used, or taken as starting point for other implementations by users. These user-level scripts can in turn be shared, opening the way to the development of new functionality with faster turnaround that the traditional approach that needs a careful integration into the program’s code base. As a specific showcase of the power of the Lua embedding, we have developed a number of variations of the nudgedelastic band method (NEB)123,124 for transition-state search. Previously proposed implementations in SIESTA involved significant, 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, but a Double Nudged Elastic Band (DNEB)125 variation, and also another version which treats atomic coordinates and lattice variables on an equal footing (the variable-cell NEB, or VC-NEB, method126). The integration of Lua functionality in SIESTA has been made possible by the development of an intermediate layer, FLOOK127, (for “fortran-Lua-hook”), which provides wrappers for access to Fortran data structures and subroutines. 2. AiiDA plugins and workflows The AiiDA framework128–130 provides support for highthroughput computations in materials science, keeping full provenance of the calculations and facilitating data handling and sharing. The framework is open-source, written in Python, and designed to support arbitrary codes via a plugin interface. A plugin for SIESTA has been implemented and is distributed as the open-source package aiida-siesta131. The plugin provides the basic operations of preparing the input files for a calculation using AiiDA-specific input objects, and parsing the results and generating AiiDA output objects. The AiiDA data are stored in a graph database that keeps a permanent record of the inputs and outputs of the calculation, and is fully searchable for, e.g. data analytics purposes. AiiDA also provides robust support for the creation of workflows that incorporate all the necessary steps in the calculation 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 standard materials properties, such as band structures. Fig. 15 shows the execution graph of a workflow designed to generate a synthetic STM image from a given structure. Work is ongoing to implement more complex ones. In addition to an interface to the computational capabilities of the SIESTA code via the plugin and workflows, the aiida-siesta package also provides an implementation of basic objects representing pseudopotential files, notably one for PSML. Families of pseudopotentials can be uploaded to an AiiDA database and shared via the provided mechanisms for data export and import, facilitating the interoperability of Code (1603) 'siesta-4.0-new' WorkCalculation (3150) code Bool (3151) CREATE Int (3152) CREATE ParameterData (3153) CREATE ParameterData (3154) CREATE KpointsData (3155) 2x2x2 (+0.0,0.0,0.0) CREATE ParameterData (3156) CREATE ParameterData (3157) CREATE WorkCalculation (3158) CALL STMCalculation (3170) FINISHED CALL ArrayData (3174) stm_array StructureData (3166) CH output_structure Str (3148) protocol Code (2845) 'plstm-4.0' stm_code StructureData (3149) CH structure clean_workdir max_iterationsparameters basis kpoints SiestaCalculation (3162) FINISHED kpoints settings options output_structure CALL ParameterData (3165) output_parameters RemoteData (3163) remote_folder FolderData (3164) retrieved stm_array RemoteData (3171) remote_folder FolderData (3172) retrieved ParameterData (3173) output_parameters output_structure output_parameters remote_folderretrieved ArrayData (3167) output_array parent_calc_folder STM image of benzene molecule FIG. 15. Automatically generated graph for the execution of an AiiDA workflow for simulation of STM images. different codes. N. Utilities for post-processing and supplementary features SIESTA offers several features beyond the core functionality 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 basis set enables the use of a very intuitive suite of analysis tools, which take advantage of the fact that most of the concepts relating to chemical bonding use the language of atomic orbitals. The (partial) density of states, atomic and orbital populations, and other useful output can be obtained directly from the program. The SIESTA distribution includes also several tools in the Util directory for band-structure and wavefunction plotting, bonding analysis, etc. Beyond these, special tool packages that implement a specific feature that extends the functionality of the main program, or that provide extra options for visualization or post-processing in general, are available in alternate distribution points. We describe in what follows the most relevant developments. 1. Updates to core utilities A number of improvements, enhancements, and additions have been made to the core utilities shipped with the SIESTA distribution. There is now a “fat-bands” feature, by which bands can be decorated with information about the relative weight of given orbitals in each state. The wave-function-related analysis tools have been extended to the non-collinear and spinorbit case. This includes the COOP/COHP bonding analysis, band-structures, and a new tool spin-texture calculation. There have been also improvements to band-structure plotting SIESTA 18 utilities and to the visualization of charge densities, potentials, and other magnitudes represented in a real-space grid. A band unfolding utility has been added. Based on the Fourier decomposition of the Bloch wave functions, it allows to perform a “full unfolding” even for non-periodic systems (e. g. liquids) calculated with a large simulation cell. By refolding the fully unfolded bands, from the reciprocal supercell of a perturbed or defective crystal, into the reciprocal unit cell of the primitive crystal, one recovers the conventional “unfolded” bands132. 2. sisl SISL is a Python toolbox that was initially conceived to handle and manipulate SIESTA/TRANSIESTA output103. It has since been extended to support other DFT codes, with the aim of offering equivalent operations for them. By reading the LCAO outputs from SIESTA one can postprocess the Hamiltonian and calculate e.g. Brillouin zone integrated DOS, wave functions expanded on grids, eigenvalues, band velocities and many more. SISL can process nearly all the SIESTA output files. In particular, it is also able to post-process data on the real-space grid. Its command line interface allows data format changes, e.g. conversion of SIESTA XV files to xyz/xsf files or SIESTA binary grid data (VH,VT, ...) to cube/xsf files. As it can process density matrices from SIESTA, one can also use SISL to prepare an input electronic-structure for new calculations, which may be helpful to reduce initial SCF steps. SISL also allows creation of custom tight-binding models (both orthogonal and non-orthogonal), and since it extracts the DFT Hamiltonian matrix one can manipulate the Hamiltonian to retain certain band-structure features and thus perform large-scale simulations133. This allows calculating far-field currents using reduced basis-sets with very little loss of accuracy. The Atomic Simulation Environment (ASE)134 and SISL have a certain degree of overlap in terms of geometry handling functionality. One can easily convert to and from ASE objects in SISL, thus allowing seamless interaction. 3. Other post-processing and visualization utilities The body of utilities contributed by non-core developers and other SIESTA users has continued to expand. In particular, we feature in this section two suites of utilities, one dealing with alternate visualization tools for some SIESTA results, and another one specifically dealing with lattice dynamics. For structures, the xv2xsf and xv2vesta converters process data from the SIESTA .XV file into the native formats of XCrySDen135 and VESTA136, respectively. Each of these two codes offers many options of graphical representation of structures, adding translations, clipping fragments etc. Threedimensional spatial functions (e.g., charge density, local density of states integrated throughout the chosen energy range), computed by SIESTA on a real-space grid. Tools are provided for interpolating the data from the SIESTA output grid (fixed by the unit cell dimensions and the MeshCutoff parameter) onto an arbitrarily cut (and possibly rotated or resampled) parallelepipedic box. XCrysDen provides a number of display options, including contour lines over grid planes, or isosurfaces. A special feature available in XCrySDen is plotting the Fermi surfaces. A special script, eig2bxsf, serves to analyze the list of k-points handled by SIESTA, expanding it onto a regular sequence, and writing the respective band energies in the necessary format. The tools concerning the lattice dynamics have been developed having in mind the Γphonons calculated for a large enough supercell, that is a typical case in a simulation of molecular crystals or disordered substitutional alloys. For visualization, vib2xsf and vib2vesta place arrows at the atoms according to the vibration pattern stored in the eigenvectors file (.vectors), produced by the core Vibra utility, and can also be used to make animations (sequences of snapshots) of selected vibration modes. Both vib2xsf and vib2vesta tools allow the selection of a part of the system to be exposed. The phdos tool is designed for analyzing zone-center vibration results. As the system is supposed to be large (e.g., a supercell chosen for a periodic crystal), the (artificially broadened, for convenience) discrete spectrum may serve as a fair approximation to the total density of modes, and if weighted with (squared) components of eigenvectors at different atoms – provide a decomposition into contributions of different atoms in the total density of vibration modes. A more sophisticated option is the projection of different eigenvectors according to various criteria. The typical system under study is a supercell in which e.g. an alloying, or some kind of deformation, breaks the underlying perfect periodicity. Still, some trends related to the latter can be revealed by appropriate projections. The two obvious cases are the projections onto (1) q-vectors of the underlying lattice and (2) irreducible representations of the space group of the underlying lattice; the corresponding formulas and some results can be found in Ref. 137. The first type of projection, if done for a sequence of qvalues, helps to reveal “phonon dispersions”, obviously blurred by the broken periodicity, also making distinction between transversal and longitudinal modes – see Ref. 138 for an example of use. To make the trends more pronounced, the supercell needs to be sufficiently long in the direction concerned – see, e.g., Fig. 16. The simplest case, a projection onto a single q=0 value, may also be of interest, since it enhances the modes which are expected to dominate the infrared or Raman spectra, and thus facilitates their comparison with experiment. The symmetry projection may help to isolate in a possibly complex spectrum those modes which are expected to dominate according to a given selection rule, again in view of their verification against the experiments. The group-symmetry information needed for the projections is available e.g. from the Bilbao Crystallographic Server,139 and the technical details are explained in the documentation included in the tools. The vibent tool performs a straightforward calculation (see, e.g., Sec. II.C in Ref. 140, or Sec. 5.3 in Ref. 141) of SIESTA 19 Be Zn Se 0 100 200 300 400 wave number (cm-1) 0.19 qz (Bohr-1) 0.18 0.17 0.16 0.15 0.14 0.13 0.12 0.11 0.10 0.08 0.09 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 002-reduced Be32Zn64Se96 ; Zn contribution projected vibration DOS for q = ( 0 0 qz) TO LO FIG. 16. Left panel: a 192-at. quasirandom supercell representative for the Be1/3Zn2/3Se solid solution; right panel: density of modes within the frequency range of Zn-Se vibrations, extracted with phdos and projected onto different values of qzand 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. 0 50 100 150 200 250 300 350 400 Arb. units, per supercell Σ (64at.) Cu (16 at.) Zn (8 at.) Sn (8 at.) S (32 at.) Cu2ZnSnS4: 64-at. supercell, total density of modes 0 50 100 150 200 250 300 350 400 Frequency (cm-1) Difference from perfect ZnCu(Sn-pl.) +0 |_ _|Cu(Sn-pl.) +0 |_ _|Cu(Zn-pl.) +0 |_ _|Zn +0 |_ _|S +0 ZnCu(Zn-pl.) +0 CuZn +0 ZnSn +0 0 200 400 600 800 1000 Temperature (K) -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 Vibrational free energy (eV) Vibrational properties of point defects in CZTS: doped minus perfect 0 200 400 600 800 1000 Temperature (K) -0.4 -0.2 0 0.2 0.4 Vibrational entropy (meV/K) ZnCu(Sn-pl.) +0 ZnCu(Zn-pl.) +0 CuZn +0 ZnSn +0 |_ _|Cu(Sn-pl.) +0 |_ _|Cu(Zn-pl.) +0 |_ _|Zn +0 |_ _|S +0 FIG. 17. Vibration properties of Cu2ZnSnS4(CZTS) with substitutional impurities, used in Ref. 140. Left panel: densities of modes (extracted with phdos); right panel: vibration contributions to the free energy and entropy (calculated with vibent). Adapted from Fig. 5.3 and 5.4 of Ref. 141. temperature-dependent vibration contributions to the free energy and entropy – see Fig. 17 as an example. The necessary input information is the vibration spectrum, originating from the Vibra frozen phonon calculation on a sufficiently large system. The velcf tool calculates the velocity autocorrelation function and its Fourier transform from a (presumably sufficiently long) molecular dynamics (MD) history, recorded in the .MD or .ANI file. This technique142 can be used to obtain phonon frequencies, and was applied along with a SIESTA calculation in Ref. 143. An example of such simulation (1000 0 500 1000 1500 2000 2500 3000 3500 4000 Frequency (cm-1) Ni H Mo O "Ni4": frozen phonon calculation 0 100 200 300 400 500 600 700 800 900 1000 MD steps -0.05 0 0.05 0 500 1000 1500 2000 2500 3000 3500 4000 Frequency (cm-1) 0 0.002 0.004 0.006 cos - VCF "Ni4": molecular dynamics at 600 K Velocity autocorrelation function FIG. 18. Vibration properties of “Ni4” molecular magnet, [Mo12O30(µ2-OH)10H2{Ni(H2O)3}4]·14H2O. Left panel: density of modes from frozen phonon calculation; right panel: velocity autocorrelation function, its Fourier transform and hence resulting density of vibration modes. MD steps at 600 K) is shown in Fig. 18 in comparison with frozen phonon results, revealing similarities of the spectra obtained. 4. Optical properties of finite systems: linear response TDDFT starting from Siesta orbitals The SIESTA package offers at least two ways of obtaining optical properties of finite systems. The first way uses real-time TD-DFT propagation by applying an external electric field with a simple time dependence (e.g., a Heaviside step-function)74. The second way is by computing the noninteracting dielectric function1,2. Both methods are implemented in SIESTA and can be employed without any external tools. However, they are limited in different aspects. The non-interacting dielectric function often underestimates the HOMO-LUMO gaps and calls for the use of the phenomenological scissor-shift operator. Real-time propagation makes cumbersome the analysis of the optical response properties in the frequency domain. Furthermore, the frequency resolution scales with the duration of the real-time simulation. Thus, accurate spectra require long simulations. Fortunately, there are two efficient implementations of linear-response TDDFT that use the Kohn-Sham orbitals from SIESTA as a starting point and are available for the opensource community144,145. In both packages, the linear density response δn(r,ω)is obtained directly in the frequency domain which makes straightforward the analysis of derived properties. However, there are differences between both implementations on the construction of the auxiliary basis necessary to expand the orbital products. These differences can severely affect the computational cost of the calculation. The linear-response TDDFT is built on the concept of the induced electronic density δn(r,ω)in response to a small perturbation of the external potential δVext(r,ω). The integral operator connecting δn(r,ω)to δVext(r,ω)is the interacting density response function χ(r,r0,ω). By virtue of the KS equations, χ(r,r0,ω)can be connected to the non-interacting density response function χ0(r,r0,ω)146 with a Dyson equation SIESTA 20 χ(ω) = χ0(ω)+ χ0(ω)Kχ(ω),(29) where the interaction kernel K(r,r0)contains the bare Coulomb interaction and the so-called exchange and correlation kernel Kxc, which is a known operator for simple functionals like LDA and GGA. The non-interacting response function χ0(r,r0,ω)can be expressed as a sum over electronhole excitations within the basis formed by the KS orbitals Ψn(r)145–147 χ0(r,r0,ω) = ∑ nm (fn−fm)Ψn(r)Ψm(r)Ψm(r0)Ψn(r0) ω−Em+En , (30) where fnare occupations of the KS orbitals and Enare their energies. The optical polarizability tensor α(ω)is related to the induced density by α(ω) = Rrδn(r,ω)dror alternatively α(ω) = ZZ rχ0(r,r0,ω)δVs(r0,ω)drdr0,(31) where due to Eq. (29) and using the dipole approximation for the electron-photon coupling, the screened effective perturbation δVs(r0,ω)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 χ(ω) from Eq. (29). Once δVs(ω)is known, Eqs. (30) and (31) allow the computation of the optical properties of the system. 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 to existence of the sum over the electron-hole pairs in Eq. (30). Similarly, one can achieve other types of Mulliken-like analysis145,147,149 of the optical polarizability tensor α(ω)or the induced density δn(r,ω). The Python implementation of linear response TDDFT in the PySCF-NAO package as described in Ref. 145 is convenient to use and rather potent. It is capable of computing the optical properties of compact metallic objects containing up to several hundreds of atoms147,150,151. For example, we were able to track down the different size-dependence of the plasmon resonance in sodium and silver clusters due to the screening effect of silver d-orbitals in the latter case152. In those calculations, using an optimized version that incorporates some additional memory-saving features not present in the currently distributed version of PySCF-NAO, icosahedral silver and sodium clusters containing up to 5043 atoms were studied. In Figures 19 and 20, we show the photo-absorption cross sections of a series of compact silver clusters147 and the real part of induced density change in the cluster Ag147 close to its surface-plasmon frequency (3.4 eV), respectively. 0 0.4 0.8 1.2 1.6 2 2 3 4 5 6 7 8 σopt (nm2) ω (eV) Clusters Ag561 Ag309 Ag147 Ag55 Ag13 FIG. 19. The absorption cross sections of silver clusters of icosahedral shape. One can recognize sharp surface-plasmon resonances around 3–4 eV and a broad resonance at 6–7 eV. FIG. 20. The isosurfaces of density change Re(δn(r,ω)) of the Ag147 cluster close to the frequency of the surface-plasmon resonance of the cluster (3.4 eV). 5. Thermal transport by the AEMD method The approach to equilibrium molecular dynamics (AEMD) method153 has been implemented to obtain the thermal conductivity. In the first stage of the method, the system is decomposed in two different regions, each one equilibrated to a different initial temperature (canonical run with Bose, or Anneal MD). Then, a microcanonical run (Verlet) is carried out for the whole system, and the average temperature of each subsystem is monitored. This temperature transient regime is then used to extract the thermal conductivity from the exact solution of the heat transport equation.154 6. Core level shifts Core-level shifts can serve to analyze changes in the local and chemical environment of atoms of a given species. Density-functional-theory calculations have proved to be quite useful in complementing the experimental information, which is sometimes hard to interpret. Two schemes have been SIESTA 21 implemented in SIESTA for the calculation of core-level shifts within a pseudopotential approach155. In the so-called initial-state approximation the electronic relaxation in the presence of the core hole is neglected, and the photo-electron’s binding energy is directly related to the eigenvalue of the core level. A pseudopotential calculation obviously cannot compute the latter, but differences in core eigenvalues in different environments can be estimated by the changes in the expectation value of the crystal potential using the core state’s atomic wavefunctions ψlm nat different sites. These can be extracted from the matrix elements Vmm0=Zd3r(ψlm n(~r−~ τ))∗V(~r)ψlm0 n(~r−~ τ)(33) with a further step of averaging to remove the splittings stemming from the loss of spherical symmetry. In the final-state approximation, the relaxation is explicitly taken into account, and the experimental shifts (measured via the kinetic energy of an exiting electron) are correlated with the differences in the energy of the crystal with a “core-hole” in different sites. For this, a special pseudopotential with a missing core electron has to be generated, and a full SIESTA calculation is needed for each different site. The implemented methodology has been used to study, for example, the shifts induced by hydrogen bonding in organic molecules156. O. Software-engineering advances and partnerships The traditional development model for scientific codes in academic settings has been typically based on multiple contributions with various levels of programming competence, and with very little time to plan ahead in the face of pressing scientific demands. SIESTA has been no exception, and has grown in features and complexity over the years. It is very important to keep complexity under control, or else a project becomes un-maintainable and cannot survive. It is not simple, however, to balance the need of incorporation of new features, and the need to increase the computing performance in a landscape of constantly evolving hardware and programming models. One essential route is modularization, which allows the separation of concerns at various levels. In the context of a code like SIESTA , this means that the scientific ideas and algorithms should be handled at a high level, calling on lower-level modules for specific functionality (domainspecific libraries, mathematical libraries, communication protocols, etc). These lower-level modules can hopefully be reused by different codes and, most importantly, can be focused on by highly-skilled programmers for optimization on relevant architectures. Another important method of taming complexity involves the streamlining of the data structures of the code. This is an ongoing process (see Sect. V), but has already taken a very significant step by the introduction of reference-counted data structures. They build on a well-known and not particularly advanced technique of memory-handling157, but in SIESTA they have enabled a much simpler bookkeeping of the data structures needed for a richer control of molecular-mechanics and scf iterations. Regarding performance-oriented developments, in the recent past we have implemented a mixed MPI/OpenMP programming model, which allows, for suitable systems, to better balance arithmetic intensity and communications needs. The deployment of this model is more advanced in the TranSiesta module, and significant speedups have been obtained for large systems. Some of the above software-engineering developments have been enabled and strengthened by the participation of SIESTA in a number of international partnerships, notably the MaX (Materials at the eXascale) EU center of excellence158 and the Electronic Structure Library initiative159. The “separation of concerns” described above in the context of modularization is an example of the so-called “open-innovation” paradigm, at the foundation of the ESL strategy for code reusability, and is also a cornerstone of MaX’s efforts to achieve exascale-readiness for its flagship materials science codes (with SIESTA among them): performance-enhancement efforts are to be focused on relevant domain-specific modules. A number of modules from SIESTA have been turned into stand-alone libraries which now feature in the ESL: libGridXC for exchange and correlation calculations, libPSML as a handler of PSML files, xmlf90 for general purpose handling of XML files, etc. Conversely, SIESTA uses some of the libraries offered by the ESL, notably the ELSI library of electronic-structure solvers mentioned in Sect. IIIG 2, whose development, including its API design and internal data organization, has been in turn influenced by contributions and feedback from the SIESTA project, among others. There are also plans to incorporate the PSolver library160 for the solution of the Poisson problem, a contribution to the ESL from the BigDFT project. We should mention that the renewed dynamism of SIESTA development and the advances made possible by the interaction with community initiatives are both a blessing and a challenge. It is non-trivial, for example, to handle the building process of a code that relies on a number of different external libraries, programming models, and special features such as the embedded Lua interpreter. Luckily, as will be discussed in Sect. V, these are issues that are being addressed in wider contexts, and SIESTA is well placed to take advantage of it. IV. APPLICATIONS We present here a few showcase applications that illustrate the capabilities of SIESTA, in breadth, efficiency, and accuracy. A. 4 terminal NEGF on germanium surface Breakthrough simulations using the new multi-terminal implementation on TRANSIESTA were fundamental to elucidate the electronic transport mechanism on a novel and complex SIESTA 22 FIG. 21. First-principles transport simulations for the two-probe experiments. a) Representation of the four-terminal setup. The electrode regions are highlighted by blue boxes, two of them located at each Ge(001)-c(4×2) slab terminations (leads left and right) and the 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 Ref. 101. experiment.101 For the first time a two-probe scanning tunneling microscopy/spectroscopy (STM/STS) with probes operating in tunneling conditions over the same atomic-scale system was used to extract detailed information of in-plane electronic transport. The addressed system was the reconstructed (001) surface of germanium, where electrons injected from one STM tip at a position determined with atomic precision were collected at the same Ge dimer row at a distance as short as 30 nm. The experiment was theoretically modeled by a system composed of a twelve-layer Ge(001)-c(4×2) slab contacted by Au tips oriented along the (100) direction (Fig. 21). On this self-consistent 4-terminal treatment, two Ge electrodes were connected at each slab termination and other two at the Au model tips. The whole system was defined by 4924 atoms (36442 atomic orbitals), in a super-cell of dimensions ∼32×160×80Å3, and where 5 different tipto-sample distances were considered. Besides the large dimensions of the system, another important challenge of such simulation was the level alignment between the metallic and semiconducting leads and the scattering region, for which a method had to be devised. A remarkable agreement was found between the calculated transmission function and the experimental transconductance spectra, allowing the identification and assignment of the observed resonances to transport channels existing along the surface Ge dimer rows. Moreover, the simulations elucidated the transport directionality of the injected hot electrons, revealing a transition from 2D to quasi1D coherent transport regime as a function of the carrier’s energy. This work shows that complex experiment setups combined with advanced calculations can provide new insights into transport properties at the nanoscale. B. Novel topological phases in ferroelectric materials In material systems with several interacting degrees of freedom (such as spin, charge and lattice distortions), the complex interplay between these factors can give rise to exotic phases. A prototypical example are the superlattices of alternating lead titanate and strontium titanate layers. Simulations on such PbTiO3/SrTiO3heterostructures, consisting on nunit cells of PbTiO3and nunit cells of SrTiO3stacked along the [001] direction were carried out with SIESTA. As a function of the periodicity, the superlattices undergo a phase transition from a monodomain configuration (small periodicity, n.3−4) with a normal component of the polarization that is preserved throughout the structure, to a multidomain configuration (large periodicity, n&3−4) with alternating up and down domains.161 In order to further reduce the electrostatic energy costs, the local dipoles within the PbTiO3layer continuously rotate forming a sequence of clock-wise/counterclockwise array of vortices along the [100] direction. The theoretical predictions, done with SIESTA162 after the relaxation of supercells of up to 1000 atoms, were experimentally confirmed five years later by atomic-scale mapping of the polar 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. C. 1D and 2D systems SIESTA is particularly well suited to study low dimensional nanostructures, such as 1D and 2D systems where a large vacuum region is needed within the simulation cell. When, in addition, a large number of atoms is required to study particular physical effects is where SIESTA could excel with respect to other methods. There is extensive literature on simulations of graphene and other exfoliated materials, where the properties of point defects, edges or grain boundaries are of much relevance. To list a few examples, the magnetic properties of impurities,165,166 and edges167, but also electronic properties, including transport characteristics, in grain boundaries168,169, ribbons170, nanoporous graphene171, large graphene flakes68,172, or the effect of substrates173. Other materials, such as monoand multi-layered dichalcogenides174,175 or phosphorene176,177, are also being widely studied, including optical properties in nanoflakes with up to a few thousand atoms.178. 1. CDWs A number of recent studies on charge density waves (CDW) in low dimensional materials illustrates the impressive accuracy that can be obtained with SIESTA for systems with very subtle electronic structures.179 For example, in 2H-NbSe2 SIESTA calculations were able to predict the existence of six different atomic structures within a narrow energy range of a SIESTA 23 12 INVESTIGACIÓN Y CIENCIA, agosto 2019 Panorama Todos los aparatos electrónicos que nos rodean se basan en el uso de transistores: diminutos dispositivos de mecanismo aparentemente simple. Un transistor puede entenderse como un pequeño interruptor con dos posiciones: «encendido» (permite la circulación de corriente) y «apagado» (no la permite). Una combinación adecuada de ellos permite implementar cualquier operación lógica. Su descubrimiento, efectuado en los años cuarenta del siglo pasado por John Bardeen, Walter Brattain y William Shockley en los Laboratorios Bell, fue reconocido en 1956 con el premio Nobel de física y supuso el inicio de la revolución tecnológica de la que disfrutamos hoy. Desde sus comienzos, la industria microelectrónica se ha obsesionado con integrar un número cada vez mayor de transistores en los circuitos. Esta carrera por la miniaturización se resume en la conocida ley de Moore, enunciada en 1965 por Gordon Moore, cofundador de Intel, y según la cual el número de transistores en un circuito se duplica aproximadamente cada dos años. En la actualidad, las dimensiones de estos componentes electrónicos se sitúan entre los 15 y los 50 nanómetros, una escala en la que los efectos cuánticos comienzan a ser relevantes. Como consecuencia, un procesador moderno puede llegar a albergar miles de millones de transistores. Sin embargo, ese aumento en el número de componentes por circuito tiene ELECTRÓNICA Condensadores con capacidad negativa Un estudio encuentra el origen microscópico de la capacidad negativa, una exótica propiedad electrónica que aparece en ciertos materiales. El hallazgo augura el diseño de nuevos transistores más eficientes Pablo García-Fernández y Javier Junquera CORTESÍA DE RAmAmOORThy RAmESh, UNIVERSIDAD DE CALIFORNIA EN BERKELEy REMOLINOS DE POLARIZACIÓN: Mapa microscópico de la polarización (flechas amarillas) en capas de titanato de plomo y titanato de estroncio. La inusual reacción de estos materiales ante un campo externo genera zonas donde la capacidad eléctrica es negativa, un fenómeno considerado imposible hasta hace pocos años. FIG. 22. Top panel: local polarization profile of polydomain structures in (PbTiO3)n/(SrTiO3)nwith n=6 obtained from an atomic relaxation with SIESTA. The PbTiO3and SrTiO3are depicted as grey and white regions respectively. Clockwise and counterclockwise vortices within the PbTiO3are clearly visible. Red dashed square in the SrTiO3layers mark the position where antivortices are formed. Reprinted with permission from Aguado-Puente and Junquera162 Phys. Rev. B 85, 184105 (2012). Bottom panel: experimental observation of vortex–antivortex structures in a crosssectional high-resolution scanning transmission electron microscopy image with an overlay of the polar displacement vectors for a (SrTiO3)10/(PbTiO3)10 superlattice, showing that an array of vortex–antivortex pairs is present in each PbTiO3layer. Courtesy of R. Ramesh, adapted from Ref. 163. few meV, all of them compatible with the experimental 3×3 CDW modulation. Careful analysis of theoretical and experimental STM images for different bias potentials allowed to identify two of these structures that can coexist in the same image.180 In a different work,181 the temperature dependency of the electronic Lindhard response function in blue bronze K0.3MoO3was studied. This system has a rather complex monoclinic structure, with twenty formula units per unit cell where MoO6octahedra form chains along one direction (baxis). The Lindhard function shows well decoupled sharp responses that correspond to intraand interband Fermi surface nesting. By fitting these peaks one can obtain the coherence length of the fluctuating 1D electron-hole pair (that determines the length scale of the experimental intrachain CDW correlations), and the intrachain modulation of the response (that determines the shape of the Kohn anomaly measured in experiments), providing, for the first time, a quantitative evidence of the weak electron-phonon coupling scenario for the Peierls transition. D. Siesta in biology: pilin proteins as conductors SIESTA’s efficiency and the clear band gaps of biomolecules in general have made molecular biology a very suitable field for SIESTA 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 hairpin 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, in solid-state parlance, represents a polarization, with clear electrostatic implications. Indeed, a DFT calculation of the molecule in vacuum shows a well defined electrostatic potential ramp along the protein, which tends to close the effective band gap. The question is then, how does an aqueous environment affect this depolarizing field. Long molecular mechanics (MM) simulations were performed for the protein in a suitable solution of NaCl at a concentration of 0.1 M. The protein’s residues had charge states corresponding to pH =7, and the MM field was validated with SIESTAcalculations in vacuum (944-atom dynamic relaxation in a 104.43Å3box). The wet system contained 4580 atoms, and the statistical average of the electrostatic potential around the molecule (see Fig. 23) was obtained from a sample of full SIESTA calculations of statistically independent snapshots, taken every 50 ps during the last 0.5 ns of the simulation. FIG. 23. Colour coded electrostatic potential on a plane cutting along the main axis of the geobacter sulfurreducens pilin molecule in wet conditions. A perspective ball rendering of the atomic strucuture of the protein is superposed. For the meaning and details on this Figure see Ref. 185 (Figure courtesy of Gustavo T. Feliciano). SIESTA 24 Fig. 23 shows how the aqueous environment kills the quite homogeneous potential ramp along the protein axis that appears in vacuum and replaces it with long-wave-length slow, but quite significant fluctuations. The gap remains sizeable, and coherent transport is not likely. However, the frontier orbitals evolve in a very suggestive way for enhanced diffusive electron transport.185 E. Use of Siesta in other fields Although an exhaustive summary of all the recent results obtained with SIESTA is out of the scope of this work, we would like to point the attention of the reader to a sample of recent reviews in various fields in which the program is featured. These cover biological sciences186 (including interaction between organic and inorganic materials187,188), geology and materials under high-pressure189, isotopic fractionation predictions for Martian geochemistry190, the engineering of typical core structural materials used in nuclear reactors,191 or even in astrophysical and atmospheric systems192. The reactivity of metallic nanoparticles for catalysis was treated by Viñes, Gomes, and Illas193, and the role of SIESTA in the computation of the kinetic and dynamics of catalytic reaction at surfaces (including adsorption and desorption of reactants or products) was explored in Chapter 8 of Ref. 194 by Catapan and coworkers. V. FUTURE EVOLUTION Work on enhancing SIESTA’s capabilities, performance, and robustness is continuing, driven by a good number of developers and collaborators. A mature and flexible development platform and practices are essential to keep them productive. Our recent platform changes have forced developers to shift workflows twice in the past four years. Through the changes we have learned a lot but also spent a significant amount of time on ensuring SIESTA’s continuous development. At the current state we believe we have stabilized the development platform on GitLab while we will add more integrated development features in the coming years, e.g. continuous integration (CI) and source code checks. Using CI will also enable easier code-style checks to conform to coding standards. We hope that our open-platform initiative will keep external contributions coming into the program. Our basic-development plans include also refactoring, apparently unexciting but essential to streamline the code base to enable further implementations. Also, we foresee a change in the release model, moving away from coexisting long-lived release branches whose maintenance takes up a lot of time, and offering instead more frequent and short-maintenance releases. We plan to exploit the idea of modularization, continuing the abstraction of relevant reusable pieces, but also dealing with a higher-level, exposing the core electronic-structure capabilities of SIESTA to other programs. It will be necessary to redesign some of the internal data structures to remove global variables and encapsulate them into objects or derived types associated to particular configurations and stages of the calculations. This encapsulation will be matched by a streamlining of the input/output operations. This work will open the door to the creation of complex workflows leveraging the strengths of various codes. Accelerated hybrid architectures (including, for example, GPUs) are very likely going to feature prominently in the upcoming exascale machines. In the case of SIESTA, the data indirection associated to the handling of sparse matrices limits the acceleration possibilities of the section of the code that builds the Hamiltonian and overlap matrices, but the solver stage is more amenable to porting, and in fact several solver libraries used by SIESTA are being enhanced to offer GPU support, as mentioned in Sec. IIIG 2. Modularization and the use of new programming models cause an increase in the complexity of the building and deployment of the code. We will leverage the ESL bundle, created to facilitate the use of the modules in the ESL collection, to streamline SIESTA’s building process, and explore containerization as an option for deployment of the code. The “pseudopotential barrier to entry” has been lowered by the availability of curated databases supporting the PSML format. Basis sets are a perennial challenge, but new tools and ideas are being explored to provide users with appropriate basis sets: High-throughput workflows for optimization; "tiers" of quality/cost, but perhaps not just of a simple “periodic table” form, as offered by other codes (e.g., FHI-aims13), but with a possible dependence on an approximate characterization of the chemical environment in which a given atom finds itself. Complementary to the underlying basis-set optimization that focuses on providing an adequate variational freedom, an on-the-fly contraction of the basis set, which results in a set of lower-cardinality adapted to the description of the occupied subspace can be exploited for increased efficiency. This is particularly relevant for FOE methods (see Sect. IIIG 1, in which the number of polynomial terms depends on the extent of the spectrum. The original claim to fame of SIESTA was based on its linear-scaling solver. We are in the process of a re-design of the O(N)code with a new, more efficient backend, based on the DBCSR library for handling distributed block-sparse matrices195,196 with the MatrixSwitch library82 acting as an intermediary interface between it and high-level physical ideas and algorithms. A connection between the internal SIESTA formats and MatrixSwitch itself has been recently provided, using initially the cubic-scaling libOMM library197 as a test bed, hence still using a dense coefficient matrix, as it corresponds to the case without localization constraints in the solution of the electronic-structure problem. The implementation of a sparse coefficient matrix will make it possible to perform efficient O(N)calculations. The computational effort can be further reduced through the analysis of sparsity of the Hamiltonian and overlap matrices and their re-organization in the block-compressed sparse form. Other developments in the pipeline are linear-response calculations for arbitrary distortions, electronic transport calcu- SIESTA 25 lations with spin-orbit coupling, thermal transport with the Green-Kubo formalism, as described in Ref. 198, a redesign of the molecular dynamics subsystem, and the development of workflows for the generation of data for SCALE-UP. ACKNOWLEDGMENTS SIESTA development has been historically supported by different Spanish National Plan projects: MEC-DGESPB95-0202, MCyT-BFM2000-1312, MEC-BFM2003-03372, FIS2006-12117, FIS2009-12721, FIS2012-37549, FIS201564886-P, and RTC-2016-5681-7, the latter one together with Simune Atomistics Ltd. Currently, we thank financial support from the Spanish Ministry of Science, Innovation and Universities through the grant No. PGC2018-096955-B. We acknowledge the Severo Ochoa Centers of Excellence Program under Grants No. SEV-2015-0496 (ICMAB), and SEV-2017-0706 (ICN2), the GenCat Grant No. 2017SGR1506, and the European Union MaX Center of Excellence (EU-H2020 Grant No. 824143). P.G.-F. acknowledges support from Ramón y Cajal Grant No. RyC-2013-12515. J.I.C acknowledges RTI2018-097895B-C41. R.C. acknowledges to the European Union’s Horizon 2020 research and innovation program under the Marie Skłodoswka–Curie grant agreement no. 665919. D.S.P, P.K, and P.B acknowledge MAT2016-78293-C6, FET-Open No. 863098, and UPV-EHU Grant IT1246-19. V. Yu was supported by a MolSSI fellowship (U.S. NSF award 1547580), and the ELSI development (V.B.,V.Yu) by NSF award 1450280. We also acknowledge Honghui Shang and Xinming Qin for giving us access to the HONPAS code, where a preliminary version of the hybrid functionals support described here was implemented. We are indebted to other contributors to the SIESTA project, whose names can be seen in the file in the Docs/Contributors.txt file of the SIESTA distribution, and we thank those, too many to list, contributing fixes, comments, clarifications, and documentation for the code. The data that support the findings of this study are available from the corresponding author upon reasonable request. 1J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, “The SIESTA method for ab initio order-n materials simulation,” J. Phys.: Condens. Matter 14, 2745–2779 (2002). 2E. Artacho, E. Anglada, O. Diéguez, J. D. Gale, A. García, J. Junquera, R. M. Martin, P. Ordejón, J. M. Pruneda, D. Sánchez-Portal, and J. M. Soler, “The SIESTA method; developments and applicability,” J. Phys.: Condens. Matter 20, 064208 (2008). 3G. Galli, “Linear scaling methods for electronic structure calculations and quantum molecular dynamics simulations,” Curr. Opin. Solid State Mater. Sci. 1, 864 – 874 (1996). 4S. Goedecker, “Linear scaling electronic structure methods,” Rev. Mod. Phys. 71, 1085–1123 (1999). 5P. Ordejón, E. Artacho, and J. M. Soler, “Self-consistent order-ndensityfunctional calculations for very large systems,” Phys. Rev. B 53, R10441– R10444 (1996). 6D. Sánchez-Portal, P. Ordejón, E. Artacho, and J. M. Soler, “Densityfunctional method for very large systems with lcao basis sets,” Int. J. Quantum Chem. 65, 453–461 (1997). 7O. F. Sankey and D. J. Niklewski, “Ab initio multicenter tight-binding model for molecular-dynamics simulations and other applications in covalent systems,” Phys. Rev. B 40, 3979–3995 (1989). 8E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler, “Linear-scaling ab-initio calculations for large and complex systems,” Phys. Status Solidi (b) 215, 809–817 (1999). 9J. Junquera, O. Paz, D. Sánchez-Portal, and E. Artacho, “Numerical atomic orbitals for linear-scaling calculations,” Phys. Rev. B 64, 235111 (2001). 10E. Anglada, J. M. Soler, J. Junquera, and E. Artacho, “Systematic generation of finite-range atomic basis sets for linear-scaling calculations,” Phys. Rev. B 66, 205101 (2002). 11http://www.openmx-square.org/. 12https://www.synopsys.com/silicon/quantumatk.html. 13V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler, “Ab initio molecular simulations with numeric atomcentered orbitals,” Comput. Phys. Commun. 180, 2175 – 2196 (2009). 14A. Carreras, S. Conejeros, A. Camón, A. García, N. Casañ-Pastor, P. Alemany, and E. Canadell, “Charge delocalization, oxidation states, and silver mobility in the mixed silver–copper oxide agcuo2,” Inorg. Chem. 58, 7026–7035 (2019). 15M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, “Density-functional method for nonequilibrium electron transport,” Phys. Rev. B 65, 165401 (2002). 16F. Corsetti, E. Artacho, J. M. Soler, S. S. Alexandre, and M.-V. FernándezSerra, “Room temperature compressibility and diffusivity of liquid water from first principles,” J. Chem. Phys. 139, 194502 (2013). 17X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, T. Deutsch, L. Genovese, P. Ghosez, M. Giantomassi, S. Goedecker, D. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, M. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, D. Sangalli, R. Shaltaf, M. Torrent, M. Verstraete, G. Zerah, and J. Zwanziger, “Abinit: First-principles approach to material and nanosystem properties,” Comput. Phys. Commun. 180, 2582 – 2615 (2009). 18See: https://www.gnu.org/licenses/gpl-3.0.html. 19See: https://launchpad.net/siesta. 20See: https://gitlab.com/siesta-project. 21A. García, M. J. Verstraete, Y. Pouillon, and J. Junquera, “The psml format and library for norm-conserving pseudopotential data curation and interoperability,” Comput. Phys. Commun. 227, 51 – 71 (2018). 22See: https://siesta-project.github.io/psml-docs, accessed November 2019. 23D. R. Hamann, “Optimized norm-conserving vanderbilt pseudopotentials,” Phys. Rev. B 88, 085117 (2013). 24ATOM code for the generation of norm-conserving pseudopotentials. The version maintained by the SIESTA project can be accessed at http: //icmab.es/siesta/Pseudopotentials/index.html. An alternative version is available at http://bohr.inesc-mn.pt/~jlm/pseudo. html. (Accessed July 2017). 25M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, “The PSEUDODOJO: Training and grading a 85 element optimized norm-conserving pseudopotential table,” Comput. Phys. Commun. 226, 39–54 (2018). 26See: http://www.pseudo-dojo.org. 27X. Gonze, F. Jollet, F. Abreu Araujo, D. Adams, B. Amadon, T. Applencourt, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, E. Bousquet, F. Bruneval, D. Caliste, M. Côté, F. Dahm, F. Da Pieve, M. Delaveau, M. Di Gennaro, B. Dorado, C. Espejo, G. Geneste, L. Genovese, A. Gerossier, M. Giantomassi, Y. Gillet, D. Hamann, L. He, G. Jomard, J. Laflamme Janssen, S. Le Roux, A. Levitt, A. Lherbier, F. Liu, I. Lukaˇ cevi´ c, A. Martin, C. Martins, M. Oliveira, S. Poncé, Y. Pouillon, T. Rangel, G.-M. Rignanese, A. Romero, B. Rousseau, O. Rubel, A. Shukri, M. Stankovski, M. Torrent, M. Van Setten, B. Van Troeye, M. Verstraete, D. Waroquiers, J. Wiktor, B. Xu, A. Zhou, and J. Zwanziger, “Recent developments in the ABINIT software package,” Comput. Phys. Commun. 205, 106–131 (2016). 28V. I. Anisimov, J. Zaanen, and O. K. Andersen, “Band theory and mott insulators: Hubbard u instead of stoner i,” Phys. Rev. B 44, 943–954 (1991). 29B. Himmetoglu, A. Floris, S. de Gironcoli, and M. Cococcioni, “Hubbard-