This is the submitted version of the article: Romero A.H., Allan D.C., Amadon B., Antonius G., Applencourt T., Baguet L., Bieder J., Bottin F., Bouchet J., Bousquet E., Bruneval F., Brunin G., Caliste D., Côté M., Denier J., Dreyer C., Ghosez P., Giantomassi M., Gillet Y., Gingras O., Hamann D.R., Ha. ABINIT: Overview and focus on selected capabilities. Journal of Chemical Physics, (2020). 152. 124102: - . 10.1063/1.5144261. Available at: https://dx.doi.org/10.1063/1.5144261 ABINIT: Overview, and focus on selected capabilities Aldo H. Romero,1 Douglas C. Allan,2 Bernard Amadon,3 Gabriel Antonius,4 Thomas Applencourt,3 Lucas Baguet,3, 5 Jordan Bieder,6, 3 François Bottin,3 Johann Bouchet,3 Eric Bousquet,6 Fabien Bruneval,7 Guillaume Brunin,8 Damien Caliste,9 Michel Côté,10 Jules Denier,3 Cyrus Dreyer,11, 12, 13 Philippe Ghosez,6 Matteo Giantomassi,8, 14 Yannick Gillet,8 Olivier Gingras,10 Donald R. Hamann,11, 15 Geoffroy Hautier,8 François Jollet,3 Gérald Jomard,16 Alexandre Martin,3, 6 Henrique P. C. Miranda,8 Francesco Naccarato,8 Guido Petretto,8 Nicholas A. Pike,17, 14 Valentin Planes,3 Sergei Prokhorenko,6 Tonatiuh Rangel,3 Fabio Ricci,6 Gian-Marco Rignanese,8, 14 Miquel Royo,18 Massimiliano Stengel,18, 19 Marc Torrent,3 Michiel J. van Setten,8, 14, 20 Benoit Van Troeye,8, 21 Matthieu J. Verstraete,17, 22, 14 Julia Wiktor,23, 16 Josef W. Zwanziger,24 and Xavier Gonze8, 14, 25 1) Physics and Astronomy Department, West Virginia University, Morgantown, USA. 2) Corning Incorporated, SP-FR-05, Sullivan Park, Corning, NY 14831, USA 3) CEA DAM-DIF, F-91297 Arpajon, France. 4) Dépt. de Chimie, Biochimie et Physique, Institut de recherche sur l’hydrogène, U. du Québec à Trois-Rivières, Trois-Rivières, Canada. 5) IMPMC, UMR 7590 of Sorbonne Université/CNRS/MNHN/IRD, Paris, France. 6) Theoretical Materials Physics/Q-Mat/CESAM, Université de Liège (B5), B-4000 Liège, Belgium. 7) DEN, Service de Recherches de Métallurgie Physique, CEA, Université Paris-Saclay, F-91191 Gif-sur Yvette, France 8) Institute of Condensed Matter and Nanoscience, UCLouvain, B-1348 Louvain-la-Neuve, Belgium. 9) IRIG-MEM, L-SIM , University Grenoble Alpes, CEA , F-38000 Grenoble , France. 10) Dépt. de Physique, U. de Montréal, C.P. 6128, Succursale Centre-Ville, Montréal H3C 3J7, Canada. 11) Dept. of Physics and Astronomy, Rutgers U., Piscataway, NJ 08854-8019, USA. 12) Dept. of Physics and Astronomy, Stony Brook U., Stony Brook, NY 11794-3800, USA. 13) Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA. 14) European Theoretical Spectroscopy Facility, www.etsf.eu. 15) Mat-Sim Research LLC, P.O. Box 742, Murray Hill, NJ, 07974, USA. 16) CEA, DEN, DEC, Cadarache, F-13108 Saint-Paul-Lez-Durance, France. 17) NanoMat/Q-Mat/CESAM, Université de Liège (B5), B-4000 Liège, Belgium. 18) Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain. 19) ICREA - Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona, Spain. 20) Imec, Kapeldreef 75, B-3001 Leuven, Belgium. 21) Dept. of Physics, Applied Physics, and Astronomy, Rensselaer Polytechnic Institute, Troy, New York 12180, USA. 22) Catalan Institute of Nanoscience and Nanotechnology (ICN2), Campus UAB, Bellaterra, 08193 Barcelona, Spain. 23) Dept. of Physics, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden. 24) Department of Chemistry, Dalhousie Univeristy, Halifax, NS, B3H 4R2 Canada. 25) Skolkovo Institute of Science and Technology, Moscow, Russia. (Dated: 6 July 2020) A BINIT is probably the first electronic-structure package to have been released under an open-source license, about twenty years ago. It implements density functional theory (DFT), density-functional perturbation theory (DFPT), many-body perturbation theory (GW approximation and Bethe-Salpether equation), and more specific or advanced formalisms, like dynamical mean-field theory (DMFT) and the "temperature-dependent effective potential" (TDEP) approach for anharmonic effects. Relying on planewaves for the representation of wavefunctions, density and other space-dependent quantities, with pseudopotentials or projector-augmented waves (PAW), it is well suited for the study of periodic materials, although nanostructures and molecules can be treated with the supercell technique. The present article starts with a brief description of the project, a summary of the theories upon which A BINIT relies, and a list of the associated capabilities. It then focuses on selected capabilities that might not be present in the majority of electronic structure packages, either among planewave codes, or in general: treatment of strongly correlated materials using DMFT; materials under finite electric fields; properties at nuclei (electric field gradient, Mössbauer shifts, orbital magnetization); positron annihilation; Raman intensities and electro-optic effect; DFPT calculations of response to strain perturbation (elastic constants, piezoelectricity), spatial dispersion (flexoelectricity), electronic mobility, temperature dependence of the gap, spin-magnetic-field perturbation. The A BINIT DFPT implementation is very general, including systems with van der Waals interaction, or with noncollinear magnetism. Community projects are also described: generation of pseudopotential and PAW data sets, high-throughput calculations (databases of phonon 2 band structure, second-harmonic generation, and GW computations of band gaps), and the library LIBPAW. A BINIT has strong links with many other software projects, that are briefly mentioned. I. INTRODUCTION II. THE PROJECT: HISTORY, COMMUNITY AND IMPACT Since the introduction of density functional theory (DFT) The A BINIT project can be traced back to an initial effort in in 19641 , the field of electronic structure calculations has the late eighties by D.C. Allan, supervised by Michael P. Teter, changed profoundly. This theory became the most popular at Corning Incorporated and Cornell University. The code electronic structure method used to characterize materials at used norm-conserving pseudopotentials and a planewave rep- the atomic scale and has given rise to different research ef- resentation of the Kohn-Sham orbitals, with the local-density forts that have been pushed further by the use and applications approximation (LDA) as the only exchange-correlation func- of DFT. One of the reasons why this theory has become the tional available, and was written in FORTRAN 77. In 1990, X. workhorse of material characterization is the wide distribu- Gonze joined at Cornell, and they implemented density func- tion of computational packages, where solutions of the Kohn- tional perturbation theory (DFPT) on top of what was then Sham equations are implemented. Every package has a differ- dubbed the C ORNING code, in a separate application, named ent scheme and philosophy, but all of them have in common R ESPFN (Corning Incorporated is a large American company the generation of a software that is user-friendly and solves focused on glass and materials applications). After a few years the DFT equations. In that respect, A BINIT 2–5 was one of the in the hands of Biosym inc., where the name was changed to first free licensed solid-state electronic structure packages on P LANE _ WAVE and the code was sold commercially, the de- the market. velopment of C ORNING was stopped in 1996. At that point, Gonze and Allan decided to create a free license code that should be available to the community, with the critical recog- In this publication, we give first a global overview of nition that worldwide collaboration was necessary to develop A BINIT , at the level of its history, community and impact, the code. The previous codes (P LANE _ WAVE , R ESPFN , as well as its main capabilities. Then, we focus on selected and C ORNING ) became the pillars of this new implementa- capabilities, that are perceived as rather specific among the tion. Corning Incorporated agreed to release the source of set of available first-principle packages. P LANE _ WAVE to support this effort and agreed not to enforce a patent they held on the preconditioned conjugate gradient After describing the historical development of the project algorithm.6 The code was rewritten in Fortran 90 with par- and its impact (Sec. II), we present an overview of A BINIT allel features, originally named DFT2000 but soon changed capabilities through a collection of keywords and concepts, to A BINIT (September 1998). with entry points in the documentation (Sec. III). This should The first version of A BINIT was made public in March allow users to know whether A BINIT is able to deliver some 1999, primarily to beta testers. Initially, A BINIT was only specific property of materials. able to find the total energy, electronic charge density, and the electronic structures of periodic systems, using pseudopoten- The two next sections, Sec. IV and V, single out some tials and a planewave basis. In July 1999, the full response specifics of A BINIT . Sec. IV focuses on ground-state and function capability was also available. This implementation electronic properties: correlated materials (Dynamical Mean- was a significant step, as it was one of the only codes that Field Theory - DMFT) Sec. IV A; treatment of finite electric allow users to calculate phonons, dielectric constants, Born field, Sec. IV B; properties at the atomic nuclei, Sec. IV C; and effective charges, etc. The first publicly available release of positron annihilation, Sec. IV D). Sec. V focuses on response A BINIT was made in December 1999, under the GNU Gen- properties: Raman intensities, Sec. IV A; responses to strain, eral Public License. In June 2000, an international advisory Sec. V B; responses including van der Waals interactions, committee was selected, to help with the strategy, support and Sec. V C; spatial dispersion (flexoelectricity and dynamical management of the code. The first A BINIT developer work- quadrupoles), Sec. V D; electron-phonon coupling, Sec. V E; shop took place in Louvain-la-Neuve in 2002. Since then, the temperature-dependent properties of materials in the anhar- developer Abinit meeting takes place every two years. The 9th monic regime, Sec. V F; responses of solids presenting non- developer workshop took place in Louvain-la-Neuve in 2019, collinear magnetism, Sec. V G; response to spin-magnetic per- with the participation of around 60 speakers. turbation, Sec. V H; and temperature-dependent optical spec- The A BINIT spirit is not only to offer the community a free tra, Sec. V I. Finally, we spend a few words on our recent license code but also to encourage users and developers to community efforts related to accurate and efficient Norm- reuse parts of the computational package. The source code is Conserving Pseudopotentials as well as Projector-Augmented always available, and if a user is interested in introducing spe- Waves (PAW) atomic data, Sec. VI A, to High-Throughput cific implementations into the package, an account on GitHub calculations (phonon band structure, second-harmonic gener- suffices, such that changes can be merged and tested in the ation and accurate GW band gap calculations), Sec. VI B, and official version. The A BINIT community has created a series to the LIBPAW library, Sec. VI C. of methods and tests to guarantee the stability of the code, 3 years, where the most recent research efforts are described. With this method, A BINIT tries to acknowledge the devel- opers’ achievements and use these publications to describe how the new implementations are part of the broader pack- age. The A BINIT attitude towards publications is based on research independence, but also collaborative efforts that con- tinue to grow the computational package. The research interest of the developers group (around 50 people at the time of writing) are respected by the A BINIT community and overseen by an international advisory com- mittee. This committee consists of around ten senior scien- tists, which have been involved in the development of A BINIT FIG. 1. A BINIT citation count by adding the four references 2–5 and have been supportive of the philosophy and efforts of the and deleting all record repetitions. See details in text. community. A metric of A BINIT ’s impact can be obtained by using bib- liometric analysis accessed from the Web of Science (data precisely to avoid the introduction of errors by new imple- taken on December 26th , 2019). This can be assessed by mentations. A BINIT is a very well documented code, with adding all citations of all the main papers reporting A BINIT extensive description of installation, details of the input vari- developments, meaning Refs. 2–5. We found a total of 5125 ables and their different dependencies, a very established tuto- citations since 2002, of which 3995 are unique, with an av- rial, with specific applications and a full set of examples (more erage number above 250 total citations per year since 2011 than 800). To support general users, A BINIT has created a Fo- (see Figure 1). Using this database of unique entries, we can rum, which is linked to the A BINIT web site, and allows users address the impact of A BINIT on different fields of science to ask questions related to theory, implementations or use of by counting some of the papers reported in high impact jour- A BINIT . The Forum is a primary resource for users, and both nals. For example, there are 160 papers from Physical Review developers and users answer questions and offer advice. Letters, 109 Journal of Chemical Physics, 108 Computational To spread the use of A BINIT and to train young scien- Materials Science, 32 Nano Letters, 33 Scientific Reports, 18 tists in the field of electronic structure, the A BINIT develop- Nature Communications, 8 Nature Materials, 7 Nature, etc. ers have been supporting schools around the world. While As it is clear from this analysis, A BINIT is a healthy package the first A BINIT workshop was in 2002, there are usually be- that impacts several scientific fields and supports the work of tween 2-3 schools and workshops per year in many different many agencies and institutions, with authors of 84 different places. There is a regular presence at the March Meeting of countries and from more than 2000 different institutions. the American Physical Society, where half day tutorials are or- ganized by the A BINIT community, and at CECAM. Some of the CECAM events are theory oriented workshops, with dis- III. OVERVIEW OF FORMALISMS AND PROPERTIES cussions of the formalism and algorithms behind A BINIT , and examples of calculations are discussed. On the other hand, An overview of the capabilities of a large scientific code can there are also regular CECAM hands-on schools, which al- be structured in different ways. The present work targets users low users to gain experience in the use of A BINIT , and also who simply wish to know whether A BINIT is able to compute learn the theoretical background. These schools are very well a particular property of a material or nanosystem, and possi- attended, but depend on specific locations with sufficient ac- bly to represent it graphically. In this respect, Tables I and II cess to computational facilities. In all cases, the presentations present an alphabetic list of keywords (or concepts), that best of the A BINIT developer workshops, schools and conferences represent a capability of A BINIT , with related documenta- are stored in the A BINIT web page and are accessible glob- tion: possibly some section of the present publication, and/or ally. These presentations offer the users details of code but are reference to some "topics" of the on-line A BINIT Web doc- also an archive of the different implementations and how they umentation https://docs.abinit.org/topics/features, and/or ref- have been tested and used. erence to some scripts and illustrations in the A BI P Y gallery The A BINIT community champions free software devel- http://abinit.github.io/abipy/gallery/index.html. Highlighted opment, but also wants to keep the integrity and modularity (bold) keywords refer to specific sections of the present pub- of the computational package. All new developments need lication. The emphasis on these keywords stems from the be- to be approved by the community, to be sure to keep track lief that such characteristics/capabilities of A BINIT are not of the different efforts, and conduct a minimal screening on commonly found among first-principles codes, and can be put the dependability and ethical behavior of people involved in forward as examples of what makes A BINIT unique. the code development. The driving idea is to preserve a con- In the "topics" tags of the online A BINIT Web documenta- structive and supporting environment for all scientists willing tion, we try to address the generic challenge of software docu- to collaborate. All major implementation efforts in A BINIT mentation, by providing for each topic a hub to the underlying are connected to individual publications but they are also in- theory (including bibliographical references), to the related troduced through specific publications, about every three-five A BINIT input variables, to the example A BINIT input files 4 Keyword/Concept Documentation Atomic mean square displacement (phonons) topic:Temperature, plot_phonons_msqd.html Bader atomic charges topic:Bader Born effective charges Sec. VI B, topic:DFPT Born effective charges, non-collinear magnetism case Sec. V G Correlated electronic state (dynamic - DMFT, static - DFT+U) Sec. IV A, topic:DMFT Conductivity, electrical (electron-phonon coupling) topic:ElPhonTransport Conductivity, thermal (electron-phonon coupling) topic:ElPhonTransport Debye-Waller temperature factors (phonons) topic:Temperature, plot_phonons_msqd.html Dielectric permittivity Sec. VI B, topic:DFPT, topic:Phonons Dielectric permittivity, non-collinear magnetism case Sec. V G Dielectric function (optical - frequency dependent) topic:Susceptibility, plot_mdf.html, plot_scr.html, plot_multiple_mdf.html, plot_optic.html Dielectric function (infrared - frequency dependent) topic:Phonons, plot_phonons_infrared.html Dynamical matrices topic:Phonons Dynamical quadrupoles Sec. V D Effective mass topic:EffectiveMass Elastic tensor Sec. V B, Sec. V C, topic:DFPT, topic:Elastic Elastic tensor, temperature-dependent Sec. V F Electric field (finite, using Berry phase) Sec. IV B, topic:Berry Electric field gradient (at nuclei) Sec. IV C, topic:EFG Electro-optic effect Sec. V A, topic:nonlinear Electric polarization (Berry phase) Sec. IV B, topic:Berry Electron self-energy (GW) topic:SelfEnergy Electron spectral function (GW) topic:GW, plot_gw_spectral_functions.html Electron spectral function (DMFT) Sec. IV A, topic:DMFT Electron-phonon coupling strength topic:ElPhonInt, plot_a2f.html Electronic band gap (GW) topic:GW, plot_qps.html Electronic band gap, temperature-dependent Sec. V E, topic:TDepES Electronic band structure Sec. VI B, topic:ElecBandStructure, topic:GW, topic:ldaminushalf, plot_ebands_edos.html, plot_ebands.html, plot_kpath_from_ibz.html, plot_ejdos.html, plot_qpbands_with_scissors.html Electronic band structure (spin-polarized) topic:ElecBandStructure, plot_ebands_spin.html Electronic Density Of States (DOS) topic:ElecDOS, plot_edos.html, plot_ebands_edos.html Electronic fat band structure topic:ElecDOS, plot_efatbands.html, plot_efatbands_spin.html (with angular momentum weights) plot_ejdos.html Electronic quasiparticles (SCGW) topic:GW, plot_scqpgw.html Eliashberg function (electron-phonon) topic:ElPhonInt, plot_a2f.html Entropy (phonons) topice:Temperature, plot_phthermo.html Excited states (charged excitations) topic:GW, topic:Coulomb Excited states (neutral excitations) topic:DeltaSCF, topic:BSE, topic:TDDFT Fermi surface topic:ElecBandStructure, plot_fermisurface.html Flexoelectricity Sec. V D Free energy (phonons) Sec. VI B, topic:Temperature, plot_phthermo.html Geometry optimization topic:GeoOpt, topic:ForcesStresses, plot_hist.html Grüneisen parameters topic:Temperature, plot_gruneisen.html GW corrections topic:GW, plot_qpbands_with_interpolation.html TABLE I. Overview of A BINIT capabilities (first part, A-G). For each keyword, ordered alphabetically, one or several links are provided. Each highlighted keyword emphasizes a capability of A BINIT that is not present in many first-principles packages, or that is particularly strong (e.g. high-throughput calculations have been performed). A specific section is dedicated to such keywords in the present publication. In addition, other links to documentation are pointed out. “topic:" refers to the A BINIT topic online documentation, available at https://docs.abinit.org/topics/. “plot_" refers to the A BI P Y gallery of plotting scripts, available at http: //abinit.github.io/abipy/gallery/plot_. and to the A BINIT tutorials, also mentioning possible restric- ory (TD-DFT)11 , many-body perturbation theory (MBPT)12 tions in the implementation, or different levels of accuracy of and DMFT13,14 . Each of these formalisms allows one to ad- the implementation. dress different properties at different levels of accuracy. As A BINIT implements different first-principles formalisms, an example, the electronic band gap can be computed from i.e. DFT7 , DFPT3,8–10 , time-dependent density functional the- DFT, but the predicted value is notoriously inaccurate. In- 5 Keyword/Concept Documentation Infrared reflectivity topic:Phonons Interatomic Force Constants topic:Phonons Internal energy (phonons) Sec. VI B, topic:Temperature, plot_phthermo.html Inverse Dielectric function (optical - frequency dependent) topic:Susceptibility, plot_scr_matrix.html Joint Density Of States (electronic) topic:ElecBandStructure, plot_ejdos.html Macroscopic average topic:MacroAve Magnetic field (finite) topic:MagField Magnetic moments topic:MagMom Magnetic susceptibility Sec. V H Molecular dynamics topic:MolecularDynamics, topic:DynamicsMultibinit, topic:PIMD Mössbauer isomer shift Sec. IV C, topic:EFG Optical absorption Sec. V I, topic:BSE, topic:Optic Optical response, temperature-dependent Sec. V I, topic:BSE, topic:Optic Phonon bands Sec. V C, Sec. VI B, topic:Phonons, topic:DFPT, plot_phonons_lo_to.html, plot_phonons.html, plot_ddb_asr.html, topic:Band2eps Phonon bands, temperature-dependent Sec. V F, topic:Tdep Phonon bands, non-collinear magnetic case Sec. V G, topic:Phonons, topic:DFPT Phonon fat bands topic:Phonons, plot_phonon_fatbands.html, plot_phbands_and_dos.html Phonon Density Of States topic:Phonons, plot_phonon_pjdos.html, plot_phbands_and_dos.html Phonon Density Of States, temperature-dependent Sec. V F, topic:Tdep Phonon linewidth (electron-phonon coupling) topic:PhononWidth Piezoelectric tensor Sec. V B, Sec. V C, topic:DFPT, topic:Elastic Positron annihilation Sec. IV D, topic:positron Projected Phonon Density Of States topic:Phonons, plot_phonon_pjdos.html, plot_phonons_msqd.html Quasiparticle energies (GW) topics:GW, plot_qps.html Raman cross section/intensities Sec. V A, topic:nonlinear, topic:Phonons Refraction index Sec. V I, topics:DFPT, topics:Optic, topic:BSE Resistivity (electron-phonon coupling) Sec. V E, topic:ElPhonTransport Scanning Tunneling Microscopy map topic:STM Stopping power of charged particles topic:RandStopPow Second Harmonic Generation Sec. VI B, topic:nonlinear, topic:Optic, plot_optic.html Sound velocity Sec. VI B, topic:PhononBands, plot_speed_of_sound.html Specific heat (phonons) Sec. VI B, topic:Temperature, plot_phthermo.html Structural relaxation topic:GeoOpt, topic:ForcesStresses, plot_hist.html Superconducting transition temperature Sec. V E, topic:ElPhonTransport Thermal Expansion topic:Temperature Thermodynamic properties Sec. VI B, topic:Temperature Thermodynamic properties including anharmonicities Sec. V F, topic:Tdep Transition states, transition paths topic:CrossingBarriers, topic:TransPath Two-phonon DOS, sum and difference spectra topic:PhononBands Unfolding supercell band structure topic:Unfolding, plot_fold2bloch.html Wannier interpolation topic:ElecBandStructure, topic:Wannier, plot_wannier90_abiwan.html Zero-point renormalization of band gap Sec. V E, topic:TDepES TABLE II. Overview of A BINIT capabilities (second part, I-Z). For each keyword, ordered alphabetically, one or several links are provided. Each highlighted keyword emphasizes a capability of A BINIT that is not present in many first-principles packages, or that is particularly strong (e.g. high-throughput calculations have been performed). A specific section is dedicated to such keywords in the present publication. In addition, other links to documentation are pointed out. “topic:" refers to the A BINIT topic online documentation, available at https://docs.abinit.org/topics/. “plot_" refers to the A BI P Y gallery of plotting scripts, available at http: //abinit.github.io/abipy/gallery/plot_. stead, the band gap is better computed using the GW ap- rected functionals. It will not be the purpose of the present proximation within MBPT, also implemented in A BINIT . section to specify for each property the different levels of ac- Moreover, even within MBPT-GW, different levels of ac- curacy. By contrast, such information might be found thanks curacy are available: the simple one-shot G0 W0 approach, to the on-line A BINIT topics. with different plasmon pole models, or more sophisticated Similarly, within some formalism, the A BINIT implemen- self-consistent flavors of GW (SCGW). Similarly, different tation might not be compatible with both the norm-conserving families of exchange-correlations functionals can be used pseudopotential (NCPSP) approach and the PAW approach. in DFT: LDA, generalized-gradient approximation (GGA), Also, the treatment of non-spin-polarized systems, collinear meta-GGA, hybrid functionals, van der Waals (vdW) cor- magnetic systems, collinear antiferromagnetic systems, and 6 non-collinear magnetic systems might not all be available. Variable Description In Ref. [4], we provided such a more detailed description of B Magnetic field the formalism+implementation level for many of the A BINIT c speed of light capabilities (see Tables 1 - 4 of Ref. [4]). A BINIT has evolved C0 Clamped-ion elastic tensor C6,IJ Dispersion coefficient since then, and the best available source of information is now the on-line topics. D/b f q Dynamic matrix e electronic charge Also, A BINIT has several graphical post-processors, in- e0 Clamped-ion piezoelectric tensor cluding A BI P Y and AGATE . A gallery of A BI P Y examples E Total energy (with associated Jupyter notebooks) is available at http:// E Electric Field abinit.github.io/abipy/gallery/index.html. The links to these ε∞ Optical dielectric tensor examples are also mentioned in Tables I and II. Eγ Energy of a nuclear transition e−p Finally, note that some of the capabilities of A BINIT are Ec Electron-positron correlation functional presently not yet available/documented in the latest public re- f band occupancy lease of A BINIT (v8), but will be so in the forthcoming one fmk+q Fermi-Dirac occupation factor (v9). Some of these are mentioned in a section of the present |ii Polarization vector along direction i i, j Reduced coordinates paper, but there is no associated topic or A BI P Y example. n(ω), nqν Bose occupation factor Some of the common variables used within this text are found nα Refractive index along α in Table III. P Polarization re Classic electron radius r linear electro-optic tensor IV. GROUND STATE AND ELECTRONIC PROPERTIES U Onsite Coulomb interaction energy |unk i Bloch wavefunction um Eigen-displacement vector A. Correlated materials: Dynamical-Mean Field Theory U0 Ground-state energy hv2 i Mean square velocity of the nucleus Systems with localized orbitals, such as transition metals, vext External Potential vnk Particle velocity lanthanides or actinides exhibit strong correlation effects and Z Atomic number are difficult to describe using DFT with currently available Z∗ Born effective charge tensor functionals14 . One way to improve their description is to ex- α, β Cartesian coordinates plicitly include the onsite Coulomb interaction U between the αm mode dependent Raman tensor electrons in the correlated orbitals of the system (e.g the d or Γ Mode independent line-width f shells). γ Lattice-strain coupling tensor The DFT+U method, available in A BINIT 15–17 , treats this δ Mössbauer isomer shift interaction statically and is most efficient on magnetic Mott ∆ hr2 i change in size of the nucleus insulators. To describe in a more coherent way systems ηαβ Strain component with various interaction strengths, the method of choice is µ̄ II Frozen ion Flexelectricity DMFT13 . This method solves the local many-body problem τ Positron Lifetime for a given correlated atom in the effective field created by χ (1) First order electric susceptibility the other atoms. This field is self-consistently related to the χ (2) Second order susceptibility (p) solution of the atomic impurity problem. The combination of Φ Interatomic force constants DFT with DMFT14,18 enables the description of realistic sys- |Ψi Planewave basis functions tems with both correlated and non-correlated electrons. The |Φi All electron wavefunctions DFT+DMFT method has been helpful to improve in particular |Φ̃i Pseudized atomic wavefunctions spectral functions (describing both Hubbard bands and quasi- ωl Laser light frequency particle peaks in the spectral function), total energy (for e.g. Ω0 Unit cell volume iron systems, actinides or lanthanides), and magnetic proper- ties (see e.g Ref. 14). TABLE III. Summary of commonly used variables and there mean- ings found within the text. Tensor and vector quantities are shown in The method has been implemented in A BINIT 5,19,21 and bold. is presented in a tutorial. The key points of the available DFT+DMFT implementation are the following: • Correlated orbitals are defined as Projected Local Or- • Impurity solvers directly available in A BINIT are the bitals Wannier functions21 , their localization can easily Hubbard I method and the Continuous Time Quan- be changed by modifying the associated energy win- tum Monte Carlo22 (using either a diagonal or a gen- dows. eral hybridization function) in the simplified but ef- ficient density-density approximation19,23 . Spin–orbit • The method is fully self-consistent with respect to the coupling calculations are possible, using a real valued electronic density (see Fig. 2)19 . imaginary time hybridization function. 7 DFT Define Wannier orbitals Use KS orbitals and energies Diagonalize H KS DMFT Loop KS Hamiltonian H KS [n(r)] Lattice Green’s function Hybridization function Green’s function Impurity Solver (CTQMC) Self-energy Electronic density n(r) FIG. 2. The DFT+DMFT scheme as implemented in A BINIT com- bines the DMFT self-consistent loop which computes the local self-energy, and a DFT self-consistent loop which uses the DMFT Green’s function to compute the DFT+DMFT electronic density. KS stands for "Kohn Sham". See also Refs. 19 and 20. • Internal energy and electronic entropy can be com- FIG. 3. Spectral function of α-cerium computed in LDA+DMFT puted19,23 . along symmetry lines of the Brillouin zone (with U = 6 eV and J = 0.7 eV). Two visible effects brought by LDA+DMFT are the • Impurity and k-resolved spectral functions can be com- upper Hubbard band at 4 eV and the renormalization of the LDA puted using analytical continuation of the Green’s func- band structures near the Fermi level. tion or the self-energy, with an external code (e.g OmegaMaxent24 ). Fig. 3 shows the LDA+DMFT spec- tral function of α-cerium, as computed in A BINIT . derivative is the electrical polarization P, which may be non- zero even in the absence of external fields, e.g. in ferroelectric • Several parallelization schemes can be used17 , allowing materials. Computing this term appears formally equivalent to application of the method to large systems. R computing dr rn(r), that is, the dipole moment density, but this expression is ill-defined in an extended system. This prob- • A BINIT is coupled to the CTHYB impurity solver of lem is solved by the modern theory of polarization35 , in which the TRIQS library25,26 through a C++ interface. This al- the key insight is the recognition that while the polarization is lows one to solve the impurity problem in the fully rota- not well-defined, its derivative with respect to a change in the tionally invariant formulation of the interacting Hamil- system is, and so the polarization may be computed (up to a tonian. constant) through integration. The key formula is • A BINIT is built with a python invocation scheme that fe Z allows any personalised python script to solve the im- P= dk hunk | i∇k |unk i . (1) purity problem. Thus, an experienced user can invoke (2π)3 ∑ n their favourite solver within a DFT+DMFT calculation from A BINIT . Eq. (1) is implemented in A BINIT using both NCPSP and the PAW formalism. In both cases the wavefunction deriva- DFT+DMFT uses effective interactions parameters U and J tives are computed using a finite-difference scheme, so that a as input. They can be computed in A BINIT 27 on the same cor- coherent phase relationship between the k-points may be en- related Wannier orbitals, using the constrained Random Phase sured. From this approach, the polarization may be computed. Approximation (cRPA) method28,29 . In cRPA29 , as the screen- With a scheme in place to compute the polarization, the re- ing arising from electronic transitions among the correlated sponse to a finite electric field may be computed by adding a shell is a correlation effect already described in the CTQMC term −P · E to the total energy36 . Then, its effect is included solution of the impurity problem, only other electronic transi- in the self-consistent energy minimization cycle, through the tions are used to compute the screening29 . The cRPA method gradient δ (−P · E )/δ hunk |37,38 . In this way, Kohn-Sham in A BINIT is documented by a tutorial. states that minimize the total energy including the electric These schemes were recently used in f electron systems to field are found (the polarization from Eq. (1) is updated at compute the f effective interactions parameters20,27,30 and to each step as the Kohn-Sham states evolve), provided that the improve the description of their spectral or structural proper- field is not so strong that the insulating gap breaks down. ties16,20,23,31 and improve our understanding of the supercon- Using this approach we have computed a number of re- ducting symmetries in Sr2 RuO4 32 . sponses, including linear terms that permit validation against DFPT, and nonlinear responses that are only available using finite field calculations38 . B. Finite electric fields A BINIT can be used to compute the response to electric C. Properties at the atomic nuclei fields in several ways. From DFPT, the susceptibility can be calculated33,34 as the second-order derivative of the total en- A BINIT provides several features for computing properties ergy at zero field ∂ 2 E/∂ Ei ∂ EJ . Alternatively, the first order that arise from the electronic structure very near the atomic 8 nuclei. These features are of use for computing and interpret- in reciprocal space and then evaluated at the nuclear site of in- ing a variety of experimental probes, in particular Mössbauer terest by Fourier-transformation; the PAW-sphere corrections spectroscopy39 , and nuclear magnetic and nuclear quadrupole are computed in real space, noting that the use of the PAW resonance spectroscopy40,41 . compensation charge ensures that no inter-sphere contribu- Such observables involve the overlap of the electronic tions need be considered; and the gradient due to the fixed ions wavefunctions with the nucleus, in the case of the Mössbauer is computed using an Ewald-sum method46 . As an example of isomer shift, and the gradient of the electric field at the nu- the application of this feature, we showed41 the sensitivity of cleus, observed in both Mössbauer spectra and magnetic res- the electric field gradient response to strong correlation effects onance experiments. As both observables depend sensitively in LaTiO3 , by combining the field gradient calculation with on details of the electronic structure very close to the nuclear the DFT+U method, also available in A BINIT , see Sec. IV A. position, these are cases where a planewave-only treatment with NCPSP is quite inaccurate and it is necessary to use the PAW formalism. D. Positron annihilation The Mössbauer isomer shift δ , in velocity units, is deter- mined by the overlap of the electronic density with a nucleus One of the unique features of A BINIT is the fully self- undergoing a nuclear state transition, and is given by consistent implementation of two-component density func- tional theory (TCDFT)47,48 within the PAW formalism44 . This c 2πZe2 δ= [ρA (0) − ρS (0)] ∆hr2 i, (2) technique allows one to accurately compute various proper- Eγ 3 ties of annihilating electron-positron pairs, which in turn can be used to interpret positron annihilation spectroscopy (PAS) for the electron density of the absorber (A) and source (S) at measurements49 . the nucleus (here located at 0). In this formula, Eγ , the energy Within TCDFT, the total energy of interacting positrons and of a nuclear transition, and ∆hr2 i, the change in size of the electrons is written as: nucleus, are nuclear properties that must be supplied by the user, while A BINIT can determine the electronic contributions E[n+ , n− ] = E[n+ ] + E[n− ] ρA and ρS . To do this, the PAW formalism is used, in which Z an observable hAi is computed as42–44 + dr vext (r)[n− (r) − n+ (r)] (4) n− (r)n+ (r) Z Z hAi = hΨ| A |Ψi + (3) − dr dr0 + Ece−p [n+ , n− ],  |r − r0 | ∑ hΨ̃| p̃i i h p̃ j |Ψ̃i hΦi | A |Φ j i − hΦ̃i | A |Φ̃ j i . ij where n+ and n− are positron and electron densities, E[n+ ] and E[n− ] are one-component functionals for positron and For properly constructed PAW data, all-electron accuracy may − (r)n+ (r) electrons, and dr dr0 n |r−r R R be recovered in this formalism. For computation of ρ(0) ap- 0| corresponds to the Hartree pearing in Eq. (2), we use as the observable the 3-D delta- interaction. Various approximations can be used to calculate e−p function δ (0)39 . the Ec term, both within LDA48,50 and GGA.51,52 Using the above formalism, isomer shifts may also be com- The TCDFT implementation in A BINIT is based on a uni- puted, given available data for Eγ and ∆ hr2 i for the transi- fied formalism, in which both electron and positron wave- tion in question. To validate the method, a range of shifts functions are expressed using the same mixed PAW basis were computed and then ∆ hr2 i extracted by comparison to (planewaves and atomic orbitals)53 . This means that the elec- known shifts; we found excellent agreement for a variety of tronic and positronic energies and forces can be calculated nuclei, including tin, germanium, and zinc39 . The zinc case self-consistently. It allows, for example, for full geometry op- is a very stringent test because the isomer shift range is so timisation of systems containing positrons, an effect which small and the transition so sharp; accurate modeling requires has been shown to be critical in determining reliable annihila- inclusion of the secondary Doppler shift, which we computed tion features54 . It is worth noting, however, that the positron by first principles from the phonon dispersion relations. The calculations within the PAW method are sensitive to the com- secondary Doppler shift takes the form SD = − hv2 i /2c, and pleteness of the PAW dataset. To achieve accurate description is obtainable once the full phonon dispersion curves are com- of positron densities and wavefunctions, it is often necessary puted. We found that even in this case the PAW-derived ∆ hr2 i to include semicore electrons in the PAW dataset53 . was in excellent agreement with experimental values39 . Various properties of annihilating electron-positron pairs The electric field gradient calculation proceeds along sim- can be calculated within A BINIT . First, positronic wavefunc- ilar lines, except the observable of interest is now the second tions and densities in the direct space can be accessed and derivative of the potential n(r)/|r − r0 |, for density n(r). The visualized to inspect the localization of the particle in a given implementation in A BINIT 40 is similar to that in other codes, system. An example is given in Fig. 4, where the isodensity of such as Q UANTUM E SPRESSO 45 . As usual in PAW, the den- a positron localized in a silicon vacancy in SiC is shown. Sec- sity is decomposed into a planewave-based part, and PAW- ond, total energies of a system containing a positron can be sphere corrections. In addition, there is a charge density due used to calculate binding energies or affinities. Third, based to the ions. The planewave-based field gradient is computed on the electron and positron densities, the positron lifetime τ 9 ∂ ∂ ∂ ∂τ ∂E ∂η ∂ ∂τ Dq Z∗ γ ∂ ∂E Z∗ ε ∞ e0 ∂ ∂η γ e0 C0 TABLE IV. Array of second derivatives of the appropriate energy functional with respect to atomic displacement (τ), electric field (E ) and strain (η) perturbations. In addition to these three types of per- turbations, the magnetic field perturbation is also implemented, with corresponding cross derivatives, see Sec. V H. strain (η) perturbations33,34,63 , which are related to dynam- FIG. 4. Isodensity of a positron localized inside a silicon vacancy in ical matrices at any q-vector (Dq ), optical dielectric tensor silicon carbide. The figure was generated using VMD62 . (ε ∞ ), Born effective charge tensors (Z∗ ), clamped-ion elas- tic tensors (C0 ), clamped-ion piezoelectric tensors (e0 ) and lattice-strain coupling tensor (γ) as summarized in Table IV. can be calculated, as the inverse of the annihilation rate λ : Relying on these bare quantities A BINIT also provides ac- cess to additional properties like the full phonon dispersion 1 Z λ= = πre2 c dr n− (r)n+ (r)g(n− , n+ ), (5) curves (interpolated by separating short-range and long-range τ dipolar interactions), interatomic force constants in real space where g(n− , n+ ) is an enhancement factor, corresponding to (short-range and dipolar contributions), infra-red intensities, the increase in annihilation due to the screening of the positron static and infra-red dielectric constants, relaxed ion elastic and by the electrons. Several parametrizations for g(n− , n+ ) can piezoelectric tensors. All of these are accessible both with be used in A BINIT . Finally, the electron and positron wave- NCPSP and in PAW, using LDA or GGA functionals. The functions can be used to calculate momentum distributions specific case of the strain perturbation63 that is a unique fea- ρ(p) of annihilating pairs, following47 ture of ABINIT is further discussed in Section V.B Going further, A BINIT also implements some third-order Z 2 energy derivatives and responses to additional perturbations, e-p ρ(p) = πre2 c ∑ dr e−ip·r Ψi (r) , (6) providing access to an even much broader set of properties as i discussed below. e-p where Ψi is the two-particle wavefunction in the state i, and p is the given momentum. The TCDFT implementation in A BINIT has already been A. Raman spectroscopy and the electro-optic effect applied to various systems. It has been used to calculate anni- hilation properties of bulk metals and semiconductors53 . Cal- culations for vacancy-type defects in systems such as SiC54–56 Beyond the development of DFPT up to second-order34 , UO2 ,57,58 GaN,59 and yttria-stabilized zirconia,60 have been non-linear properties can be accessed using third-order DFPT. performed and used to analyze experimental data. Positron Thanks to the 2n + 1 theorem9,64,65 , third-order energy deriva- lifetimes and binding energies have also been investigated on tives can be obtained from wavefunction derivatives up to Fe(001) surfaces with and without adatoms61 . first order, as already done in second-order response calcula- tions. So, the additional computational cost to access such third energy derivatives is typically negligible compared to V. RESPONSE PROPERTIES: SPECTROSCOPY, that needed for second-order energy derivatives. There are, VIBRATIONS, DIELECTRIC RESPONSE, TEMPERATURE however, some linear and non-linear properties in which an DEPENDENCE, SPIN-MAGNETIC FIELD COUPLING analytic treatment of the electric field perturbation requires supplementary wavefunction derivatives33 . Many physical properties of materials can be formulated as In A BINIT , the accessible non-linear properties are (i) the derivatives of the energy. Although such derivatives might be second-order optical susceptibilities (χ (2) ), (ii) the Raman computed from finite differences techniques, they can also be tensors (α m ), and (iii) the electro-optic tensor (r). conveniently determined from perturbation theory. A BINIT The second-order nonlinear optical susceptibility is a third- does not only compute many energy derivatives from DFPT rank tensor related to the response of the electrons of the sys- but also gives readily access to various related physical quan- tem to optical electric fields, which are a priori frequency de- tities in Cartesian coordinates and conventional units. pendent. In the present A BINIT context of the 2n + 1 theorem As a basic feature, A BINIT routinely determines standard applied to (static) DFT, we neglect the dispersion of χ (2) and second energy derivatives with respect to atomic displacement compute it as the purely electronic response of the system at (τ q ), homogeneous static electric field (E ) and homogeneous zero frequency. Within these conditions, χ (2) is formulated as 10 a third-order derivative of a field-dependent energy functional: 1.0 0.03 0.8 (2) 1 ∂ 3E 0.02 χi jk = − (7) 0.6 Intensity 2Ω0 ∂ Ei ∂ E j ∂ Ek 0.01 0.4 The χ (2)tensor delivered by A BINIT is related to the non- 0.2 0.00 (2) linear optic “d-tensor” as di jk = (1/2)χi jk and can be calcu- 800 1000 1200 0.0 lated with an additional scissor correction on the electron band 200 400 600 800 1000 1200 energies. Such second-order optical susceptibility has been frequency (cm 1) recently used for high-throughput Second-Harmonic Genera- tion calculations, see Sec. VI B 0 b. FIG. 5. Raman spectra of polycrystalline quartz, α-SiO2 . The green The Raman tensor describes the change of linear optical line is obtained from DFPT and PAW pseudopotentials, where the susceptibility, χi j , produced by an atomic displacement τκβ of peak width is arbitrary set to 3 cm−1 . The black dotted line shows atom κ in direction β . It is therefore related to the following experimental data from Ref. 73. third energy derivative66 : (1) discussed above) was first implemented in A BINIT relying on ∂ χi j 1 ∂ 3E =− (8) the "Berry phase" formalism66 ; this implementation is avail- ∂ τκβ Ω0 ∂ τκβ ∂ Ei ∂ E j able using a NCPSP approach within the LDA and for unpo- larized systems. Following the scheme of Gonze33 , Miwa72 Using Placzek’s approximation67 , the Raman intensity αimj of used an alternative approach, where the electric field pertur- a Stokes process associated with a phonon mode m depends bation is treated analytically, leading to a derivative operator on the Raman tensors as66,68,69 : with respect to k vectors (∂ /∂kα ). At the non-linear level, (1) this derivative operator is applied to first-order derivatives of ∂ χi j wavefunctions, leading to second-order wavefunction deriva- αimj =∑ um (κβ ). (9) κ,β ∂ τκβ tives. This way, even in the context of the 2n + 1 theorem, a (non self-consistent) second-order Sternheimer equation has Depending on the type of phonon mode, i.e. transverse optic to be solved. Our second implementation follows that scheme, (TO) or longitudinal optic (LO), the derivative of the electric and is valid for both NCPSP and PAW approaches. It is avail- susceptibility with respect to an atomic displacement must be able with LDA functionals, and is valid for unpolarized or written as66,70 : collinear spin-polarized systems. This second implementa-   tion also converges faster than the original one in terms of the (1) (1) ∂ χi j ∂ χi j k-point sampling. The extension to the PAW+U formalism is  = (10) under development.  ∂ τκβ ∂ τκβ E =0 A BINIT delivers all the tensors given above in Cartesian T O ∗ q coordinates and with clear units (SI or conventional literature  (1) (1) ∑ Zκβ ,l l ∂ χi j  = ∂ χi j 8π l (2) choices). These quantities can then be combined to compute − χi jl ql (11) ∂ τκβ E =0 Ω0 ∑ ql εll∞0 q0l ∑  ∂ τκβ l the polarization dependent Raman spectra of single crystals LO l,l 0 or even Raman spectra of powders. In either case, the Raman where q is the direction along which the Brillouin zone center intensities can be obtained as74 : is approached. Both quantities are provided by A BINIT . Iimj (ω) = Cm (ω)| hi| α m | ji |2 , The electro-optic tensor is related to the change of linear optical susceptibility produced by a static field. As such, I m (ω) = 2πCm (ω) ((10Gm m m m 0 + 4G2 ) + (5G1 + 3G2 )) , (13) it combines a priori electronic, lattice and strain contribu- tions. The so-called clamped (fixed strain) electro-optic tensor where G0 , G1 , and G2 are functions of the components of combining electronic and lattice contributions is directly pro- Eq. (9)75 and a Lorentzian function is used for the broaden- vided by A BINIT from the knowledge of previously discussed ing: quantities66 : (ωm − ωl )4 Γ Cm (ω) = [n(ωm ) + 1] , (14) 2ωm c4 (ω − ωm )2 + Γ2 " # −4π (2) 1 αimj ∗ ri jγ = 2 2 2χi jγ + √ ∑ 2 ∑ Zκβ ,γ um (κβ ) . (12) ni n j Ω0 m ωm κ,β Post-processing scripts are provided with A BINIT to extract the Raman tensor and print the powder-average and polariza- The unclamped electro-optic tensor includes an additional tion dependent Raman intensity to a file. An example of a piezoelectric contribution (strain response) that can be com- powder-averaged Raman spectra for polycrystalline quartz is puted but requires an independent calculation of the elasto- shown in Fig. 5. optic coefficients by finite differences71 . Using the non-linear tensors computed by A BINIT , theo- The computation of the third-order energy derivatives in retical Raman spectra have been successfully compared with Eq. (7) and Eq. (8) (providing access to the different tensors experiment for many materials including oxides76–78 , organic 11 compounds79 and transition metal dichalcogenides80,81 (see terms of the real- and reciprocal-space primitive lattice vectors also the WURM database of Raman spectra82 ). Electro-optic as tensors have also been computed for various compounds in- cluding ferroelectric oxides71,78 , multiferroics83 and mono- Ξi j = ∑ RPαi RPα j , ϒi j = ∑ GPαi GPα j . (17) α α layer transition metal dichalcogenides81 . Going further, re- lying on these quantities obtained in DFT at zero Kelvin, non- Dot products of real- and reciprocal-space vectors can be ex- linear optical properties at finite temperature were also suc- pressed in terms of their strain-invariant reduced-space coun- cessfully computed making use of a first-principles based ef- terparts using these tensors. For example, reciprocal-space fective Hamiltonian approach84 . vectors K = k + G represented by their reduced counterpart K̃ have dot products B. The strain perturbation in density-functional perturbation K0 · K = ∑ K̃i0 ϒi j K̃ j . (18) theory ij First and second strain derivatives of the metric tensors follow One of the unique capabilities of A BINIT is the direct cal- from the derivatives of real and reciprocal space vectors with culation of the elastic and piezoelectric properties of materi- respect to strains als by DFPT. This involves calculating the responses to three types of perturbations: atomic displacements, electric fields, ∂ Xγ ∂ Kγ = δαγ Xβ , = −δαγ Kβ . (19) and strains and appropriately combining the results85 . The ∂ ηαβ ∂ ηαβ general structure of DFPT is based upon the systematic expan- The only simplifying aspect of all this is that the dot product sion of the variational expression for the DFT total energy in of a real and a reciprocal lattice vector, typically appearing in powers of a parameter λ characterizing some dependence of phase factors, are merely 2π times their reduced counterparts, the energy functional. Such parameters as internal atomic co- and therefore strain-independent63 . ordinates and the macroscopic electric field could be handled All the terms in the total energy expression in reduced coor- in this framework in a conceptually straightforward manner. dinates can be expressed in terms of dot products. Therefore Treating macroscopic strain as a parameter within this formal- the evaluation of Eq. (15) is in principle straightforward while ism, however, is less straightforward. The approach developed in practice exceedingly tedious. This formalism was initially for A BINIT was based on its existing overall formulation of developed for multi-projector NCPSP and LDA exchange- the DFT energy expression in reduced coordinates. This in- correlation functionals. There are special considerations for troduces real- and reciprocal-space metric tensors into every GGA functionals which were added87 . Extension to PAW po- expression. The underlying reduced-coordinate lattice struc- tentials within LDA was recently added88 . ture is invariant – all unit cells are unit cubes. Strain affects The terms described so far, derivatives with respect to two only the metric tensors, and the λ derivatives of every term strain components (the elastic tensor), one strain component in the Hamiltonian can be developed from derivatives of these and a uniform electric field (the piezoelectric tensor) and one tensors63 . strain component and one atomic displacement (the so-called The second derivatives of the energy with respect to a pair internal strain) are all based on strain-independent atomic po- of perturbations are evaluated using so-called non-stationary sitions in reduced coordinates. In fact, infinitesimal strains do expressions produce infinitesimal reduced-coordinate changes in atomic occ D (λ ) (λ ) (λ1 ) (0) E positions. Taking this into account modifies the clamped-atom Eelλ1 λ2 = ∑ ψα 2 (T (λ1 ) +Vext 1 + HHxc0 ) ψα elastic and piezoelectric tensors calculated so far. occ D (0) α E 2 , Combining the clamped-ion derivatives with the atomic (λ λ ) (0) + ∑ ψα (T (λ1 λ2 ) +Vext 1 2 ) ψα + 12 ∂∂ λE∂Hxc λ (0) position contributions is reminiscent of the chain rule for α 1 2 n (15) derivatives85 . The necessary intermediate derivatives (inter- where (λ ) ψα 2 is first-order wavefunction obtained from the atomic force constants and Born effective charges) were al- self-consistent Sternheimer equation for perturbation λ2 , the ready available in A BINIT . The DFPT calculations of the parenthesized superscripts denote partial derivatives, and the clamped-ion tensors and these additional derivatives are per- notation for the various energy terms and the ground-state formed separately in A BINIT , and the stored results combined wavefunctions should be self-explanatory33 . The derivative in the auxiliary program A NADDB to produce the final, phys- with respect to strain and the electric field, necessary for the ical relaxed-atom quantities. Both the NCPSP and PAW for- piezoelectric tensor, has the alternative form malism and codes have been verified to a high degree of accu- racy, by comparison with finite-difference calculations. Z occ D ∂ 2 Eel Ω (k̃ ) (η ) E =2 ∑ i ψkmj ψkmαβ dk , (16) ∂ E j ∂ ηαβ ˜ (2π)3 m C. Density-functional perturbation theory including van der BZ Waals interaction: phonons and strain perturbation where ψkm (k̃ j ) is the first-order wavefunction for the ∂ /∂ k perturbation86 . The metric tensors involved in this context are Most DFT functionals lack a proper treatment of the long- simple and constant throughout space. They are expressed in range e− - e− correlation. This term is crucial to properly 12 describe the structure of weakly-bound compounds and, con- lighted the importance of the dispersion contributions to in- sequently, their underlying properties (phonons, elastic con- teratomic force constants and to phonons80,95 . Such meth- stants, ...). Several methods have been developed in the past ods generally lead to a better agreement with the experiments in order to include some long-range e− - e− correlation in DFT than ordinary GGA-PBE. Similarly to the general DFPT im- computations. Notably, one could mention Grimme’s DFT-D plementation, not only the zone-center phonons are accessi- methods89,90 , the vdW-DF functionals91,92 as well as the TS- ble, but also those associated to an arbitrary wavevector, al- vdW and MBD methods93,94 . All of these approaches add lowing the computations of phonon band structures in a fully- new contributions to the energy and its derivatives, that have consistent manner. In addition, the present implementation to be properly taken into account in DFPT computations for can also be straightforwardly used to study thermal expansion self-consistency. within the quasi-harmonic approximation (see Ref. 97 for an A BINIT allows one to carry out such computations for the example application) or electron-phonon coupling. Grimme’s DFT-D2, -D3 and -D3(BJ) methods95,96 . These methods add to the DFT energy a contribution that depends only on the atomic positions: D. Spatial dispersion: flexoelectricity and dynamical quadrupoles Edisp = ∑ C6,IJ ({R}) f (RIJ ), (20) I,J The flexoelectric (FxE) effect, where a strain gradient de- formation results in a macroscopic polarization, is a challeng- where I and J are atomic indices, C6,IJ is the dispersion coeffi- ing electromechanical coupling to simulate computationally. cient that depends on the whole set of atomic positions R, and The main reason is the a priori incompatibility of a spatially f (RIJ ) is a function that depends on the chosen Grimme fla- varying strain perturbation with periodic boundary conditions. vor (∼ R−6 IJ at long range) and on RIJ , the difference in atomic A recent set of developments in A BINIT has overcome this position between the I and J atoms. This dispersion contri- difficulty by adapting the long-wave method introduced by bution is not compatible with all ad hoc exchange-correlation Born and Huang98 in the early days of quantum mechanics functionals (see Refs. 89 and 90 for the current list of compat- to the modern tools of the DFPT. A BINIT now has the unique ible functionals). The common GGA-based functionals, like capability to calculate the bulk FxE tensor for any insulating GGA-PBE or revPBE, are nonetheless available. Due to the material. way these DFT-D methods were constructed, they may not be Achieving this result has taken nearly a decade of contin- the most suited to deal with metals. uing efforts, in order to settle the remarkable number of for- Nearly all response functions available in A BINIT can be mal subtleties. Based on the seminal work of Resta,99 early computed in DFPT with the inclusion of DFT-D dispersion attempts100,101 made use of supercell geometries to calculate corrections, with the only exception of flexoelectricity. Note the relevant ingredients discussed below Later, a long-wave that, since these corrections do not rely on the knowledge of framework in reciprocal space102 was established, together the electronic density, the (total) energy derivatives with re- with a curvilinear-coordinates formulation,103 which led to spect to electric fields remain untouched by the dispersion the calculation of the full flexoelectric response of SrTiO3 scheme. This is both an advantage and a drawback of the slabs104 . (See Ref. 105 for a summary of the pre-2015 de- DFT-D method, since it eases considerably its implementa- velopments.) En route towards a practical implementation, tion for response functions, but it prevents at the same time additional technical and formal issues were addressed, regard- the method to capture any direct electronic contribution (di- ing e.g. the generalization of the uniform strain response to electric constants, Raman susceptibility, etc.). The lack of a "metric-wave" perturbation106,107 , and the proper definition literature on this topic does not allow us to rule on the im- of the current-density operator in presence of nonlocal pseu- portance of dispersion contributions in the proper description dopotentials108 . These efforts culminated with the present of these properties. The reason why the DFT-D dispersion long-wave method,109 which paves the way towards the com- corrections are not available for flexoelectricity comes from putation of many spatial dispersion properties with a com- the indirect dependence of this property on the spatial deriva- putational burden that is comparable to conventional linear- tive of the dynamical matrix, as mentioned in Sec. V D. Such response calculations. spatial derivative is not yet implemented. An ab initio calculation the FxE tensor (µ) requires the All response functions related to the energy derivatives with sum of an electronic (clamped-ion), a lattice-mediated and a respect to atomic displacements or to strains are directly af- mixed contribution. This intricate structure can be formulated fected by the dispersion corrections. Thanks to their pair-wise as102,105 form, the DFT-D strain-related quantities can be related to derivatives with respect to atomic positions II II (1,λ ) κ 1 ∗ κ µαλ ,β γ = µ̄αλ ,β γ − Pα,κρ Γρβ γ + Zκ,αρ Lρλ ,β γ , (22) | {z } | {z } |Ω0 {z ∂ ∂ elec. mix. } = ∑ Rκβ , (21) latt. ∂ ηαβ κ ∂ Rκβ where κ labels the atomic sublattice and the II superscript in- where κ is an atomic index. dicates that the strain gradient tensor is defined as the gradient The investigation of several weakly-bound systems (ben- of the symmetric strain, i.e., the so-called type-II definition is assumed. The quantities on the right-hand side are: µ̄αλ II zene, MoS2 , etc.) using DFT-D methods has already high- ,β γ 13 (1,λ ) the clamped-ion FxE tensor, Pα,κρ the first spatial moment of faster convergence of the phonon bands calculation with re- the polarization field induced by an atomic displacement (also spect to the density of q points and, in some materials, repre- known as P(1) tensor), Γκρβ γ the piezoelectric internal strain sents the only route to obtain the correct sound velocities. tensor, and Lρλ κ The A BINIT implementation of the long-wave DFPT has ,β γ the flexoelectric internal strain tensor. ∗ been used to compute the clamped-ion FxE tensor of Si and Whereas Z and Γ are well known quantities which can be SrTiO3 obtaining an excellent agreement with existing results obtained by means of the linear response DFPT with standard in the literature109 . The computed spatial dispersion proper- perturbations such as atomic displacements, electric field and ties accurately reproduce the sum rules predicted in Refs. 110 strain33,34,63 , Eq. (22) demands the calculation of up to four II (1,λ ) and 102 that relate them with specific quantities from uniform new tensor quantities. These include the µ̄αλ ,β γ and Pα,κρ perturbation theories. tensors introduced above, but also the first spatial moments of Currently, the implementation is restricted to the use of (1,λ ) the dynamical matrix (Φκα,κ 0 ρ ) and of the piezoelectric force- NCPSP without non-linear core corrections, and the LDA response tensor (C̄αλ κ κ functional. ,β γ ), both required to build Lρλ ,β γ fol- lowing the prescriptions in Ref. 102. These new tensors repre- sent four spatial dispersion properties, i.e., they are only sen- E. Electron-phonon interaction sitive to a perturbation gradient, not to a uniform perturbation. In principle, however, it is possible to obtain them as momen- tum derivatives of counterpart tensors from uniform perturba- With A BINIT , it is possible to compute many physi- tions, when the corresponding DFPT second order energies cal properties related to electron-phonon (e-ph) interaction. are generalized to finite momentum q. Thus, for example, the In metallic systems, for instance, one can study conven- clamped-ion FxE tensor is the spatial dispersion counterpart tional superconducting properties within the isotropic Migdal- of the clamped-ion piezoelectric tensor (ēα,β γ )63 and could be Eliashberg formalism111 and compute transport properties in calculated as the momentum derivative of a finite-q version of the normal state by solving the linearized Boltzmann equa- ēα,β γ . tion within the LOVA approach112–114 . In this section, how- This is precisely the procedure adopted in the long-wave ever, we focus on temperature-dependent band structures and DFPT of Ref. 109 that has been implemented in A BINIT . In the zero-point renormalization of the band gap in semiconduc- order to obtain the finite-q generalization of the second order tors115,116 , a subject that in recent years has received increased energies, the approach reformulates the electric-field and the attention within the electronic structure community117–120 . As strain perturbation problems as the time derivative of a vector discussed in the review paper by Giustino111 , the renormal- potential and as the gradient of a metric wave,107 respectively. ization of the electron state due to the e-ph interaction is de- In this way, both perturbations are generalized to finite q, as scribed by the self-energy Σe-ph (ω) = ΣFM (ω) + ΣDW . The is already the case for atomic displacements. This enables us diagonal matrix elements of the Fan-Migdal self-energy in the to carry out an analytical third order derivative of the energy KS basis set are defined by with respect to two of the standard perturbations, and to the Z dq momentum q, which directly provides the sought-after spatial ΣFM nk (ω) = ∑ |gmnν (k, q)|2 × (24) dispersion tensors. Remarkably, by virtue of the 2n + 1 the- m,ν Ω BZ BZ orem9 , the third-order energies are computed in one shot us-   nqν + fmk+q nqν + 1 − fmk+q ing precalculated first-order response functions to the standard + , ω − εmk+q + ωqν + iη ω − εmk+q − ωqν + iη perturbations, without the necessity of self-consistently com- puting any response function to a perturbation gradient. After where η is a positive real infinitesimal. The e-ph matrix ele- execution, the long-wave DFPT routines generate a derivative ments gmnν (k, q) are given by database that is subsequently used by post-processing tools implemented in A NADDB to compute and print the different gmnν (k, q) = hψmk+q |∆qν V KS |ψnk i , (25) contributions to the FxE tensor. The implementation also provides access to another spatial with ∆qν V KS the first-order variation of the self-consistent KS dispersion property, the dynamical quadrupoles, which can be potential that can be computed with DFPT33,121 . The static obtained as the symmetrized sum of the P(1) tensor,109 Debye-Waller (DW) term involves the second order derivative   of the KS potential with respect to the nuclear displacements. (2,γδ ) (1,γ) (1,δ ) Qκβ = Ω0 Pδ ,κβ + Pγ,κβ . (23) State-of-the-art implementations approximate the DW contri- bution with The dynamical quadrupoles are the spatial dispersion coun- terparts of Z∗ and can be used in lattice dynamics calcula- g2,DW mnν (k, q) ΣDW nk = ∑ (2nqν + 1) , (26) tions to improve the prevalent dipole-dipole treatment of the qνm εnk − εmk long-range interactions. The A NADDB routines that carry 2,DW out the process of interpolating the dynamical matrix fol- where gmnν (k, q) is an effective matrix element that, within lowing Ref. 34 have been adapted to incorporate the dipole- the rigid-ion approximation, can be expressed in terms of the quadrupole and quadrupole-quadrupole electrostatic interac- gmnν (k, q) matrix elements111 . In principle, the quasi-particle tions derived in Ref. 102. This new functionality results in a (QP) excitations are defined by the solution(s) in the complex 14 e-ph plane of the equation z = εnk + Σnk (z). In practice, the prob- For the computation of electron lifetimes, A BINIT imple- lem is usually simplified by seeking approximated solutions ments advanced integration techniques that take advantage of along the real axis following two different approaches. In the the linear tetrahedron method44 and double-grid techniques on-the-mass-shell approximation, the QP energy is given by for the q-space integration. Last but not least, A BINIT pro- the real part of the self-energy evaluated at the bare KS eigen- vides a specialized driver to interpolate the DFPT scattering value potentials, compute e-ph matrix elements on arbitrarily dense QP e-ph k-meshes and save the results to NETCDF files. This driver can εnk = εnk + ℜ Σnk (εnk ). (27) be used by externals codes such as B ERKELEY GW 132 to treat The second approach linearizes the self-energy near the KS e-ph interactions at the GW level as discussed in Ref. 133. eigenvalue and evaluates the QP correction using QP e-ph εnk = εnk + Znk ℜΣnk (εnk ) (28) F. A-TDEP: Temperature dependent thermodynamic properties using the TDEP approach with the renormalization factor Z given by −1 The features and the temperature dependence of the phonon  " e-ph # ∂ Σnk spectrum determine a large number of thermodynamic proper- Znk = 1 − ℜ  . (29) ∂ε ties of crystals, such as dynamic and thermodynamic stability, ε=εnk elastic properties, and heat transport within the material. The Both approaches are implemented in A BINIT although it theory of lattice dynamics at finite temperature has long been should be noted that, according to recent works, the on-the- a central field of research in condensed matter, and its appli- mass-shell approach provides results that are closer to those cations are numerous in geophysics, material science, astro- obtained with more advanced techniques based on the cu- physics... mulant expansion122 or self-energy calculations employing an The calculation of the phonon spectrum is available in eigenvalue–self-consistent cycle123 . A BINIT since the 90’s and is based on DFPT. In this frame- Accurate calculations of e-ph renormalization are still chal- work, the ground state undergoes a small perturbation around lenging even by present standards, because the e-ph self- its 0 K equilibrium positions in order to describe the shape energy is quite sensitive to the q-point sampling. Moreover of the potential energy surface (PES) very near the minimum. a large number of empty states m is usually required to con- Thereafter, in order to capture the temperature effects, the har- verge the real part of the self-energy. In order to avoid the ex- monic approximation (HA) or quasi-harmonic approximation plicit computation of the DFPT scattering potentials on dense (QHA) can be applied. The latter approximation is handled by q-grids, A BINIT generalizes the Fourier-based interpolation assuming that the temperature dependency of the phonon fre- scheme proposed by Eiguren et al.124 to account for the non- quencies can be taken into account implicitly through a vari- analytical behaviour associated to the long-range interactions ation of the volume. This treatment generally gives excellent present in semiconductors125 . The dipolar potentials gen- agreement with experimental measurements of the thermal ex- erated by the Born effective charges in polar materials are pansion and the elastic constants134 . treated using a generalized Fröhlich model126,127 . Further de- However, the QHA fails or cannot be applied in various par- tails concerning the implementation are given in Ref. 17. In ticular cases. For example, when the system is close to a phase order to accelerate the convergence with the number of empty transition, when the phase of interest is not stable at 0 K, or states, A BINIT replaces the contributions given by the high- when one wants properties at high temperature. In these situ- energy states above a certain band index M with the solu- ations, an implicit treatment of the temperature through a vol- tion of a non-self-consistent Sternheimer equation in which ume variation is no longer sufficient and a dedicated calcula- only the first M states are required. The methodology, pro- tion, taking into account the temperature explicitly, is needed. posed in Ref. 128, is based on a quasi-static approximation in One can formalize and synthesize this statement as follows: which the phonon frequencies in the denominator of Eq. (24)         are neglected and the frequency dependence of Σ is approx- ∂ω ∂ω ∂ω ∂V = + . (31) imated with the value computed at ω = εnk . This approxi- ∂T p ∂T V ∂V T ∂ T p mation is justified when the bands above M are sufficiently high in energy with respect to the nk states that must be cor- The second term of the right hand side of Eq. (31) is in- rected. Furthermore, this upper-bands contribution to Σ con- cluded in the QHA. However, the first term on the right side verges quickly with respect to q-points sampling, and it can of Eq. (31) is only included when the temperature dependence be safely computed on a coarse q-grid123,129 . is explicit. The code can compute QP corrections and lifetimes due to Since the beginning of the 60’s, a large number of theo- e-ph scattering as well as spectral functions. The lifetimes retical studies have been carried out in order to go beyond obtained from the imaginary part of Eq. (24) can be used to the simple harmonic crystal and to deal with anharmonic ef- compute carrier mobilities within the self-energy relaxation fects135–141 . More recently, there has been a revival of such time approximation111,130,131 : developments. During the last decade, several groups pro- posed theoretical frameworks coupled to computational meth- −1 dk ∂f Z µe,αβ = vnk,α vnk,β τnk , (30) ods able to capture anharmonic effects in crystals142–148 . One Ω0 ne ∑ n ΩBZ ∂ε εnk way, chosen and implemented in A BINIT , has been proposed 15 by Esfarjani and Stokes149 then developed by Hellman et Once the IFC coefficients are obtained, many dynamic, al.150–152 . In the framework of TDEP, the anharmonic terms elastic and thermodynamic properties can be evaluated: are treated in an effective manner. Let us define a 3-dimensional crystal and consider that a • As usual, the phonon spectrum is obtained after the di- ground state energy U0 is obtained when the atoms are in their agonalization of the dynamical matrix. equilibrium positions. The potential energy of this system can be rewritten using a Taylor expansion around equilibrium po- • Then, the phonon density of states is built and used to sitions as compute various thermodynamic properties153 : the free energy, the specific heat, the entropy... (p) p 1 α1 ...α p αk U = U0 + ∑ i1 ...i p ∏ uik . (32) p ! α1∑ Φ • By computing the elastic constants using the second- p≥1 ...α p k=1 i1 ...i p order effective IFC135,141 , it is possible to obtain the bulk and shear moduli. (p) α ...α with Φ i11...i p p the interatomic force constants (IFC) at the pth • The Grüneisen parameter is obtained using the third or- order and uαi the displacement of atom i along Cartesian di- der effective IFC154,155 which leads to other important rection α. quantities: the thermal expansion, the isentropic com- By performing ab initio molecular dynamics (MD) simula- pressibility, the constant pressure specific heat, sound α (t) and atomic displacements uα tions, the forces Fi,MD i,MD (t) velocities... are acquired at each time step t. If we differentiate Eq. (32) with respect to atomic displacements and insert the MD quan- Since the effective IFC include an explicit temperature depen- tities, we obtain the following system of equations : dency, all the quantities listed above do as well. p The A BINIT implementation of the TDEP algorithm, re- P (p) α1 1 α α2 ...α p αk ferred to as A - TDEP , has produced a number of applications FMD,i (t) = − ∑ ∑ Θ i11i2 ...i ∏ uMD,i (t), 1 p≥2 (p − 1) ! α2 ...α p p k=2 k in the last five years156–160 . We recommend the user to read i2 ...i p the A - TDEP documentation distributed in the A BINIT package (33) before any calculations, and also to read the article dedicated to this implementation161 . (p) αβ ...δ with Θ i j...l the effective IFC at the pth order and P the max- imum order of the expansion. Eq. (33) is non-linear with re- G. Density functional perturbation theory within spect to the atomic displacements and can be rewritten as: non-collinear magnetism α α Fi,MD (t) = ∑ fi,pλ (uMD (t))θ pλ , (34) pλ The necessity to treat non-collinear magnetism (NCM) in DFT codes is more and more important for spintronic appli- where θ pλ being the λ th are now the effective IFC at the pth cations studies. NCM can arise naturally, from geometrical α (u order and fi,pλ MD (t)) a function gathering all the contri- frustration of anti-ferromagnetic interactions, e.g. in triangu- butions coming from the atomic displacements. This sys- lar lattices, by the magnetic anisotropy generated by the pres- tem of equations is now linear as a function of θ pλ (the un- ence of a preferred direction of magnetization, or by the com- known variables of the system) and can be solved using a petition between exchange interactions. However, the gen- least-squares method. The solution is given by Θ = f† .FMD , eralization of DFPT for non-collinear magnetic simulations with FMD ≡ Fi,MDα (t), Θ ≡ θ † pλ and f the pseudoinverse of has only recently been attempted. Within the non-collinear f ≡ fi,pλ (uMD (t)). α DFT scheme, the density is described by a four-component (p) (p) matrix instead of two (in collinear magnetism) which is re- Note the difference between the true Φ and effective Θ quired to account for the spin-orbit coupling (SOC), except IFC. The latter includes the effects of all the IFC not taken when time-reversal symmetry induces vanishing of the mag- into account in the expansion (i.e. beyond the Pth order). For netization everywhere162 . This last term, in current DFT code instance, if P = 2, the 2nd order effective IFC incorporates all implementations, implies a partial loss of symmetry163–165 , higher anharmonic terms, in an effective way. Consequently, which worsens the calculation time and explains the scarcity the effective IFC acquire a dependency on the temperature, at of such calculations. Additionally, the treatment of the elec- odds with the true IFC. tronic exchange-correlation (XC) term is more complicated, The number of time steps needed to solve Eq. (34) has to since the functionals are usually built in a collinear frame- be large with respect to the number of unknown variables θ pλ . work. We have proposed166 three different approaches to per- By taking into account properties of the system (translation form the energy derivatives in the non-collinear regime treat- and rotation invariances, crystal symmetries...) the number ing opportunely the XC term, and have implemented these for of independent and non-zero IFC coefficients of a cubic sys- the atomic displacement (phonons) and electric field perturba- tem with one hundred atoms in the supercell can be reduced tions. to around ten at the second order, tens at the third order and The first method exploits the local transformation from around a hundred at the fourth order. non-collinear to collinear magnetization: at each point of 16 space, we align the local system of reference to the local quan- response calculations within the DFPT formalism has been es- tization axis of the magnetization, then we perform the deriva- tablished several decades ago168 , practical implementations of tives in this locally collinear framework and finally we return this method are still rare in DFT codes169 . to the lab frame of reference. The unitary transformation di- We have implemented and tested the variational DFPT agonalizing the 2 × 2 density matrix ρ̂ (0) , and its derivative method for computing the linear spin response to a small, q- ρ̂ (1) , have to be determined at the zeroth (U (0) ) and first order dependent, magnetic field of the form B = beiq·r + b∗ e−iq·r . (U (1) ) of perturbation, to fully characterize the problem (see The first derivatives of the density and magnetization are ob- Ref. 166). (1) tained by considering the first order external potential vext = The second method consists in writing the U (0) and U (1) 1 2 σα , where the Cartesian component α of the Pauli matrix analytically in terms of the Pauli matrices σαβ , i. e., the gen- corresponds to the direction of the external field. The pe- erators of the Lie group SU(2). riodic part of the first order (spin) magnetization is repre- The third method evaluates explicitly the expression of the (1) (1) † (1) sented as mq = ∑nk u†nk σ unk+q + unk−q σ unk , with the sum first order XC potential, specifically for the LSDA expression running over the occupied Kohn-Sham orbitals n. For an of the XC energy, which is easy to derive. arbitrary wavevector q 6= 0 it is thus necessary to find the The three methods have been tested on two particular sys- (1) tems Cr2 O3 and RuCl3 166 . Cr2 O3 has been used to test Bloch-periodic parts of the first-order wavefunctions unk−q (1) the implementations because it is one of the simplest and and unk+q at both k + q and k − q points. This is done by best known collinear antiferromagnetic and magneto-electric performing two state-by-state minimizations of the second or- systems. Thanks to these tests we have demonstrated that der energy functional at each SCF iteration, using the self- the two first methods are closely equivalent in terms of sec- (1) (1) consistent potential vscf,q and its Hermitian conjugate vscf,−q , ond derivatives quality and timing performance (calculations respectively. of phonons, electronic dielectric constant and Born effective The implementation of the algorithm has been tested by charges). However, the third one requires twice as many self comparing the spin susceptibility obtained from DFPT cal- consistent iterations to give the same quantities. On RuCl3 , culations with the finite-difference derivatives estimated from we have applied our methods to predict the change in phonon ground state DFT calculations with finite magnetic field. frequencies induced by the non-collinearity of the magnetic Specifically, for q = 0, we have performed calculations of moments. Comparing the frequencies estimated in collinear the longitudinal magnetization response of ferromagnetic bcc and non-collinear frameworks, we have shown the signature Fe as well as longitudinal and transverse spin susceptibili- of non-negligible NCM: certain phonon modes of specific ir- ties of antiferromagnetic Cr2 O3 . These two examples were reducible representations have their frequency modified by up chosen to test both metallic and insulating occupations. To to 10 cm−1 depending on the chosen direction (see Ref. 166 verify the validity of q 6= 0 implementation, we have further for more details). performed ground-state DFT calculations with harmonically The different methods of magnetization rotation can be set varying Zeeman field B = b cos(qr) with q = (1/2, 0, 0) and through the ixcrot A BINIT input flag. The default value is q = (1/4, 0, 0) for bcc Fe in corresponding supercells. The method 1 (ixcrot=1). The method 2 or 3 can be used by setting periodic parts of the first order density and magnetization ob- ixcrot to 2 or 3.167 Although in the corresponding A BINIT tained using the first-order finite-difference scheme were then version of Ref. 166 we implemented the NCM DFPT formal- compared to the corresponding results obtained from single ism at the Γ point only, further developments have been com- unit cell DFPT run. All considered test cases for q = 0 and pleted in order to account for the full q dependence in the q 6= 0 showed good agreement between DFPT and finite field phonon spectra (available in the next production version, only results, even in the case of coarse Brillouin zone sampling, with ixcrot=3). The implementation is valid for NCPSP only e.g. 2x2x2 and 4x4x4 k-point grids. Furthermore, to verify the but work is ongoing to extend it to the PAW formalism. accuracy of the approach we have performed the calculation of the Heisenberg exchange parameters J for bcc Fe. These are derived from the transverse spin susceptibility χ⊥ (q) and H. Spin-magnetic-field perturbation the corresponding Fourier components of the super-exchange parameters J(q) = χ(q)(−1) , yielding couplings of 18.9 meV The linear response of the magnetization to static and dy- (9.8 meV) for the first (second) shell of nearest neighbors, in namical magnetic fields (spin magnetic susceptibility χm ), is good agreement with previous literature170 . one of the key characteristics of magnetic materials. This All the test simulations described above were performed quantity is intrinsically linked to the spin fluctuation spec- within the LSDA approximation, which is currently the only trum and hence not only allows one to identify and charac- exchange-correlation functional implemented for linear spin terize the relevant magnetic excitations, but also to under- response calculations in A BINIT . Furthermore, for the q 6= 0, stand their character and estimate the strength of spin inter- case, it is imperative to use the explicit evaluation of the first- actions. The latter case appears to be particularly important order XC functional (ixcrot=3 A BINIT input flag). Further- in the view of the ever-growing need for accurate lattice spin more, at present and similarly to atomic displacement pertur- models to predict finite-temperature properties of spin sys- bations with non-collinear magnetism, the implementation of tems, and understand the physics of magnetic and topologi- zeeman magnetic field perturbation is limited to NCPSP for- cal materials. While the theoretical framework for linear spin malism. 17 VI. COMMUNITY PROJECTS: PSEUDOPOTENTIALS, HIGH-THROUGHPUT, AND THE LIBPAW LIBRARY A. Pseudopotentials and PAW data sets Planewave calculations are almost invariably performed in conjunction with some sort of pseudization scheme that freezes the inner core electrons and replaces their sharply varying wavefunctions in the region around the nucleus with smoother orbitals, which are easier to describe in Fourier space with a finite basis set7 . On the one hand, this approach allows one to reduce significantly the computational cost and enjoy the advantages of the planewave basis set: orthogonal- ity, systematic convergence, ease of implementation and ef- ficient Fast-Fourier-Transforms. On the other hand, as these pseudopotential operators are supposed to mimic the KS po- tential felt by the valence electrons in an all-electron DFT cal- culation, they imply a trade-off between accuracy and compu- FIG. 6. Comparison of the imaginary part of the dielectric function of silicon, directly related with the optical absorption, obtained in tational efficiency. Generating reliable and accurate tables of the present work (red and cyan curves) with Ref. 175 (“Marini”,blue pseudopotentials therefore represents a highly nontrivial task, curves), and experimental measurements (“Exp”, green curves)176 , especially for end-users who are not familiar with the pseu- for different temperatures. dopotential formalism and all its intricacies. For this reason, in the last few years, a significant effort has been made by the A BINIT community in order to provide users with a recom- mended set of pseudopotential tables that have been carefully I. Optical response with the Bethe-Salpeter equation, crafted and validated against ground-state reference results including temperature dependence obtained with all-electron codes177,178 . In what follows, we briefly describe the main features of the two official PAW and Within MBPT, one can address the computation of the ener- norm-conserving (NC) tables provided by the A BINIT group, gies and related characteristics of charged excitations as well and their design principles. as those of neutral excitations11 . For the former, the GW ap- a. JTH table To perform calculations in the frame of proximation is the state-of-the-art, while for the latter, one the PAW44 methodology, atomic data are needed, which cor- relies on the Bethe-Salpeter equation (BSE) which includes respond to the pseudopotentials in the norm-conserving ap- excitonic effects. Such approaches produce optical absorp- proach. We deliver a Mendeleev table of PAW data - the JTH tion spectra which are much more accurate than independent- table179 - on the A BINIT website. This table has been gener- particle approaches, such as the simple sum-over-states. ated thanks to the atomic code ATOMPAW 180 for 86 elements, The implementation of the BSE in A BINIT has been previ- from H to Rn. It is available within the scalar relativistic ap- ously described in Refs. 5 and 171. An efficient interpolation proximation, both for LDA and PBE XC functionals. The files allows on to cut down significantly the CPU time needed to are provided in a standard XML format for use in any PAW obtain optical spectra converged with respect to the Brillouin electronic structure code, following the specifications given in Zone sampling. Ref. 181. The XML format allows one to introduce new tags for special developments: for instance, pre-calculated atomic As a salient example, the A BINIT implementation has been matrix elements are given, to be used for calculations with used to obtain frequency-dependent Raman intensities172–174 , hybrid functionals. In v2.0 version of the JTH table, data by finite difference of optical spectra, even for the second- enabling calculations within the LDA-1/2182 approximation order Raman spectrum. are provided for some elements. Each file contains suggested The MBPT calculations can be combined with temperature planewave cut-offs, according to a requested accuracy. The dependent electronic structure presented above, Sec. V E, to input file of the ATOMPAW code is also provided for each ele- compute temperature-dependent optical spectra. The effect ment to allow the user to modify the parameters used for the of the electron-phonon interaction is included in the diagonal generation of the data, if needed. part of the BSE Hamiltonian according to: Hvck,v0 c0 k0 (T ) = b. PseudoDojo The P SEUDO D OJO project183 provides FA FA is the Hvck,v 0 c0 k0 +[∆εck (T ) − ∆εvk (T )] δvv0 δcc0 δkk0 , where H two predefined set of NC pseudopotentials: standard and “frozen-atom” expression for the BSE Hamiltonian intro- stringent. The standard set is designed for conventional duced in Ref. 5. Fig. 6 illustrates the agreement observed DFT and DFPT applications while the stringent version con- for silicon, when compared with experimental data, as well tains pseudopotentials with more electrons in valence, smaller as with Ref. 175. Note that the Bethe-Salpether Hamiltonian pseudization radii and improved scattering properties at high is not Hermitian, but can be tackled by an iterative Bi-Lanczos energies. The stringent set is recommended for performing algorithm5,174 . DFT/DFPT calculations in which high precision is needed, or 18 for many-body applications such as GW that are quite sen- code; the A BINIT group provides this framework in the sitive to the inclusion of semi-core states and to the quality A BI P Y project190 . These features have enabled the develop- of the logarithmic derivative in the empty region. The ma- ment of the following high-throughput projects. jority of the pseudopotentials include non-linear core correc- a. DFPT phonons The vibrational properties of a ma- tions184 with model core charges generated following Teter’s terial represent a fundamental ingredient for understanding a approach185 , which produces smooth core charges in real variety of physical phenomena and A BINIT provides accu- space with reasonably fast decay in reciprocal space. The in- rate and efficient algorithms based on the DFPT formalism clusion of the non-linear core correction improves the trans- (Sec. V). In order to fully exploit the features of the code and ferability of the pseudopotential and turned out to be crucial to handle a large number of calculations, a new python project avoid unphysical oscillations in the local part of the potential has been developed, A BI F LOWS . The package contains all in the region around the atom. These oscillations worsen the the functionalities required to execute high-throughput work- convergence rate, as well as the fulfillment of the acoustic sum flows with the F IRE W ORKS framework191 and to store the re- rule in phonon calculations, especially when GGA functionals sults in a MongoDB database in a standardized format for are employed183 . later queries and analysis. A BI F LOWS heavily relies on the At present, the P SEUDO D OJO provides NC tables gener- API implemented in A BI P Y , not only for the aforemen- ated with three different XC functionals: LDA, PBE and PBE- tioned properties, but also for the automatic generation of sol. For each XC flavor, one can opt for the scalar-relativistic A BINIT inputs for different kinds of calculations, and the or for the fully relativistic version, which includes the addi- post-processing of results. In order to provide sensible input tional projectors required to treat spin-orbit coupling186 . Fi- configurations for phonon calculations and avoid expensive nally, a specialized set of pseudopotentials for the lanthanides convergence studies for each material, a global study has been is also available, with f -electrons frozen in the core. These performed and a set of heuristic optimal parameters has been pseudopotentials are supposed to be used for calculations in determined192 . These parameters have been used to compute which the lanthanide is in the 3+ oxidation state. For each in an automated way the dynamical matrix, the phonon disper- pseudopotential, three different hints (low, normal and high) sion and the density of states for more than 2000 materials193 . for the planewave cutoff energy are suggested on the basis Taking advantage of the post-processing tools implemented of the convergence studies performed during the validation in A NADDB , a set of derived quantities has been obtained tests183 . These cutoff hints are employed in high-throughput for these materials. These include the dielectric permittivity studies, to implement machine-friendly workflows, but can tensors, Born effective charges, temperature dependent ther- also be used as a starting point in more conventional conver- modynamic properties (entropy, heat capacity, Helmholtz free gence studies. The pseudopotentials can be downloaded from energy, vibrational internal energy), frequency-dependent di- the official website or alternatively from a GitHub repository. electric tensor, speed of sound along high symmetry direc- The data is available in three different formats: the original tions, and the Debye-Waller tensor153 . All these quantities are psp8 file format implemented by Hamann, the UPF format freely available on the M ATERIALS P ROJECT website194 . and the psml format187 . The graphical interface of the of- b. DFPT SHG The workflow for the calculation of vi- ficial website allows users to select an element of the peri- brational properties discussed in the previous paragraph has odic table, the XC functional as well as the accuracy level and been extended to compute the optical dielectric tensor ε ∞ and the relativistic version. Jupyter notebooks with pre-generated the second-order susceptibility χ (2) (see Sec. V). This gives figures showing the convergence of physical properties as a access to optical properties beyond the linear regime with function of the cutoff energy and the results produced by the minimum effort. For instance, the second-order susceptibility ONCVPSP pseudopotential generator are also provided. The leads to the nonlinear Second Harmonic Generation (SHG) same interface can be used to download the PAW atomic data coefficient deff 195 . SHG processes play an important role in of the JTH table. modern optics, especially in laser-related science and technol- ogy196 . Only non-centrosymmetric materials are SHG-active and the calculation of the deff is possible in A BINIT through B. High-throughput: phonons, second-harmonic generation, the 2n + 1 theorem9,66 . and GW By selecting the non-centrosymmetric materials from databases197,198 in which the optical dielectric tensors were al- Using a first-principles package for high-throughput ap- ready available, we built up a new set of candidates including plications places more stringent demands on the implemen- more than 800 compounds aiming at calculating their second- tation details than the traditional approach based on single- order susceptibility. At this stage, we have computed the non- step operations. A BINIT implements two important fea- linear properties for more than 400 materials employing the tures to accommodate these needs: portable machine-readable above-mentioned A BI F LOWS workflow. The results of this output files in binary NETCDF format188,189 to communi- study will be presented in more detail in a future publication, cate results to external software, as well as automatic algo- and many more materials are being added to the list. rithms to select optimal parallelization settings at runtime. c. GW Computing the excited-state properties of ma- Besides these features directly implemented in the code it- terials (such as band gaps) implies going beyond standard self, efficient high-throughput calculations require a high- DFT and its fundamental shortcomings. MBPT provides a level library to programmatically interact with the Fortran rigorous and well-established formalism for accurate band- 19 structure calculations11 . The main ingredients of MBPT are BandStructureWork (w0) the Green’s function G, the screened Coulomb interaction W w0_t0 and the electronic self-energy Σ that is usually approximated ScfTask with the GW method199 . Automating GW calculations at the DEN high-throughput level represents, however, a significant chal- lenge. Compared to the DFT electron density in which only w0_t1 NscfTask occupied states are needed, the Green’s function requires the knowledge of all unoccupied states of the systems. In practi- WFK WFK WFK WFK WFK WFK cal implementations, the sum over an infinite number of bands Work (w1) is avoided by introducing a cutoff (nband) beyond which the contribution of the high-energy states is considered negligi- w1_t0 w1_t1 ble and the Fourier expansion of the two-point function W ScrTask ScrTask is restricted to G-vectors lying inside a sphere of kinetic en- ergy ecuteps. The required computational resources increase SCR SCR SCR SCR steeply when these two cutoffs are increased. An additional problem is that these two convergence parameters are linked, w2_t0 w3_t0 SigmaTask SigmaTask i.e. the convergence with respect to one cutoff is determined by the value of the other parameter with a rate that is system- w2_t1 w2_t2 w3_t1 w3_t2 dependent. The level of coupling, moreover, varies strongly SigmaTask SigmaTask SigmaTask SigmaTask from one material to another. Determining reliable input pa- SigmaConvWork (w2) SigmaConvWork (w3) rameters is, hence, of crucial importance to achieve precise and efficient MBPT calculations. FIG. 7. Schematic representation of a typical GW workflow. First To enable high-throughput GW calculations, an auto- (w0), the SCF density is calculated (DEN) followed by a non-SCF mated convergence framework has been developed using the calculation of a large number of unoccupied states (WFK). Then A BI P Y package200,201 . A typical workflow is illustrated in (w1), the latter are used to determine the screening (SCR). Finally Fig. 7. At a low density k-grid, the two interdependent cut- (w2), the GW corrections are computed for different k-points. In the offs are treated in the same loop ensuring convergence of both. graph, two screening calculations are represented corresponding to The convergence in k-space is treated separately since it was different sets of values for the two cutoff parameters. Each of them found to be unconnected. In the last step, the convergence is followed by three GW calculations corresponding to different k- rates of the screening cutoff and the unoccupied states cutoff points. at high k-space sampling are compared to those at low sam- pling density to finally ensure converged results. This method was used to compute the GW correction for about 90 crys- cently developed Gaussian-Type Atomic Orbital based DFT talline compounds200 . code204 . The use of LIBPAW in three codes, using three differ- ent types of basis functions (planewaves, wavelets, gaussians) highlights its portability. Other projects based on the LIBPAW C. The libpaw library are in progress. The library is currently built "on the fly" from A BINIT The A BINIT implementation of PAW44 — although not the source files. A standalone and A BINIT independent package only one — is one of the most complete available in a DFT is generated and can be directly inserted in a "host" code. At code. Ground-state properties (DFT), excited states (MBPT) present this is the only way to get the LIBPAW library; in the and response functions (DFPT) can be computed within PAW. near future the LIBPAW package will be directly downloadable However, the PAW method should theoretically be imple- as a separate package. mented in the same way in all DFT codes, regardless of the As soon as a host code uses the LIBPAW library, it can basis used to represent the electronic wavefunctions. access at low cost a full PAW implementation, and take From our original PAW implementation42 , we have devel- advantage of regular library updates, new features and de- oped a portable PAW library, called LIBPAW202 in order to bugging. Some features are automatically available be- facilitate the porting of PAW in other codes. This library pro- cause they are only implemented in "on-site" PAW contri- vides a formalism for computational physics — the pseudopo- butions such as, for example, Hubbard Hamiltonian calcula- tential/PAW approach — at a relatively high implementation tions (DFT+U). The library provides a module to access PAW level (not only low-level methods, such as I/O, error-handling atomic datasets in XML format and in A BINIT proprietary or systems resolution). It is packaged in such a way that it format (legacy). The XML format181 then gives access to sev- gives access to implementations of core PAW procedures with eral PAW atomic datasets tables: JTH179 , GBRV177 , ATOM - basis-independent data interfaces and thus can potentially be PAW 205 , and GPAW 206 . combined with any other basis used in different programs. The library fits very easily into the host code. A simple At present, the LIBPAW library — initially used in libpaw.h header file has to be adapted to integrate the LIB - A BINIT — has been integrated in at least two other codes: PAW . In the libpaw.h file, one has to specify which routines the wavelet-based BIGDFT software202,203 as well as a re- are used for the I/O and the error handling, as well as the de- 20 pendencies shared with the LIBPAW (eg, the LIBXC electronic • It allows one to substantially expand the user base and exchange and correlation library207 or NETCDF interfaces). thus to reinforce user feedback. Most of the low-level dependencies can be shared with the host code. The library is made of a complete collection of • It avoids developers re-coding what already exists, and Fortran 2003 modules. Therefore it can most easily be in- concentrates human work in implementations of novel corporated in codes written in Fortran. Each Fortran mod- physical quantities. ule is related to a specific "object" of the PAW formalism and contains associated methods. In the case of A BINIT , this code sharing effort is only in Among the provided objects, the following are worth men- its early stages; other features could be easily packaged and tioning (a more complete description is available in Ref. 202): shared: the handling of crystal symmetries, sampling of the Brillouin zone, self-consistent cycle preconditioning/mixing • m_pawrad: contains all functions related to the PAW algorithms, ... radial meshes and associated derivation/integration rou- tines, for different kinds of meshes (linear or logarith- mic). VII. CONCLUSION • m_pawpsp: used to read PAW atomic datasets (XML or A BINIT legacy formats), With the increasing use of DFT for materials science appli- cations, the need to have user-friendly packages which calcu- • m_pawtab: used to define non-self-consistent tabulated late and predict structural, electronic, vibrational, and elastic PAW data either read from the atomic dataset or directly properties is essential. In this paper, we have presented the deduced from it, case of A BINIT , an electronic structure package for materials and nanosystem simulations. Though A BINIT originally fo- • m_pawxc: computes the exchange and correlation po- cused on the “simple” solution of the ground state Kohn-Sham tential/energy in the PAW augmentation regions using equations of DFT, it has broadened to include many different developments over spherical harmonics, theories with far reaching applications, including DFT, DFPT, • m_pawcprj: calculates, stores and manipulates projec- TD-DFT, MBPT, and DMFT. Similarly, two different treat- tions of pseudo-wavefunctions Ψ̃nk on PAW non-local ments of core electrons are provided: norm-conserving pseu- projectors p̃i : cink = h p̃i |Ψ̃nk i. Note that the routine ded- dopotentials and projector-augmented waves. Each of these icated to the computation of these projections has to be methodologies implement a different set of properties, which provided by the host code, are summarized in the previous sections. A list of the ca- pabilities offered by A BINIT has been presented, with par- • m_pawrhoij: computes and manipulates the PAW oc- ticular attention given to the graphical interfaces created by j the A BI P Y workflow manager and the post-processing tool cupancy matrix ρi j = ∑nk ci∗ nk cnk , especially ensuring the use of symmetries, AGATE . In this paper, we have also concentrated on the properties • m_pawdij: computes and manipulates the PAW non- which make A BINIT rather unique. We started with DMFT, local potential intensities Di j (see Ref. 42 for a com- the application of finite electric fields, probes of nuclei prop- plete expression). In this routine all the physical in- erties (such as the Mössbauer and nuclear spectroscopy), and gredients are included (eg electronic correlations, spin- positron annihilation. We then moved to response functions, orbit, hybrid functionals, etc.). which are the most developed features of A BINIT : sections on Raman spectroscopy, electro-optic effects and strain perturba- All available methods are fully compatible with the Mes- tion, van der Waals interactions, flexoelectricity, the electron- sage Passing Interface (MPI) distributed parallelism, includ- phonon interaction, anharmonic interatomic force constants, ing parallelization over atomic sites, as well as over the real non-collinear magnetism, spin magnetic field and temperature space grid sampling the augmentation regions. A MPI com- dependence of the optical response with the Bethe-Salpeter municator has to be provided. equation. The LIBPAW library is a perfect example of "code shar- After examining individual properties, we presented some ing", in the perspective of outsourcing and sharing compo- of our efforts in community projects and library generation, nents common to DFT pseudopotential-based codes. The ad- which can be used by A BINIT and many other electronic vantages of code sharing are numerous: structure codes. For example, we mention the distribution of • It enforces clean and structured object-oriented pro- reliable and well tested PAW data sets, LIBPAW, the P SEU - gramming (not obvious in the case of Fortran). Strong DO D OJO project, and high-throughput developments to calcu- links to a specific code must be removed, and strong op- late advanced properties in the MBPT and DFPT formalisms. timization can be implemented once and used widely. In conclusion, we have presented here the evolution of A BINIT , from its history to the most recent developments. We • It allows one to implement a cross-validation process, have also included links to the code documentation and tutori- thanks to the use of the code in different contexts (here als. The basic theory and frameworks were discussed, which different electronic wavefunction bases). allow users to calculate a broad set of properties. Beyond the 21 features with unique implementations in A BINIT , we encour- Horizon 2020 research and innovation program (Grant Agree- age interested users to visit the A BINIT web page for more ment No. 724529). general documentation, tutorials, and information. Ph.G. acknowledges support from FNRS Belgium through PDR grant HiT4FiT, ULiège and the Communauté française de Belgique through the ARC project AIMED, the EU and VIII. DATA AVAILABILITY FNRS through M.ERA.NET project SIOX, and the European Funds for Regional Developments (FEDER) and the Wal- A large part of the data presented in this paper is available loon Region in the framework of the operational program directly from the Abinit Web page www.abinit.org. Any other “Wallonie-2020.EU” through the project Multifunctional thin data not appearing in this web page can be provided by the films/LoCoTED). corresponding author upon reasonable request. The Flatiron Institute is a division of the Simons Founda- tion. ACKNOWLEDGMENTS 1 P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864–B871 (1964). This work (A.H.R.) has been supported by the DMREF- 2 X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rig- NSF 1434897, National Science Fondation OAC-1740111 nanese, L. Sindic, M. Verstraete, G. Zérah, F. Jollet, M. Torrent, A. Roy, and U.S. Department of Energy DE-SC0016176 and DE- M. Mikami, P. Ghosez, J.-Y. Raty, and D. C. Allan, “First-principles com- SC0019491 projects. putation of material properties: The ABINIT software project,” Computa- tional Materials Science 25, 478–492 (2002). N.A.P. and M.J.V. gratefully acknowledge funding from the 3 X. Gonze, G.-M. Rignanese, M. Verstraete, J.-M. Beuken, Y. Pouillon, Belgian Fonds National de la Recherche Scientifique (FNRS) R. Caracas, F. Jollet, M. Torrent, G. Zerah, M. Mikami, P. Ghosez, M. Vei- under grant PDR T.1077.15-1/7. M.J.V. also acknowledges a then, J.-Y. Raty, V. Olevano, F. Bruneval, L. Reining, R. Godby, G. Onida, sabbatical “OUT” grant at ICN2 Barcelona, as well as ULiège D. Hamann, and D. Allan, “A brief introduction to the ABINIT software package,” Zeitschrift für Kristallographie - Crystalline Materials 220, 558– and the Communauté Française de Belgique (ARC AIMED 562 (2005). G.A. 15/19-09). 4 X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, X.G. and M.J.V. acknowledge funding from the FNRS un- P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, T. Deutsch, der grant number T.0103.19 - ALPS. L. Genovese, P. Ghosez, M. Giantomassi, S. Goedecker, D. R. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, X.G. and G.-M. R. acknowledge support from the Com- M. J. T. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, munauté française de Belgique through the SURFASCOPE D. Sangalli, R. Shaltaf, M. Torrent, M. J. Verstraete, G. Zérah, and J. W. project (ARC 19/24-057). Zwanziger, “ABINIT: First-principles approach to material and nanosys- X.G. acknowledges the hospitality of the CEA DAM-DIF tem properties,” Comput. Phys. Commun. 180, 2582–2615 (2009). 5 X. Gonze, F. Jollet, F. Abreu Araujo, D. Adams, B. Amadon, T. Applen- during the year 2017. court, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, E. Bousquet, G.H. acknowledges support from the U.S. Department of F. Bruneval, D. Caliste, M. Côté, F. Dahm, F. Da Pieve, M. Delaveau, Energy, Office of Science, United States, Office of Basic M. Di Gennaro, B. Dorado, C. Espejo, G. Geneste, L. Genovese, Energy Sciences, United States, Materials Sciences and En- A. Gerossier, M. Giantomassi, Y. Gillet, D. R. Hamann, L. He, G. Jo- gineering Division under Contract DE-AC02-05-CH11231: mard, J. Laflamme Janssen, S. Le Roux, A. Levitt, A. Lherbier, F. Liu, I. Lukačević, A. Martin, C. Martins, M. J. T. Oliveira, S. Poncé, Y. Pouil- Materials Project program KC23MP. lon, T. Rangel, G.-M. Rignanese, A. H. Romero, B. Rousseau, O. Rubel, The Belgian authors acknowledge computational resources A. A. Shukri, M. Stankovski, M. Torrent, M. Van Setten, B. Van Tro- from supercomputing facilities of the University of Liège, the eye, M. J. Verstraete, D. Waroquiers, J. Wiktor, B. Xu, A. Zhou, and Consortium des Equipements de Calcul Intensif (FRS-FNRS J. W. Zwanziger, “Recent developments in the ABINIT software package,” G.A. 2.5020.11) and Zenobe/CENAERO funded by the Wal- Comput. Phys. Commun. 205, 106–131 (2016). 6 M. P. Teter and M. C. Payne, “Apparatus and methods for predicting phys- loon Region under G.A. 1117545. ical and chemical properties of materials,” United States Patent , 6,106,562 M.C. and O.G. acknowledge support from the Fonds de (2000). 7 R. M. Martin, Electronic Structure (Cambridge University Press, 2004). Recherche du Québec Nature et Technologie (FRQ-NT), 8 S. Baroni, P. Giannozzi, and A. Testa, “Green’s-function approach to lin- Canada, and the Natural Sciences and Engineering Research ear response in solids,” Phys. Rev. Lett. 58, 1861–1864 (1987). Council of Canada (NSERC) under grant RGPIN-2016- 9 X. Gonze and J.-P. Vigneron, “Density-functional approach to nonlinear- 06666. response coefficients of solids,” Phys. Rev. B 39, 13120–13128 (1989). 10 P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, “Ab initio cal- The implementation of the LIBPAW library (M.T., T.R., D.C.) was supported by the ANR NEWCASTLE project, culation of phonon dispersions in semiconductors,” Phys. Rev. B 43, 7231 (1991). grant ANR- 2010-COSI-005-01 of the French National Re- 11 G. Onida, L. Reining, and A. Rubio, “Electronic excitations: Density- search Agency. functional versus many-body green’s-function approaches,” Rev. Mod. M.R. and M.S. acknowledge funding from Ministerio Phys. 74, 601–659 (2002). 12 R. M. Martin, L. Reining, and D. M. Ceperley, Interacting Electrons. de Economia, Industria y Competitividad (MINECO-Spain) through Grants No. MAT2016-77100-C2-2-P and No. SEV- Theory and Computational Approaches (Cambrige University Press, Cam- bridge, 2016). 2015-0496, and of Generalitat de Catalunya (Grant No. 2017 13 A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, “Dynamical SGR1506). This work has received funding from the Euro- mean-field theory of strongly correlated fermion systems and the limit of pean Research Council (ERC) under the European Union’s infinite dimensions,” Rev. Mod. Phys. 68, 13–125 (1996). 22 14 G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and algorithm,” Phys. Rev. B 55, 10337–10354 (1997). C. A. Marianetti, “Electronic structure calculations with dynamical mean- 34 X. Gonze and C. Lee, “Dynamical matrices, Born effective charges, di- field theory,” Rev. Mod. Phys. 78, 865–951 (2006). electric permittivity tensors, and interatomic force constants from density- 15 B. Amadon, F. Jollet, and M. Torrent, “γ and β cerium: LDA + U calcu- functional perturbation theory,” Phys. Rev. B 55, 10355–10368 (1997). lations of ground-state parameters,” Phys. Rev. B 77, 155104 (2008). 35 R. Resta, “Macroscopic polarization in crystalline dielectrics: The geo- 16 B. Amadon and B. Dorado, “A unified and efficient theory for the struc- metric phase approach,” Rev. Mod. Phys. 66, 899–915 (1994). tural properties of actinides and phases of plutonium,” Journal of Physics: 36 R. W. Nunes and X. Gonze, “Berry-phase treatment of the homogeneous Condensed Matter 30, 405603 (2018). electric field perturbation in insulators,” Phys. Rev. B 63, 155107 (2001). 17 X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, 37 I. Souza, J. Íñiguez, and D. Vanderbilt, “First-principles approach to insu- J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, N. Brouwer, F. Bruneval, lators in finite electric fields,” Phys. Rev. Lett. 89, 117602 (2002). G. Brunin, T. Cavignac, J.-B. Charraud, W. Chen, M. Côté, S. Cottenier, 38 J. W. Zwanziger, J. Galbraith, Y. Kipouros, M. Torrent, M. Giantomassi, J. Denier, G. Geneste, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, and X. Gonze, “Finite homogeneous electric fields in the projector aug- D. R. Hamann, G. Hautier, X. He, N. Helbig, N. Holzwarth, Y. Jia, F. Jol- mented wave formalism: Applications to linear and nonlinear response,” let, W. Lafargue-Dit-Hauret, K. Lejaeghere, M. A. Marques, A. Martin, Computational Materials Science 58, 113–118 (2012). C. Martins, H. P. Miranda, F. Naccarato, K. Persson, G. Petretto, V. Planes, 39 J. W. Zwanziger, “Computation of Mössbauer isomer shifts from first prin- Y. Pouillon, S. Prokhorenko, F. Ricci, G.-M. Rignanese, A. H. Romero, ciples,” J. Phys.: Condens. Matter 21, 195501 (2009). M. M. Schmitt, M. Torrent, M. J. van Setten, B. V. Troeye, M. J. Verstraete, 40 J. W. Zwanziger and M. Torrent, “First-principles calculation of electric G. Zérah, and J. W. Zwanziger, “The abinit project: Impact, environment field gradients in metals, semiconductors, and insulators,” Appl Magn Re- and recent developments,” Computer Physics Communications , 107042 son 33, 447–456 (2008). (2019). 41 J. W. Zwanziger, “First-principles study of the nuclear quadrupole res- 18 A. Georges, “Strongly correlated electron materials: Dynamical mean- onance parameters and orbital ordering in LaTiO3 ,” Phys. Rev. B 79, field theory and electronic structure,” in AIP Conference Proceedings (AIP, 033112 (2009). 2004). 42 M. Torrent, F. Jollet, F. Bottin, G. Zérah, and X. Gonze, “Implementation 19 B. Amadon, “A self-consistent DFT + DMFT scheme in the projector aug- of the projector augmented-wave method in the ABINIT code: Application mented wave method: Applications to cerium, Ce2 O3 and Pu2 O3 with the to the study of iron under pressure,” Computational Materials Science 42, hubbard I solver and comparison to DFT +U,” J. Phys.: Condens. Matter 337–351 (2008). 24, 075604 (2012). 43 G. Kresse and D. Joubert, “From ultrasoft pseudopotentials to the projector 20 B. Amadon, “First-principles DFT+DMFT calculations of structural prop- augmented-wave method,” Phys. Rev. B 59, 1758–1775 (1999). erties of actinides: Role of hund’s exchange, spin-orbit coupling, and crys- 44 P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, tal structure,” Phys. Rev. B 94, 115148 (2016). 17953–17979 (1994). 21 B. Amadon, F. Lechermann, A. Georges, F. Jollet, T. O. Wehling, and 45 M. Profeta, F. Mauri, and C. J. Pickard, “Accurate first principles predic- A. I. Lichtenstein, “Plane-wave based electronic structure calculations for tion of 17 O NMR parameters in SiO2 : Assignment of the zeolite ferrierite correlated materials using dynamical mean-field theory and projected local spectrum,” J. Am. Chem. Soc. 125, 541–548 (2003). orbitals,” Phys. Rev. B 77, 205112 (2008). 46 A. Honma, “Dipolar lattice-sums with applications to the exciton bands of 22 P. Werner, A. Comanac, L. de Medici, M. Troyer, and A. J. Millis, anthracene crystal and the crystal field due to point charges,” Journal of “Continuous-time solver for quantum impurity models,” Phys. Rev. Lett. the Physical Society of Japan 42, 1129–1135 (1977). 97, 076405 (2006). 47 M. J. Puska and R. M. Nieminen, “Theory of positrons in solids and on 23 J. Bieder and B. Amadon, “Thermodynamics of the α-γ transition in solid surfaces,” Rev. Mod. Phys. 66, 841–897 (1994). cerium from first principles,” Phys. Rev. B 89, 195132 (2014). 48 E. Boroński and R. M. Nieminen, “Electron-positron density-functional 24 D. Bergeron and A.-M. S. Tremblay, “Algorithms for optimized maximum theory,” Phys. Rev. B 34, 3820–3831 (1986). entropy and diagnostic tools for analytic continuation,” Phys. Rev. E 94, 49 F. Tuomisto and I. Makkonen, “Defect identification in semiconductors 023303 (2016). with positron annihilation: Experiment and theory,” Rev. Mod. Phys. 85, 25 O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Mes- 1583–1631 (2013). sio, and P. Seth, “TRIQS: A toolbox for research on interacting quantum 50 M. J. Puska, A. P. Seitsonen, and R. M. Nieminen, “Electron-positron car- systems,” Computer Physics Communications 196, 398 – 415 (2015). parrinello methods: Self-consistent treatment of charge densities and ionic 26 P. Seth, I. Krivenko, M. Ferrero, and O. Parcollet, “TRIQS/CTHYB: relaxations,” Phys. Rev. B 52, 10947–10961 (1995). A continuous-time quantum Monte Carlo hybridisation expansion solver 51 B. Barbiellini, M. J. Puska, T. Torsti, and R. M. Nieminen, “Gradient cor- for quantum impurity problems,” Computer Physics Communications 200, rection for positron states in solids,” Phys. Rev. B 51, 7341–7344 (1995). 274 – 284 (2016). 52 B. Barbiellini, M. J. Puska, T. Korhonen, A. Harju, T. Torsti, and 27 B. Amadon, T. Applencourt, and F. Bruneval, “Screened coulomb inter- R. M. Nieminen, “Calculation of positron states and annihilation in solids: action calculations: cRPA implementation and applications to dynamical A density-gradient-correction scheme,” Phys. Rev. B 53, 16201–16213 screening and self-consistency in uranium dioxide and cerium,” Phys. Rev. (1996). B 89, 125110 (2014). 53 J. Wiktor, G. Jomard, and M. Torrent, “Two-component density func- 28 F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. tional theory within the projector augmented-wave approach: Accurate Lichtenstein, “Frequency-dependent local interactions and low-energy ef- and self-consistent computations of positron lifetimes and momentum dis- fective models from electronic structure calculations,” Phys. Rev. B 70, tributions,” Phys. Rev. B 92, 125113 (2015). 195104 (2004). 54 J. Wiktor, G. Jomard, M. Torrent, and M. Bertolus, “Electronic structure 29 F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schönberger, “Calcula- investigation of energetics and positron lifetimes of fully relaxed monova- tions of Hubbard U from first-principles,” Phys. Rev. B 74, 125106 (2006). cancies with various charge states in 3C-SiC and 6H-SiC,” Phys. Rev. B 30 J.-B. Morée and B. Amadon, “First-principles calculation of coulomb in- 87, 235207 (2013). teraction parameters for lanthanides: Role of self-consistence and screen- 55 J. Wiktor, X. Kerbiriou, G. Jomard, S. Esnouf, M.-F. Barthe, and ing processes,” Phys. Rev. B 98, 205101 (2018). M. Bertolus, “Positron annihilation spectroscopy investigation of vacancy 31 B. Amadon and A. Gerossier, “Comparative analysis of models for the clusters in silicon carbide: Combining experiments and electronic struc- α−γ phase transition in cerium: A DFT+DMFT study using wannier or- ture calculations,” Phys. Rev. B 89, 155203 (2014). bitals,” Phys. Rev. B 91, 161103 (2015). 56 J. Wiktor, G. Jomard, and M. Bertolus, “Electronic structure calculations 32 O. Gingras, R. Nourafkan, A.-M. S. Tremblay, and M. Côté, “Supercon- of positron lifetimes in SiC: Self-consistent schemes and relaxation effect,” ducting Symmetries of Sr2 RuO4 from First-Principles Electronic Struc- Nucl. Instrum. Methods Phys. Res., Sect. B 327, 63–67 (2014). ture,” Phys. Rev. Lett. 123, 217005 (2019). 57 J. Wiktor, M.-F. Barthe, G. Jomard, M. Torrent, M. Freyss, and M. Berto- 33 X. Gonze, “First-principles responses of solids to atomic displacements lus, “Coupled experimental and DFT+U investigation of positron lifetimes and homogeneous electric fields: Implementation of a conjugate-gradient in UO2 ,” Phys. Rev. B 90, 184101 (2014). 23 58 J. Wiktor, G. Jomard, M. Torrent, and M. Bertolus, “First-principles calcu- 83 D. Sando, P. Hermet, J. Allibe, J. Bourderionnet, S. Fusil, C. Carrétéro, lations of momentum distributions of annihilating electron–positron pairs E. Jacquet, J.-C. Mage, D. Dolfi, A. Barthélémy, P. Ghosez, and M. Bibes, in defects in UO2 ,” J. Phys.: Condens. Matter 29, 035503 (2016). “Linear electro-optic effect in multiferroic BiFeO3 thin films,” Phys. Rev. 59 S. Hagiwara, Y. Suzuki, and K. Watanabe, “Spin-polarized annihilation B 89, 195106 (2014). lifetime of positron of d0 ferromagnetism in gallium nitride: A two- 84 M. Veithen and P. Ghosez, “Temperature dependence of the electro-optic component density functional theory simulation,” Appl. Phys. Express 9, tensor and refractive indices of BaTiO3 from first principles,” Phys. Rev. 041001 (2016). B 71, 132101 (2005). 60 A. Marinopoulos, “Positron lifetimes of bare and hydrogenated zirconium 85 X. Wu, D. Vanderbilt, and D. R. Hamann, “Systematic treatment of dis- vacancies in cubic yttria-stabilized zirconia: an ab initio study,” J. Phys.: placements, strains, and electric fields in density-functional perturbation Condens. Matter 31, 315503 (2019). theory,” Phys. Rev. B 72, 035105 (2005). 61 S. Hagiwara and K. Watanabe, “Positron states at li-and o-adsorbed fe 86 R. D. King-Smith and D. Vanderbilt, “Theory of polarization of crystalline (001) ferromagnetic surfaces studied by two-component density functional solids,” Phys. Rev. B 47, 1651–1654 (1993). theory,” J. Phys. Soc. Japan 85, 114703 (2016). 87 D. R. Hamann, K. M. Rabe, and D. Vanderbilt, “Generalized-gradient- 62 W. Humphrey, A. Dalke, and K. Schulten, “VMD – Visual Molecular functional treatment of strain in density-functional perturbation theory,” Dynamics,” J. Mol. Graph. 14, 33–38 (1996). Phys. Rev. B 72, 033102 (2005). 63 D. R. Hamann, X. Wu, K. M. Rabe, and D. Vanderbilt, “Metric tensor 88 A. Martin, M. Torrent, and R. Caracas, “Projector augmented-wave for- formulation of strain in density-functional perturbation theory,” Phys. Rev. mulation of response to strain and electric-field perturbation within density B 71, 035117 (2005). functional perturbation theory,” Phys. Rev. B 99, 094112 (2019). 64 X. Gonze, “Adiabatic density-functional perturbation theory,” Phys. Rev. 89 S. Grimme, “Semiempirical GGA-type density functional constructed with A 52, 1096–1114 (1995). a long-range dispersion correction,” J. Comput. Chem. 27, 1787–1799 65 X. Gonze, P. Ghosez, and R. W. Godby, “Density-polarization functional (2006). theory of the response of a periodic insulating solid to an electric field,” 90 S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and ac- Phys. Rev. Lett. 74, 4035–4038 (1995). curate ab initio parametrization of density functional dispersion correction 66 M. Veithen, X. Gonze, and P. Ghosez, “Nonlinear optical susceptibilities, (DFT-D) for the 94 elements H-Pu,” The Journal of Chemical Physics 132, Raman efficiencies, and electro-optic tensors from first-principles density 154104 (2010). functional perturbation theory,” Phys. Rev. B 71, 125107 (2005). 91 H. Rydberg, M. Dion, N. Jacobson, E. Schröder, P. Hyldgaard, S. I. Simak, 67 G. Placzek, “Marx handbuch der radiologie,” Akademische Ver- D. C. Langreth, and B. I. Lundqvist, “Van der waals density functional for lagsgesellschaft, Leipzig (1934). layered structures,” Phys. Rev. Lett. 91, 126402 (2003). 68 M. Cardona and G. Güntherodt, Light Scattering in Solids II (Springer, 92 M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, 1982). “Van der waals density functional for general geometries,” Phys. Rev. Lett. 69 P. Brüesch, Phonons: Theory and Experiments II (Springer, 1986). 92, 246401 (2004). 70 W. D. Johnston, “Nonlinear optical coefficients and the Raman scattering 93 A. Tkatchenko and M. Scheffler, “Accurate molecular van der Waals inter- efficiency of LO and TO phonons in acentric insulating crystals,” Phys. actions from ground-state electron density and free-atom reference data,” Rev. B 1, 3494–3503 (1970). Phys. Rev. Lett. 102, 073005 (2009). 71 M. Veithen, X. Gonze, and P. Ghosez, “First-principles study of the 94 A. M. Reilly and A. Tkatchenko, “van der waals dispersion interactions electro-optic effect in ferroelectric oxides,” Phys. Rev. Lett. 93, 187401 in molecular materials: beyond pairwise additivity,” Chem. Sci. 6, 3289– (2004). 3301 (2015). 72 K. Miwa, “Prediction of Raman spectra with ultrasoft pseudopotentials,” 95 B. Van Troeye, M. Torrent, and X. Gonze, “Interatomic force constants Phys. Rev. B 84, 094304 (2011). including the DFT-D dispersion contribution,” Phys. Rev. B 93, 144304 73 B. Lafuente, R. Downs, H. Yang, and N. Stone, “The power of databases: (2016). The RRUFF project,” in Highlights in Mineralogical Crystallography 96 B. Van Troeye, M. J. van Setten, M. Giantomassi, M. Torrent, G.-M. Rig- (Walter de Gruyter GmbH, 2016) pp. 1–29, http://rruff.info/ (sample id. nanese, and X. Gonze, “First-principles study of paraelectric and ferro- : R040031). electric CsH2 PO4 including dispersion forces: Stability and related vibra- 74 S. A. Prosandeev, U. Waghmare, I. Levin, and J. Maslar, “First-order tional, dielectric, and elastic properties,” Phys. Rev. B 95, 024112 (2017). 0 00 97 C. P. Romao, “Anomalous thermal expansion and chiral phonons in Raman spectra of AB1/2 B1/2 O3 double perovskites,” Phys. Rev. B 71, 214307 (2005). BiB3 O6 ,” Phys. Rev. B 100, 060302 (2019). 75 R. Caracas and R. E. Cohen, “Theoretical determination of the Raman 98 M. Born and K. Huang, Dynamical Theory of Crystal Lattices (Oxford spectra of MgSiO3 perovskite and post-perovskite at high pressure,” Geo. University Press, Oxford, 1954). 99 R. Resta, “Towards a bulk theory of flexoelectricity,” Phys. Rev. Lett. 105, Phys. Res. Lett. 33, L12S05 (2006). 76 P. Hermet, M. Veithen, and P. Ghosez, “First-principles calculations of the 127601 (2010). 100 J. Hong and D. Vanderbilt, “First-principles theory of frozen-ion flexoelec- nonlinear optical susceptibilities and Raman scattering spectra of lithium niobate,” Journal of Physics: Condensed Matter 19, 456202 (2007). tricity,” Phys. Rev. B 84, 180101(R) (2011). 77 P. Hermet, M. Veithen, and P. Ghosez, “Raman scattering intensities in 101 J. Hong and D. Vanderbilt, “First-principles theory and calculation of flex- BaTiO3 and PbTiO3 prototypical ferroelectrics from density functional oelectricity,” Phys. Rev. B 88, 174107 (2013). 102 M. Stengel, “Flexoelectricity from density-functional perturbation theory,” theory,” J. Phys.: Condens. Matter 21, 215901 (2009). 78 H. Djani, P. Hermet, and P. Ghosez, “First-principles characterization of Physical Review B 88, 174106 (2013). 103 M. Stengel, “Microscopic response to inhomogeneous deformations in the P21 ab ferroelectric phase of Aurivillius Bi2 WO6 ,” J. Phys. Chem. C 118, 13514 (2014). curvilinear coordinates,” Nature Communications 4, 2693 (2013). 79 P. Hermet and P. Ghosez, “First-principles study of the dynamical and non- 104 M. Stengel, “Surface control of flexoelectricity,” Phys. Rev. B 90, linear optical properties of urea single crystals,” Phys. Chem. Chem. Phys. 201112(R) (2014). 105 M. Stengel and D. Vanderbilt, “First-principles theory of flexoelectricity,” 12, 835–843 (2010). 80 N. A. Pike, A. Dewandre, B. Van Troeye, X. Gonze, and M. J. Verstraete, in Flexoelectricity in Solids From Theory to Applications, edited by A. K. “Vibrational and dielectric properties of the bulk transition metal dichalco- Tagantsev and P. V. Yudin (World Scientific Publishing Co., Singapore, genides,” Phys. Rev. Materials 2, 063608 (2018). 2016) Chap. 2, pp. 31–110. 81 N. A. Pike, A. Dewandre, B. Van Troeye, X. Gonze, and M. J. Verstraete, 106 M. Stengel and D. Vanderbilt, “Quantum theory of mechanical deforma- “Vibrational and dielectric properties of the monolayer transition metal tions,” Phys. Rev. B 98, 125133 (2018). 107 A. Schiaffino, C. E. Dreyer, D. Vanderbilt, and M. Stengel, “Metric-wave dichalcogenides,” Phys. Rev. Materials 3, 074009 (2019). 82 R. Caracas, “Database of computed Raman and infrared spectra and other approach to flexoelectricity within density-functional perturbation theory,” physical properties of minerals,” http://www.wurm.info, 2019-11-26. Phys. Rev. B 99, 085107 (2019). 24 108 C. E. Dreyer, M. Stengel, and D. Vanderbilt, “Current-density implemen- 132 J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and tation for calculating flexoelectric coefficients,” Phys. Rev. B 98, 075153 S. G. Louie, “BerkeleyGW: A massively parallel computer package for (2018). the calculation of the quasiparticle and optical properties of materials and 109 M. Royo and M. Stengel, “First-principles theory of spatial dispersion: nanostructures,” Computer Physics Communications 183, 1269 – 1289 Dynamical quadrupoles and flexoelectricity,” Phys. Rev. X 9, 021050 (2012). (2019). 133 Z. Li, G. Antonius, M. Wu, F. H. da Jornada, and S. G. Louie, 110 R. M. Martin, “Piezoelectricity,” Phys. Rev. B 5, 1607–1613 (1972). “Electron-phonon coupling from ab initio linear-response theory within 111 F. Giustino, “Electron-phonon interactions from first principles,” Rev. the gw method: Correlation-enhanced interactions and superconductivity Mod. Phys. 89, 015003 (2017). in Ba1−x Kx BiO3 ,” Phys. Rev. Lett. 122, 186402 (2019). 112 P. B. Allen, “New method for solving Boltzmann’s equation for electrons 134 G. M. Rignanese, J. P. Michenaud, and X. Gonze, “Ab initio study of in metals,” Phys. Rev. B 17, 3725–3734 (1978). the volume dependence of dynamical and thermodynamical properties of 113 S. Y. Savrasov and D. Y. Savrasov, “Electron-phonon interactions and re- silicon,” Phys. Rev. B 53, 4488–4497 (1996). lated physical properties of metals from linear-response theory,” Phys. Rev. 135 G. Leibfried and W. Ludwig, “Solid state physics vol 12 : Theory of an- B 54, 16487–16501 (1996). harmonic effects in crystals,” (Academic Press, New York, 1961) pp. 275– 114 B. Xu and M. J. Verstraete, “First principles explanation of the positive 444. Seebeck coefficient of lithium,” Phys. Rev. Lett. 112, 196603 (2014). 136 A. A. Maradudin and A. E. Fein, “Scattering of neutrons by anharmonic 115 P. B. Allen and V. Heine, “Theory of the temperature dependence of elec- crystal,” Phys. Rev. 128, 2589 (1962). tronic band structures,” J. Phys. C: Solid State Phys. 9, 2305–2312 (1976). 137 R. Cowley, “Anharmonics crystals,” Rep. Prog. Phys. 31, 123 (1968). 116 P. B. Allen and M. Cardona, “Temperature dependence of the direct gap of 138 T. Barron and M. Klein, “Dynamical properties of solids : Volume 1. crys- Si and Ge,” Phys. Rev. B 27, 4760–4769 (1983). talline solids, fundamentals,” (North Holland, Amsterdam, 1974) p. 391. 117 S. Poncé, G. Antonius, P. Boulanger, E. Cannuccia, A. Marini, M. Côté, 139 P. Brüesch, Phonons : Theory and Experiments I (Springer-Verlag, Berlin and X. Gonze, “Verification of first-principles codes: Comparison of total Heidelberg, 1982). energies, phonon frequencies, electron–phonon coupling and zero-point 140 D. C. Wallace, Thermodynamics of Crystals (Dover, New York, 1997). motion correction to the gap between ABINIT and QE/Yambo,” Compu- 141 D. C. Wallace, Statistical physics of crystals and liquids (World Scientific tational Materials Science 83, 341–348 (2014). Publishing Company, Singapore, 2003). 118 G. Antonius, S. Poncé, P. Boulanger, M. Côté, and X. Gonze, “Many-body 142 A. Togo, F. Oba, and I. Tanaka, “First-principles calculations of the ferroe- effects on the zero-point renormalization of the band structure,” Phys. Rev. lastic transition between rutile-type and CaCl2 -type SiO2 at high pressure,” Lett. 112, 215501 (2014). Phys. Rev. B 78, 134106 (2008). 119 G. Antonius, S. Poncé, E. Lantagne-Hurtubise, G. Auclair, X. Gonze, and 143 A. Carreras, A. Togo, and I. Tanaka, “Dynaphopy : A code for extract- M. Côté, “Dynamical and anharmonic effects on the electron-phonon cou- ing phonon quasiparticles from molecular dynamics simulations,” Comp. pling and the zero-point renormalization of the electronic structure,” Phys. Phys. Comm. 221, 221 (2017). Rev. B 92, 085137 (2015). 144 T. Tadano, Y. Gohda, and S. Tsuneyuki, “Anharmonic force constants 120 S. Poncé, Y. Gillet, J. Laflamme Janssen, A. Marini, M. Verstraete, and extracted from first-principles molecular dynamics: applications to heat X. Gonze, “Temperature dependence of the electronic structure of semi- transfer simulations,” J. Phys. Condens. Matter 26, 225012 (2014). conductors and insulators,” The Journal of Chemical Physics 143, 102813 145 I. Errea, M. Calandra, and F. Mauri, “First-principles theory of anhar- (2015). monicity and the inverse isotope effect in superconducting palladium- 121 S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, “Phonons hydride compounds,” Phys. Rev. Lett. 111, 177002 (2013). and related crystal properties from density-functional perturbation theory,” 146 L. J. Nelson, G. L. W. Hart, F. Zhou, and V. Ozolins, “Compressive sens- Rev. Mod. Phys. 73, 515–562 (2001). ing as a paradigm for building physics models,” Phys. Rev. B 87, 035125 122 J. P. Nery, P. B. Allen, G. Antonius, L. Reining, A. Miglio, and X. Gonze, (2013). “Quasiparticles and phonon satellites in spectral functions of semiconduc- 147 A. Glensk, B. Grabowski, T. Hickel, and J. Neugebauer, “Breakdown of tors and insulators: Cumulants applied to the full first-principles theory the Arrhenius law in describing vacancy formation energies: The impor- and the Fröhlich polaron,” Phys. Rev. B 97, 115145 (2018). tance of local anharmonicity revealed by ab initio thermodynamics,” Phys. 123 F. Brown-Altvater, G. Antonius, T. Rangel, M. Giantomassi, C. Draxl, Rev. X 4, 011018 (2014). X. Gonze, S. G. Louie, and J. B. Neaton, “Band gap renormalization, 148 C. Toher, J. J. Plata, O. Levy, M. de Jong, M. Asta, M. B. Nardelli, and carrier mobilities, and the electron-phonon self-energy in crystalline naph- S. Curtarolo, “High-throughpout computational screening of thermal con- thalene,” Accepted in PRB XX, YYYYYY (2020). ductivity, Debye temperature, and Grüneisen parameters using a quasihar- 124 A. Eiguren and C. Ambrosch-Draxl, “Wannier interpolation scheme for monic Debye model,” Phys. Rev. B 90, 174107 (2014). phonon-induced potentials: Application to bulk MgB2 , W, and the (1 × 1) 149 K. Esfarjani and H. T. Stokes, “Method to extract anharmonic force con- H-covered W(110) surface,” Phys. Rev. B 78, 045124 (2008). stants from first principles calculations,” Phys. Rev. B 77, 144112 (2008). 125 P. Vogl, “Microscopic theory of electron-phonon interaction in insulators 150 O. Hellman, I. A. Abrikosov, and S. I. Simak, “Lattice dynamics of anhar- or semiconductors,” Physical Review B 13, 694–704 (1976). monic solids from first principles,” Phys. Rev. B 84, 180301(R) (2011). 126 C. Verdi and F. Giustino, “Fröhlich electron-phonon vertex from first prin- 151 O. Hellman, P. Steneteg, I. A. Abrikosov, and S. I. Simak, “Temperature ciples,” Phys. Rev. Lett. 115, 176401 (2015). dependent effective potential method for accurate free energy calculations 127 J. Sjakste, N. Vast, M. Calandra, and F. Mauri, “Wannier interpolation of solids,” Phys. Rev. B 87, 104111 (2013). of the electron-phonon matrix elements in polar semiconductors: Polar- 152 O. Hellman and I. A. Abrikosov, “Temperature-dependent effective third- optical coupling in GaAs,” Physical Review B 92, 054307 (2015). order interatomic force constants from first principles,” Phys. Rev. B 88, 128 X. Gonze, P. Boulanger, and M. Côté, “Theoretical approaches to the tem- 144301 (2013). perature and zero-point motion effects on the electronic band structure,” 153 C. Lee and X. Gonze, “Ab initio calculation of the thermodynamic prop- Ann. Phys. 523, 168–178 (2010). erties and atomic temperature factors SiO2 α-quartz and stishovite,” Phys. 129 Y.-H. Chan, G. Antonius, and S. G. Louie, “Polaron spectral properties Rev. B 51, 8610–8613 (1995). in ZnO and SrTiO3 from first principles,” To be submitted XX, YYYYYY 154 J. Shiomi, K. Esfarjani, and G. Chen, “Thermal conductivity of half- (2020). heusler compounds from first-principles calculations,” Phys. Rev. B 84, 130 S. Ponce, E. R. Margine, and F. Giustino, “Towards predictive many- 104302 (2011). body calculations of phonon-limited carrier mobilities in semiconductors,” 155 C. H. Lee and C. K. Gan, “Anharmonic interatomic force constants and Phys. Rev. B 97, 121201(R) (2018). thermal conductivity from Grüneisen parameters : An application to 131 G. K. H. Madsen, J. Carrete, and M. J. Verstraete, “BoltzTraP2, a program graphene,” Phys. Rev. B 96, 035105 (2017). for interpolating band structures and calculating semi-classical transport 156 J. Bouchet and F. Bottin, “Thermal evolution of vibrational properties of coefficients,” Computer Physics Communications 231, 140–145 (2018). α-U,” Phys. Rev. B 92, 174108 (2015). 25 157 A. Dewaele, V. Stutzmann, J. Bouchet, F. Bottin, F. Occelli, and J. R. Yates, G.-X. Zhang, and S. Cottenier, “Reproducibility in density M. Mezouar, “High pressure-temperature phase diagram and equation of functional theory calculations of solids,” Science 351, aad3000 (2016), state of titanium,” Phys. Rev. B 91, 134108 (2015). https://science.sciencemag.org/content/351/6280/aad3000.full.pdf. 158 J. Bouchet and F. Bottin, “High-temperature and high-pressure phase tran- 179 F. Jollet, M. Torrent, and N. Holzwarth, “Generation of projector sitions in uranium,” Phys. Rev. B 95, 054113 (2017). augmented-wave atomic data: A 71 element validated table in the XML 159 B. Dorado, F. Bottin, and J. Bouchet, “Phonon spectra of plutonium at format,” Comput. Phys. Commun. 185, 1246–1254 (2014). high temperatures,” Phys. Rev. B 95, 104303 (2017). 180 N. A. W. Holzwarth, A. R. Tackett, and G. E. Matthews, “A projector 160 J. Bouchet, F. Bottin, V. Recoules, F. Remus, G. Morard, R. M. Bolis, and augmented wave (PAW) code for electronic structure calculations, Part I: A. Benuzzi-Mounaix, “Ab initio calculations of the B1-B2 phase transition Atompaw for generating atom-centered functions,” Comput. Phys. Com- in MgO,” Phys. Rev. B 99, 094113 (2019). mun. 135, 329–347 (2001). 161 F. Bottin, J. Bouchet, and J. Bieder, “ A -TDEP : Temperature dependent 181 PAW-XML, “XML specification for atomic PAW datasets,” https://esl. effective potential for abinit. thermodynamic properties using second and cecam.org/Paw-xml, 2019-11-26. third order interatomic force constants,” Submitted to Comput. Phys. Com- 182 L. G. Ferreira, M. Marques, and L. K. Teles, “Approximation to den- mun XX, YYYYYY (2019). sity functional theory for the calculation of band gaps of semiconductors,” 162 M. J. Verstraete, M. Torrent, F. Jollet, G. Zérah, and X. Gonze, “Den- Phys. Rev. B 78, 125116 (2008). sity functional perturbation theory with spin-orbit coupling: Phonon band 183 M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. structure of lead,” Phys. Rev. B 78, 045119 (2008). Hamann, X. Gonze, and G.-M. Rignanese, “The PseudoDojo: Training 163 J. E. Peralta, G. E. Scuseria, and M. J. Frisch, “Noncollinear magnetism and grading a 85 element optimized norm-conserving pseudopotential ta- in density functional calculations,” Phys. Rev. B 75, 125119 (2007). ble,” Comput. Phys. Commun. 226, 39–54 (2018). 164 D. Hobbs, G. Kresse, and J. Hafner, “Fully unconstrained noncollinear 184 S. G. Louie, S. Froyen, and M. L. Cohen, “Nonlinear ionic pseudopoten- magnetism within the projector augmented-wave method,” Phys. Rev. B tials in spin-density-functional calculations,” Phys. Rev. B 26, 1738–1742 62, 11556–11570 (2000). (1982). 165 J. Kubler, K. H. Hock, J. Sticht, and A. R. Williams, “Density functional 185 M. Teter, “Additional condition for transferability in pseudopotentials,” theory of non-collinear magnetism,” Journal of Physics F: Metal Physics Phys. Rev. B 48, 5031–5041 (1993). 18, 469 (1988). 186 D. R. Hamann, “Optimized norm-conserving Vanderbilt pseudopoten- 166 F. Ricci, S. Prokhorenko, M. Torrent, M. J. Verstraete, and E. Bousquet, tials,” Phys. Rev. B 88, 085117 (2013). “Density functional perturbation theory within noncollinear magnetism,” 187 A. García, M. J. Verstraete, Y. Pouillon, and J. Junquera, “The psml format Phys. Rev. B 99, 184404 (2019). and library for norm-conserving pseudopotential data curation and inter- 167 Note that the meaning of ixcrot was unfortunately inverted in Ref. 166 operability,” Computer Physics Communications 227, 51–71 (2018). with respect to the actual implementation. 188 X. Gonze, C.-O. Almbladh, A. Cucca, D. Caliste, C. Freysoldt, M. Mar- 168 S. Y. Savrasov, “Linear response calculations of spin fluctuations,” Phys. ques, V. Olevano, Y. Pouillon, and M. Verstraete, “Specification of an Rev. Lett. 81, 2570–2573 (1998). extensible and portable file format for electronic structure and crystallo- 169 K. Cao, H. Lambert, P. G. Radaelli, and F. Giustino, “Ab initio calcu- graphic data,” Computational Materials Science 43, 1056–1065 (2008). lation of spin fluctuation spectra using time-dependent density functional 189 D. Caliste, Y. Pouillon, M. J. Verstraete, V. Olevano, and X. Gonze, “Shar- perturbation theory, plane waves, and pseudopotentials,” Phys. Rev. B 97, ing electronic structure and crystallographic data with ETSF_IO,” Com- 024420 (2018). puter Physics Communications 179, 748–758 (2008). 170 M. Pajda, J. Kudrnovský, I. Turek, V. Drchal, and P. Bruno, “Ab initio 190 M. Giantomassi and et al., “Abipy project,” https://github.com/abinit/ calculations of exchange interactions, spin-wave stiffness constants, and abipy. Curie temperatures of Fe, Co, and Ni,” Phys. Rev. B 64, 174402 (2001). 191 A. Jain, S. P. Ong, W. Chen, B. Medasani, X. Qu, M. Kocher, M. Braf- 171 Y. Gillet, M. Giantomassi, and X. Gonze, “Efficient on-the-fly interpola- man, G. Petretto, G.-M. Rignanese, G. Hautier, D. Gunter, and K. A. tion technique for Bethe-Salpeter calculations of optical spectra,” Comput. Persson, “FireWorks: A dynamic workflow system designed for high- Phys. Commun. 203, 83–93 (2016). throughput applications,” Concurrency Computat.: Pract. Exper. 27, 5037– 172 Y. Gillet, M. Giantomassi, and X. Gonze, “First-principles study of exci- 5059 (2015). tonic effects in Raman intensities,” Phys. Rev. B 88, 094305 (2013). 192 G. Petretto, X. Gonze, G. Hautier, and G.-M. Rignanese, “Convergence 173 Y. Gillet, S. Kontur, M. Giantomassi, C. Draxl, and X. Gonze, “Ab initio and pitfalls of density functional perturbation theory phonons calculations approach to second-order resonant Raman scattering including exciton- from a high-throughput perspective,” Computational Materials Science phonon interaction,” Sci Rep 7, 7344 (2017). 144, 331–337 (2018). 174 Y. Gillet, Ab initio study of Raman and optical spectra of crystalline 193 G. Petretto, S. Dwaraknath, H. P.C. Miranda, D. Winston, M. Giantomassi, materials and their temperature dependence, Ph.D. thesis, Université M. J. van Setten, X. Gonze, K. A. Persson, G. Hautier, and G.-M. Rig- catholique de Louvain, Louvain-la-neuve, Belgium (2017). nanese, “High-throughput density-functional perturbation theory phonons 175 A. Marini, “Ab initio finite-temperature excitons,” Phys. Rev. Lett. 101, for inorganic materials,” Sci. Data 5, 180065 (2018). 106405 (2008). 194 A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cho- 176 G. E. Jellison and F. A. Modine, “Optical functions of silicon between 1.7 lia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, “The Materials and 4.7 ev at elevated temperatures,” Phys. Rev. B 27, 7466 (1983). Project: A materials genome approach to accelerating materials innova- 177 K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, “Pseudopo- tion,” APL Materials 1, 011002 (2013). tentials for high-throughput DFT calculations,” Computational Materials 195 S. K. Kurtz and T. T. Perry, “A powder technique for the evaluation of Science 81, 446–452 (2014). nonlinear optical materials,” Journal of Applied Physics 39, 3798–3813 178 K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, (1968). D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, S. de Giron- 196 J. M. Rondinelli and E. Kioupakis, “Predicting and designing optical prop- coli, T. Deutsch, J. K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, erties of inorganic materials,” Annual Review of Materials Research 45, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, 491–518 (2015). M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, 197 I. Petousis, D. Mrdjenovich, E. Ballouz, M. Liu, D. Winston, W. Chen, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, T. Graf, T. D. Schladt, K. A. Persson, and F. B. Prinz, “High-throughput D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, screening of inorganic compounds for the discovery of novel dielectric and E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Mars- optical materials,” Scientific Data 4, 160134 (2017). man, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, 198 F. Naccarato, F. Ricci, J. Suntivich, G. Hautier, L. Wirtz, and G.-M. C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, Rignanese, “Searching for materials with high refractive index and wide G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, band gap: A first-principles high-throughput study,” Phys. Rev. Materials S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, 3, 044602 (2019). D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M. Wills, 26 199 L. Hedin, “New method for calculating the one-particle Green’s function der, “Daubechies wavelets as a basis set for density functional pseudopo- with application to the electron-gas problem,” Phys. Rev. 139, A796–A823 tential calculations,” The Journal of Chemical Physics 129, 014109 (2008). (1965). 204 X.-G. Xiong and T. Yanai, “Projector augmented wave method incor- 200 M. J. van Setten, M. Giantomassi, X. Gonze, G.-M. Rignanese, and porated into Gauss-type atomic orbital based density functional theory,” G. Hautier, “Automation methodologies and large-scale validation for Journal of Chemical Theory and Computation 13, 3236–3249 (2017), GW: Towards high-throughput GW calculations,” Phys. Rev. B 96, 155207 https://doi.org/10.1021/acs.jctc.7b00404. (2017). 205 ATOMPAW, “Atompaw atomic datasets,” http://users.wfu.edu/natalie/ 201 M. Giantomassi, M. Stankovski, R. Shaltaf, M. Grüning, F. Bruneval, papers/pwpaw/PAWDatasets.html, 2019-11-26. P. Rinke, and G.-M. Rignanese, “Electronic properties of interfaces and 206 GPAW, “GPAW atomic PAW setups,” https://wiki.fysik.dtu.dk/gpaw/ defects from many-body perturbation theory: Recent developments and setups/setups.html, 2019-06-07. applications,” phys. stat. sol. (b) 248, 275–289 (2011). 207 S. Lehtola, C. Steigemann, M. J. T. Oliveira, and M. A. L. Marques, “Re- 202 T. Rangel, D. Caliste, L. Genovese, and M. Torrent, “A wavelet-based pro- cent developments in (l)ibxc — (a) comprehensive library of functionals jector augmented-wave (PAW) method: Reaching frozen-core all-electron for density functional theory,” SoftwareX 7, 1–5 (2018). precision with a systematic, adaptive and localized wavelet basis set,” Comput. Phys. Commun. 208, 1–8 (2016). 203 L. Genovese, A. Neelov, S. Goedecker, T. Deutsch, S. A. Ghasemi, A. Wil- land, D. Caliste, O. Zilberberg, M. Rayson, A. Bergman, and R. Schnei-