abinit is probably the first electronic-structure package to have been released under an open-source license about 20 years ago. It implements density functional theory, density-functional perturbation theory (DFPT), many-body perturbation theory (*GW* approximation and Bethe–Salpeter equation), and more specific or advanced formalisms, such as dynamical mean-field theory (DMFT) and the “temperature-dependent effective potential” approach for anharmonic effects. Relying on planewaves for the representation of wavefunctions, density, and other space-dependent quantities, with pseudopotentials or projector-augmented waves (PAWs), 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 abinit 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, and orbital magnetization); positron annihilation; Raman intensities and electro-optic effect; and DFPT calculations of response to strain perturbation (elastic constants and piezoelectricity), spatial dispersion (flexoelectricity), electronic mobility, temperature dependence of the gap, and spin-magnetic-field perturbation. The abinit 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 datasets, high-throughput calculations (databases of phonon band structure, second-harmonic generation, and *GW* computations of bandgaps), and the library libpaw. abinit has strong links with many other software projects that are briefly mentioned.

## I. INTRODUCTION

Since the introduction of density functional theory (DFT) in 1964,^{1} the field of electronic structure calculations has changed profoundly. This theory became the most popular electronic structure method used to characterize materials at the atomic scale and has given rise to different research efforts that have been pushed further by the use and applications of DFT. One of the reasons why this theory has become the workhorse of material characterization is the wide distribution of computational packages, where solutions of the Kohn–Sham equations are implemented. Every package has a different scheme and philosophy, but all of them have in common the generation of a software that is user-friendly and solves the DFT equations. In this respect, abinit^{2–5} was one of the first free licensed solid-state electronic structure packages on the market.

In this paper, we give first a global overview of abinit at the level of its history, community, and impact, as well as its main capabilities. Then, we focus on selected capabilities that are perceived as rather specific among the set of available first-principles packages.

After describing the historical development of the project and its impact (Sec. II), we present an overview of abinit capabilities through a collection of keywords and concepts, with entry points in the documentation (Sec. III). This should allow users to know whether abinit is capable of delivering some specific property of materials.

Sections IV and V single out some specifics of abinit. Section IV focuses on ground-state and electronic properties: correlated materials [Dynamical Mean-Field Theory (DMFT)], Sec. IV A; treatment of finite electric field, Sec. IV B; properties at the atomic nuclei, Sec. IV C; and positron annihilation, Sec. IV D. Section V focuses on response properties: Raman intensities, Sec. V A; responses to strain, Sec. V B; responses including van der Waals interactions, Sec. V C; spatial dispersion (flexoelectricity and dynamical quadrupoles), Sec. V D; electron–phonon coupling, Sec. V E; temperature-dependent properties of materials in the anharmonic regime, Sec. V F; responses of solids presenting non-collinear magnetism (NCM), Sec. V G; response to spin-magnetic perturbation, Sec. V H; and temperature-dependent optical spectra, Sec. VI. Finally, we spend a few words on our recent community efforts related to accurate and efficient Norm-Conserving Pseudopotential (NCPSP) as well as Projector-Augmented Wave (PAW) atomic data (Sec. VI A), to high-throughput calculations (phonon band structure, second-harmonic generation, and accurate GW bandgap calculations; Sec. VI B), and to the libpaw library (Sec. VI C).

## II. THE PROJECT: HISTORY, COMMUNITY, AND IMPACT

The abinit project can be traced back to an initial effort in the late 1980s by Allan, supervised by Michael P. Teter, at Corning Incorporated and Cornell University. The code used norm-conserving pseudopotentials and a planewave representation of the Kohn–Sham orbitals, with the local-density approximation (LDA) as the only exchange–correlation functional available, and was written in FORTRAN 77. In 1990, Gonze joined at Cornell, and they implemented density functional perturbation theory (DFPT) on top of what was then dubbed the corning code in a separate application, named respfn (Corning Incorporated is a large American company focused on glass and materials applications). After a few years in the hands of Biosym, Inc., where the name was changed to plane_wave and the code was sold commercially, the development of corning 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 recognition that worldwide collaboration was necessary to develop the code. The previous codes (plane_wave, respfn, and corning) became the pillars of this new implementation. Corning Incorporated agreed to release the source of plane_wave to support this effort and agreed not to enforce a patent they held on the preconditioned conjugate gradient algorithm.^{6} The code was rewritten in Fortran 90 with parallel features, originally named dft2000 but soon changed to abinit (September 1998).

The first version of abinit was made public in March 1999, primarily to beta testers. Initially, abinit was only capable of finding the total energy, the electronic charge density, and the electronic structures of periodic systems using pseudopotentials and a planewave basis. In July 1999, the full response function capability was also available. This implementation was a significant step, as it was one of the only codes that allow users to calculate phonons, dielectric constants, Born effective charges, etc. The first publicly available release of abinit was made in December 1999 under the GNU General Public License. In June 2000, an international advisory committee was selected to help with the strategy, support, and management of the code. The first abinit developer workshop took place in Louvain-la-Neuve in 2002. Since then, the developer abinit meeting takes place every two years. The 9th developer workshop took place in Louvain-la-Neuve in 2019, with the participation of around 60 speakers.

The abinit spirit is not only to offer the community a free license code but also to encourage users and developers to reuse parts of the computational package. The source code is always available, and if a user is interested in introducing specific implementations into the package, an account on GitHub suffices such that changes can be merged and tested in the official version. The abinit community has created a series of methods and tests to guarantee the stability of the code, precisely to avoid the introduction of errors by new implementations. abinit is a very well documented code with an extensive description of installation, details of the input variables and their different dependencies, a very established tutorial with specific applications and a full set of examples (more than 800). To support general users, abinit has created a forum, which is linked to the abinit website, and allows users to ask questions related to theory, implementations, or use of abinit. The forum is a primary resource for users, and both developers and users answer questions and offer advice.

To spread the use of abinit and to train young scientists in the field of electronic structure, the abinit developers have been supporting schools around the world. While the first abinit workshop was in 2002, there are usually between 2 and 3 schools and workshops per year in many different places. There is a regular presence at the March Meeting of the American Physical Society, where half day tutorials are organized by the abinit community, and at CECAM. Some of the CECAM events are theory oriented workshops, with discussions of the formalism and algorithms behind abinit, and examples of calculations are discussed. On the other hand, there are also regular CECAM hands-on schools, which allow users to gain experience in the use of abinit and also learn the theoretical background. These schools are very well attended but depend on specific locations with sufficient access to computational facilities. In all cases, the presentations of the abinit developer workshops, schools, and conferences are stored in the abinit web page and are accessible globally. These presentations offer the users details of the code but are also an archive of the different implementations and how they have been tested and used.

The abinit community champions free software development but also wants to keep the integrity and modularity of the computational package. All new developments need to be approved by the community, to ensure to keep track of the different efforts, and to conduct a minimal screening on the dependability and ethical behavior of people involved in the code development. The driving idea is to preserve a constructive and supporting environment for all scientists willing to collaborate. All major implementation efforts in abinit are connected to individual publications, but they are also introduced through specific publications about every three–five years, where the most recent research efforts are described. With this method, abinit tries to acknowledge the developers’ achievements and use these publications to describe how the new implementations are part of the broader package. The abinit attitude toward publications is based on research independence, but also collaborative efforts that continue to grow the computational package.

The research interest of the developers group (around 50 people at the time of writing) are respected by the abinit community and overseen by an international advisory committee. This committee consists of around ten senior scientists, which have been involved in the development of abinit and have been supportive of the philosophy and efforts of the community.

A metric of abinit’s impact can be obtained by using bibliometric analysis accessed from the Web of Science (data taken on December 26th, 2019). This can be assessed by adding all citations of all the main papers reporting abinit developments, i.e., Refs. 2–5. We found a total of 5125 citations since 2002 of which 3995 areunique, with an average number above 250 total citations per year since 2011 (see Fig. 1). Using this database of unique entries, we can address the impact of abinit on different fields of science by counting some of the papers reported in high impact journals. For example, there are 160 papers from *Physical Review Letters*, 109 from *Journal of Chemical Physics*, 108 from *Computational Materials Science*, 32 from *Nano Letters*, 33 from *Scientific Reports*, 18 from *Nature Communications*, 8 from *Nature Materials*, and 7 from *Nature*. As is clear from this analysis, abinit is a healthy package that impacts several scientific fields and supports the work of many agencies and institutions, with authors of 84 different countries and from more than 2000 different institutions.

## III. OVERVIEW OF FORMALISMS AND PROPERTIES

An overview of the capabilities of a large scientific code can be structured in different ways. The present work targets users who simply wish to know whether abinit is capable of computing a particular property of a material or nanosystem and possibly representing it graphically. In this respect, Tables I and II give an alphabetic list of keywords (or concepts) that best represent a capability of abinit, with related documentation: possibly some section of the present paper and/or reference to some “topics” of the on-line abinit Web documentation https://docs.abinit.org/topics/features and/or reference to some scripts and illustrations in the abipy gallery http://abinit.github.io/abipy/gallery/index.html. Highlighted (bold) keywords refer to specific sections of the present paper. The emphasis on these keywords stems from the belief that such characteristics/capabilities of abinit are not commonly found among first-principles codes and can be put forward as examples of what makes abinit unique.

Keyword/concept . | Documentation . |
---|---|

Atomic mean square displacement (phonons) | topic:Temperature, plot_phonons_msqd.html |

Bader atomic charges | topic:Bader |

Born effective charges | Section VI B, topic:DFPT |

Born effective charges, non-collinear magnetism case | Section V G |

Correlated electronic state (dynamic—DMFT and static—DFT+ U) | Section 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 | Section VI B, topic:DFPT, topic:Phonons |

Dielectric permittivity, non-collinear magnetism case | Section 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 | Section V D |

Effective mass | topic:EffectiveMass |

Elastic tensor | Sections V B and V C, topic:DFPT, topic:Elastic |

Elastic tensor, temperature-dependent | Section V F |

Electric field (finite, using Berry phase) | Section IV B, topic:Berry |

Electric field gradient (at nuclei) | Section IV C, topic:EFG |

Electro-optic effect | Section V A, topic:nonlinear |

Electric polarization (Berry phase) | Section 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) | Section IV A, topic:DMFT |

Electron–phonon coupling strength | topic:ElPhonInt, plot_a2f.html |

Electronic bandgap (GW) | topic:GW, plot_qps.html |

Electronic bandgap, temperature-dependent | Section V E, topic:TDepES |

Electronic band structure | Section 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 | Section V D |

Free energy (phonons) | Section 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 |

Keyword/concept . | Documentation . |
---|---|

Atomic mean square displacement (phonons) | topic:Temperature, plot_phonons_msqd.html |

Bader atomic charges | topic:Bader |

Born effective charges | Section VI B, topic:DFPT |

Born effective charges, non-collinear magnetism case | Section V G |

Correlated electronic state (dynamic—DMFT and static—DFT+ U) | Section 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 | Section VI B, topic:DFPT, topic:Phonons |

Dielectric permittivity, non-collinear magnetism case | Section 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 | Section V D |

Effective mass | topic:EffectiveMass |

Elastic tensor | Sections V B and V C, topic:DFPT, topic:Elastic |

Elastic tensor, temperature-dependent | Section V F |

Electric field (finite, using Berry phase) | Section IV B, topic:Berry |

Electric field gradient (at nuclei) | Section IV C, topic:EFG |

Electro-optic effect | Section V A, topic:nonlinear |

Electric polarization (Berry phase) | Section 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) | Section IV A, topic:DMFT |

Electron–phonon coupling strength | topic:ElPhonInt, plot_a2f.html |

Electronic bandgap (GW) | topic:GW, plot_qps.html |

Electronic bandgap, temperature-dependent | Section V E, topic:TDepES |

Electronic band structure | Section 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 | Section V D |

Free energy (phonons) | Section 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 |

Keyword/concept . | Documentation . |
---|---|

Infrared reflectivity | topic:Phonons |

Interatomic force constants | topic:Phonons |

Internal energy (phonons) | Section 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 | Section V H |

Molecular dynamics | topic:MolecularDynamics, topic:DynamicsMultibinit, topic:PIMD |

Mössbauer isomer shift | Section IV C, topic:EFG |

Optical absorption | Section VI, topic:BSE, topic:Optic |

Optical response, temperature-dependent | Section VI, topic:BSE, topic:Optic |

Phonon bands | Sections V C and VI B, topic:Phonons, topic:DFPT, |

plot_phonons_lo_to.html, plot_phonons.html, | |

plot_ddb_asr.html, topic:Band2eps | |

Phonon bands, temperature-dependent | Section V F, topic:Tdep |

Phonon bands, non-collinear magnetic case | Section 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 | Section V F, topic:Tdep |

Phonon linewidth (electron–phonon coupling) | topic:PhononWidth |

Piezoelectric tensor | Sections V B and V C, topic:DFPT, topic:Elastic |

Positron annihilation | Section 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 | Section V A, topic:nonlinear, topic:Phonons |

Refraction index | Section VI, topics:DFPT, topics:Optic, topic:BSE |

Resistivity (electron–phonon coupling) | Section V E, topic:ElPhonTransport |

Scanning tunneling microscopy map | topic:STM |

Stopping power of charged particles | topic:RandStopPow |

Second harmonic generation | Section VI B, topic:nonlinear, topic:Optic, plot_optic.html |

Sound velocity | Section VI B, topic:PhononBands, plot_speed_of_sound.html |

Specific heat (phonons) | Section VI B, topic:Temperature, plot_phthermo.html |

Structural relaxation | topic:GeoOpt, topic:ForcesStresses, plot_hist.html |

Superconducting transition temperature | Section V E, topic:ElPhonTransport |

Thermal expansion | topic:Temperature |

Thermodynamic properties | Section VI B, topic:Temperature |

Thermodynamic properties including anharmonicities | Section 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 bandgap | Section V E, topic:TDepES |

Keyword/concept . | Documentation . |
---|---|

Infrared reflectivity | topic:Phonons |

Interatomic force constants | topic:Phonons |

Internal energy (phonons) | Section 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 | Section V H |

Molecular dynamics | topic:MolecularDynamics, topic:DynamicsMultibinit, topic:PIMD |

Mössbauer isomer shift | Section IV C, topic:EFG |

Optical absorption | Section VI, topic:BSE, topic:Optic |

Optical response, temperature-dependent | Section VI, topic:BSE, topic:Optic |

Phonon bands | Sections V C and VI B, topic:Phonons, topic:DFPT, |

plot_phonons_lo_to.html, plot_phonons.html, | |

plot_ddb_asr.html, topic:Band2eps | |

Phonon bands, temperature-dependent | Section V F, topic:Tdep |

Phonon bands, non-collinear magnetic case | Section 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 | Section V F, topic:Tdep |

Phonon linewidth (electron–phonon coupling) | topic:PhononWidth |

Piezoelectric tensor | Sections V B and V C, topic:DFPT, topic:Elastic |

Positron annihilation | Section 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 | Section V A, topic:nonlinear, topic:Phonons |

Refraction index | Section VI, topics:DFPT, topics:Optic, topic:BSE |

Resistivity (electron–phonon coupling) | Section V E, topic:ElPhonTransport |

Scanning tunneling microscopy map | topic:STM |

Stopping power of charged particles | topic:RandStopPow |

Second harmonic generation | Section VI B, topic:nonlinear, topic:Optic, plot_optic.html |

Sound velocity | Section VI B, topic:PhononBands, plot_speed_of_sound.html |

Specific heat (phonons) | Section VI B, topic:Temperature, plot_phthermo.html |

Structural relaxation | topic:GeoOpt, topic:ForcesStresses, plot_hist.html |

Superconducting transition temperature | Section V E, topic:ElPhonTransport |

Thermal expansion | topic:Temperature |

Thermodynamic properties | Section VI B, topic:Temperature |

Thermodynamic properties including anharmonicities | Section 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 bandgap | Section V E, topic:TDepES |

In the “topics” tags of the online abinit Web documentation, we try to address the generic challenge of software documentation by providing for each topic a hub to the underlying theory (including bibliographical references), to the related abinit input variables, to the example abinit input files, and to the abinit tutorials, also mentioning possible restrictions in the implementation or different levels of accuracy of the implementation.

abinit implements different first-principles formalisms, i.e., DFT,^{7} DFPT,^{3,8–10} time-dependent density functional theory (TD-DFT),^{11} many-body perturbation theory (MBPT),^{12} and DMFT.^{13,14} Each of these formalisms allows one to address different properties at different levels of accuracy. As an example, the electronic bandgap can be computed from DFT, but the predicted value is notoriously inaccurate. Instead, the bandgap is better computed using the GW approximation within MBPT, also implemented in abinit. Moreover, even within MBPT-GW, different levels of accuracy are available: the simple one-shot G_{0}W_{0} approach with different plasmon pole models or more sophisticated self-consistent flavors of GW (SCGW). Similarly, different families of exchange–correlations functionals can be used in DFT: LDA, generalized-gradient approximation (GGA), meta-GGA, hybrid functionals, and van der Waals (vdW) corrected functionals. It will not be the purpose of this section to specify for each property the different levels of accuracy. By contrast, such information might be found, thanks to the on-line abinit topics.

Similarly, within some formalism, the abinit implementation might not be compatible with both the norm-conserving pseudopotential (NCPSP) approach and the PAW approach. In addition, the treatment of non-spin-polarized systems, collinear magnetic systems, collinear antiferromagnetic systems, and non-collinear magnetic systems might not all be available.

In Ref. 4, we provided such a more detailed description of the formalism+implementation level for many of the abinit capabilities (see Tables 1–4 of Ref. 4). abinit has evolved since then, and the best available source of information is now the on-line topics.

abinit also has several graphical post-processors, including abipy and agate. A gallery of abipy examples (with associated Jupyter notebooks) is available at http://abinit.github.io/abipy/gallery/index.html. The links to these examples are also mentioned in Tables I and II.

Finally, note that some of the capabilities of abinit are presently not yet available/documented in the latest public release of abinit (v8) (in particular those in Sec. IV D, part of IV E and part of VI B) but will be so in the forthcoming one (v9). Some of these are mentioned in a section of the present paper, but there is no associated topic or abipy example. Some of the common variables used within this text are found in Table III.

Variable . | Description . |
---|---|

B | Magnetic field |

c | Speed of light |

C_{0} | Clamped-ion elastic tensor |

C_{6,IJ} | Dispersion coefficient |

D^{/bfq} | Dynamic matrix |

e | Electronic charge |

e_{0} | Clamped-ion piezoelectric tensor |

E | Total energy |

$E$ | Electric field |

ε^{∞} | Optical dielectric tensor |

E_{γ} | Energy of a nuclear transition |

$Ece\u2212p$ | Electron–positron correlation functional |

f | Band occupancy |

f_{mk+q} | Fermi–Dirac occupation factor |

|i⟩ | Polarization vector along direction i |

i, j | Reduced coordinates |

n(ω), n_{qν} | Bose occupation factor |

n_{α} | Refractive index along α |

P | Polarization |

r^{e} | Classic electron radius |

r | Linear electro-optic tensor |

U | Onsite Coulomb interaction energy |

|u_{nk}⟩ | Bloch wavefunction |

u_{m} | Eigen-displacement vector |

U_{0} | Ground-state energy |

⟨v^{2}⟩ | Mean square velocity of the nucleus |

v_{ext} | External potential |

v_{nk} | Particle velocity |

Z | Atomic number |

Z^{*} | Born effective charge tensor |

α, β | Cartesian coordinates |

α_{m} | Mode dependent Raman tensor |

Γ | Mode independent line-width |

γ | Lattice-strain coupling tensor |

δ | Mössbauer isomer shift |

Δ⟨r^{2}⟩ | Change in size of the nucleus |

η_{αβ} | Strain component |

$\mu \xafII$ | Frozen ion flexelectricity |

τ | Positron lifetime |

χ^{(1)} | First order electric susceptibility |

χ^{(2)} | Second order susceptibility |

$\Phi (p)$ | Interatomic force constants |

|Ψ⟩ | Planewave basis functions |

|Φ⟩ | All electron wavefunctions |

$|\Phi \u0303$ | Pseudized atomic wavefunctions |

ω_{l} | Laser light frequency |

Ω_{0} | Unit cell volume |

Variable . | Description . |
---|---|

B | Magnetic field |

c | Speed of light |

C_{0} | Clamped-ion elastic tensor |

C_{6,IJ} | Dispersion coefficient |

D^{/bfq} | Dynamic matrix |

e | Electronic charge |

e_{0} | Clamped-ion piezoelectric tensor |

E | Total energy |

$E$ | Electric field |

ε^{∞} | Optical dielectric tensor |

E_{γ} | Energy of a nuclear transition |

$Ece\u2212p$ | Electron–positron correlation functional |

f | Band occupancy |

f_{mk+q} | Fermi–Dirac occupation factor |

|i⟩ | Polarization vector along direction i |

i, j | Reduced coordinates |

n(ω), n_{qν} | Bose occupation factor |

n_{α} | Refractive index along α |

P | Polarization |

r^{e} | Classic electron radius |

r | Linear electro-optic tensor |

U | Onsite Coulomb interaction energy |

|u_{nk}⟩ | Bloch wavefunction |

u_{m} | Eigen-displacement vector |

U_{0} | Ground-state energy |

⟨v^{2}⟩ | Mean square velocity of the nucleus |

v_{ext} | External potential |

v_{nk} | Particle velocity |

Z | Atomic number |

Z^{*} | Born effective charge tensor |

α, β | Cartesian coordinates |

α_{m} | Mode dependent Raman tensor |

Γ | Mode independent line-width |

γ | Lattice-strain coupling tensor |

δ | Mössbauer isomer shift |

Δ⟨r^{2}⟩ | Change in size of the nucleus |

η_{αβ} | Strain component |

$\mu \xafII$ | Frozen ion flexelectricity |

τ | Positron lifetime |

χ^{(1)} | First order electric susceptibility |

χ^{(2)} | Second order susceptibility |

$\Phi (p)$ | Interatomic force constants |

|Ψ⟩ | Planewave basis functions |

|Φ⟩ | All electron wavefunctions |

$|\Phi \u0303$ | Pseudized atomic wavefunctions |

ω_{l} | Laser light frequency |

Ω_{0} | Unit cell volume |

## IV. GROUND STATE AND ELECTRONIC PROPERTIES

### A. Correlated materials: Dynamical-mean field theory

Systems with localized orbitals, such as transition metals: lanthanides or actinides, exhibit strong correlation effects and are difficult to describe using DFT with currently available functionals.^{14} One way to improve their description is to explicitly include the onsite Coulomb interaction *U* between the electrons in the correlated orbitals of the system (e.g., the *d* or *f* shells).

The DFT+*U* method, available in abinit,^{15–17} treats this interaction statically and is most efficient on magnetic Mott insulators. To describe in a more coherent way, systems with various interaction strengths, the method of choice is DMFT.^{13} This method solves the local many-body problem for a given correlated atom in the effective field created by the other atoms. This field is self-consistently related to the solution of the atomic impurity problem. The combination of DFT and DMFT^{14,18} enables the description of realistic systems with both correlated and non-correlated electrons. The DFT+DMFT method has been helpful to improve, in particular, spectral functions (describing both Hubbard bands and quasiparticle peaks in the spectral function), total energy (for iron systems, actinides, or lanthanides), and magnetic properties (see e.g., Ref. 14).

The method has been implemented in abinit^{5,19,21} and is presented in a tutorial. The key points of the available DFT+DMFT implementation are the following:

Correlated orbitals are defined as projected local orbitals Wannier functions,

^{21}and their localization can easily be changed by modifying the associated energy windows.The method is fully self-consistent with respect to the electronic density (see Fig. 2).

^{19}Impurity solvers directly available in abinit are the Hubbard I method and the continuous time quantum Monte Carlo

^{22}(using either a diagonal or a general hybridization function) in the simplified but efficient density–density approximation.^{19,23}Spin–orbit coupling calculations are possible using a real valued imaginary time hybridization function.Internal energy and electronic entropy can be computed.

^{19,23}Impurity and

*k*-resolved spectral functions can be computed using analytical continuation of the Green’s function or the self-energy with an external code (e.g., OmegaMaxent^{24}). Figure 3 shows the LDA+DMFT spectral function of*α*-cerium, as computed in abinit.Several parallelization schemes can be used,

^{17}allowing application of the method to large systems.abinit is coupled to the CTHYB impurity solver of the TRIQS library

^{25,26}through a C++ interface. This allows one to solve the impurity problem in the fully rotationally invariant formulation of the interacting Hamiltonian.abinit is built with a python invocation scheme that allows any personalized python script to solve the impurity problem. Thus, an experienced user can invoke their favorite solver within a DFT+DMFT calculation from abinit.

DFT+DMFT uses effective interactions parameters *U* and *J* as an input. They can be computed in abinit^{27} on the same correlated Wannier orbitals using the constrained Random Phase Approximation (cRPA) method.^{28,29} In cRPA,^{29} as the screening arising from electronic transitions among the correlated shell is a correlation effect already described in the CTQMC solution of the impurity problem, only other electronic transitions are used to compute the screening.^{29} The cRPA method in abinit is documented by a tutorial.

These schemes were recently used in *f* electron systems to compute the *f* effective interactions parameters,^{20,27,30} to improve the description of their spectral or structural properties,^{16,20,23,31} and to improve our understanding of the superconducting symmetries in Sr_{2}RuO_{4}.^{32}

### B. Finite electric fields

abinit can be used to compute the response to electric fields in several ways. From DFPT, the susceptibility can be calculated^{33,34} as the second-order derivative of the total energy at zero field $\u22022E/\u2202Ei\u2202EJ$. Alternatively, the first order derivative is the electrical polarization **P**, which may be non-zero even in the absence of external fields, e.g., in ferroelectric materials. Computing this term appears formally equivalent to computing ∫*d***r r***n*(**r**), that is, the dipole moment density, but this expression is ill-defined in an extended system. This problem is solved by the modern theory of polarization^{35} in which the key insight is the recognition that while the polarization is not well-defined, its derivative with respect to a change in the system is, and so the polarization may be computed (up to a constant) through integration. The key formula is

Equation (1) is implemented in abinit using both NCPSP and the PAW formalism. In both cases, the wavefunction derivatives are computed using a finite-difference scheme so that a coherent phase relationship between the **k**-points may be ensured. From this approach, the polarization may be computed.

With a scheme in place to compute the polarization, the response to a finite electric field may be computed by adding a term $\u2212P\u22c5E$ to the total energy.^{36} Then, its effect is included in the self-consistent energy minimization cycle through the gradient $\delta (\u2212P\u22c5E)/\delta \u2009unk|$.^{37,38} In this way, Kohn–Sham states that minimize the total energy including the electric field is found [the polarization from Eq. (1) is updated at each step as the Kohn–Sham states evolve], provided that the field is not so strong that the insulating gap breaks down.

Using this approach, we have computed a number of responses, including linear terms that permit validation against DFPT, and nonlinear responses that are only available using finite field calculations.^{38}

### C. Properties at the atomic nuclei

abinit provides several features for computing properties that arise from the electronic structure very near the atomic nuclei. These features are of use for computing and interpreting a variety of experimental probes, in particular Mössbauer spectroscopy,^{39} and nuclear magnetic and nuclear quadrupole resonance spectroscopy.^{40,41}

Such observables involve the overlap of the electronic wavefunctions with the nucleus, in the case of the Mössbauer isomer shift, and the gradient of the electric field at the nucleus, observed in both Mössbauer spectra and magnetic resonance experiments. As both observables depend sensitively on details of the electronic structure very close to the nuclear position, these are cases where a planewave-only treatment with NCPSP is quite inaccurate and it is necessary to use the PAW formalism.

The Mössbauer isomer shift *δ*, in velocity units, is determined by the overlap of the electronic density with a nucleus undergoing a nuclear state transition and is given by

for the electron density of the absorber (*A*) and source (*S*) at the nucleus (here located at 0). In this formula, *E*_{γ}, the energy of a nuclear transition, and Δ⟨*r*^{2}⟩, the change in size of the nucleus, are nuclear properties that must be supplied by the user, while abinit can determine the electronic contributions *ρ*_{A} and *ρ*_{S}. To do this, the PAW formalism is used in which an observable ⟨*A*⟩ is computed as^{42–44}

For properly constructed PAW data, all-electron accuracy may be recovered in this formalism. For computation of *ρ*(0) appearing in Eq. (2), we use as the observable the 3-D delta-function *δ*(0).^{39}

Using Eq. (3), isomer shifts may also be computed, given available data for *E*_{γ} and Δ⟨*r*^{2}⟩ for the transition in question. To validate the method, a range of shifts were computed and then Δ⟨*r*^{2}⟩ extracted by comparison to known shifts; we found excellent agreement for a variety of nuclei, including tin, germanium, and zinc.^{39} The zinc case is a very stringent test because the isomer shift range is so small and the transition so sharp; accurate modeling requires inclusion of the secondary Doppler shift, which we computed by first principles from the phonon dispersion relations. The secondary Doppler shift takes the form *S*_{D} = −⟨*v*^{2}⟩/2*c* and is obtainable once the full phonon dispersion curves are computed. We found that even in this case, the PAW-derived Δ⟨*r*^{2}⟩ was in excellent agreement with experimental values.^{39}

The electric field gradient calculation proceeds along similar lines, except the observable of interest is now the second derivative of the potential *n*(**r**)/|**r** − **r**′| for density *n*(**r**). The implementation in abinit^{40} is similar to that in other codes, such as quantum espresso.^{45} As usual in PAW, the density is decomposed into a planewave-based part and PAW-sphere corrections. In addition, there is a charge density due to the ions. The planewave-based field gradient is computed in reciprocal space and then evaluated at the nuclear site of interest by Fourier-transformation; the PAW-sphere corrections are computed in real space, noting that the use of the PAW compensation charge ensures that no inter-sphere contributions need be considered; and the gradient due to the fixed ions is computed using an Ewald-summation method.^{46} As an example of the application of this feature, we showed^{41} the sensitivity of the electric field gradient response to strong correlation effects in LaTiO_{3} by combining the field gradient calculation with the DFT+*U* method, also available in abinit (see Sec. IV A).

### D. Positron annihilation

One of the unique features of abinit is the fully self-consistent implementation of two-component density functional theory (TCDFT)^{47,48} within the PAW formalism.^{44} This technique allows one to accurately compute various properties of annihilating electron–positron pairs, which, in turn, can be used to interpret positron annihilation spectroscopy (PAS) measurements.^{49}

Within TCDFT, the total energy of interacting positrons and electrons is written as

where *n*^{+} and *n*^{−} are positron and electron densities, *E*[*n*^{+}] and *E*[*n*^{−}] are one-component functionals for positron and electrons, and $\u222bdr\u2009\u222bdr\u2032\u2009n\u2212(r)n+(r)|r\u2212r\u2032|$ corresponds to the Hartree interaction. Various approximations can be used to calculate the $Ece\u2212p$ term, both within LDA^{48,50} and GGA.^{51,52}

The TCDFT implementation in abinit is based on a unified formalism in which both electron and positron wavefunctions are expressed using the same mixed PAW basis (planewaves and atomic orbitals).^{53} This means that the electronic and positronic energies and forces can be calculated self-consistently. For example, it allows for full geometry optimization of systems containing positrons, an effect that has been shown to be critical in determining reliable annihilation features.^{54} However, it is worth noting that the positron calculations within the PAW method are sensitive to the completeness of the PAW dataset. To achieve an accurate description of positron densities and wavefunctions, it is often necessary to include semicore electrons in the PAW dataset.^{53}

Various properties of annihilating electron–positron pairs can be calculated within abinit. First, positronic wavefunctions and densities in the direct space can be accessed and visualized to inspect the localization of the particle in a given system. An example is given in Fig. 4, where the isodensity of a positron localized in a silicon vacancy in SiC is shown. Second, the total energies of a system containing a positron can be used to calculate binding energies or affinities. Third, based on the electron and positron densities, the positron lifetime *τ* can be calculated as the inverse of the annihilation rate *λ*,

where *g*(*n*^{−}, *n*^{+}) is an enhancement factor, corresponding to the increase in annihilation due to the screening of the positron by the electrons. Several parameterizations for *g*(*n*^{−}, *n*^{+}) can be used in abinit. Finally, the electron and positron wavefunctions can be used to calculate momentum distributions *ρ*(**p**) of annihilating pairs, following^{47}

where $\Psi ie-p$ is the two-particle wavefunction in the state *i* and **p** is the given momentum.

The TCDFT implementation in abinit has already been applied to various systems. It has been used to calculate the annihilation properties of bulk metals and semiconductors.^{53} Calculations for vacancy-type defects in systems, such as SiC,^{54–56} UO_{2},^{57,58} GaN,^{59} and yttria-stabilized zirconia,^{60} have been performed and used to analyze experimental data. Positron lifetimes and binding energies have also been investigated on Fe(001) surfaces with and without adatoms.^{61}

## V. RESPONSE PROPERTIES: SPECTROSCOPY, VIBRATIONS, DIELECTRIC RESPONSE, TEMPERATURE DEPENDENCE, AND SPIN-MAGNETIC FIELD COUPLING

Many physical properties of materials can be formulated as derivatives of the energy. Although such derivatives might be computed from finite differences techniques, they can also be conveniently determined from perturbation theory. abinit does not only compute many energy derivatives from DFPT but also gives readily access to various related physical quantities in Cartesian coordinates and conventional units.

As a basic feature, abinit routinely determines standard second energy derivatives with respect to atomic displacement (*τ*^{q}), homogeneous static electric field ($E$ ), and homogeneous strain (*η*) perturbations,^{33,34,63} which are related to dynamical matrices at any **q**-vector (*D*^{q}), optical dielectric tensor (*ε*^{∞}), Born effective charge tensors (**Z**^{*}), clamped-ion elastic tensors (*C*^{0}), clamped-ion piezoelectric tensors (*e*^{0}), and lattice-strain coupling tensor (*γ*), as summarized in Table IV. Relying on these bare quantities, abinit also provides access to additional properties such as the full phonon dispersion curves (interpolated by separating short-range and long-range dipolar interactions), interatomic force constants in real space (short-range and dipolar contributions), infrared intensities, static and infrared dielectric constants, and relaxed ion elastic and piezoelectric tensors. All of these are accessible both with NCPSP and in PAW, using LDA or GGA functionals. The specific case of the strain perturbation^{63} that is a unique feature of abinit is further discussed in Sec. V B.

. | $\u2202\u2202\tau $ . | $\u2202\u2202E$ . | $\u2202\u2202\eta $ . |
---|---|---|---|

$\u2202\u2202\tau $ | D^{q} | Z^{*} | γ |

$\u2202\u2202E$ | Z^{*} | ε^{∞} | e^{0} |

$\u2202\u2202\eta $ | γ | e^{0} | C^{0} |

. | $\u2202\u2202\tau $ . | $\u2202\u2202E$ . | $\u2202\u2202\eta $ . |
---|---|---|---|

$\u2202\u2202\tau $ | D^{q} | Z^{*} | γ |

$\u2202\u2202E$ | Z^{*} | ε^{∞} | e^{0} |

$\u2202\u2202\eta $ | γ | e^{0} | C^{0} |

Going further, abinit also implements some third-order energy derivatives and responses to additional perturbations, providing access to an even much broader set of properties, as discussed below.

### A. Raman spectroscopy and the electro-optic effect

Beyond the development of DFPT up to second-order,^{34} non-linear properties can be accessed using third-order DFPT. Thanks to the 2*n* + 1 theorem,^{9,64,65} third-order energy derivatives can be obtained from wavefunction derivatives up to first order, as already done in second-order response calculations. Therefore, the additional computational cost to access such third energy derivatives is typically negligible compared to that needed for second-order energy derivatives. There are, however, some linear and non-linear properties in which an analytic treatment of the electric field perturbation requires supplementary wavefunction derivatives.^{33}

In abinit, the accessible non-linear properties are (i) the second-order optical susceptibilities (*χ*^{(2)}), (ii) the Raman tensors (*α*^{m}), and (iii) the electro-optic tensor (*r*).

The second-order nonlinear optical susceptibility is a third-rank tensor related to the response of the electrons of the system to optical electric fields, which are *a priori* frequency dependent. In the present abinit context of the 2*n* + 1 theorem applied to (static) DFT, we neglect the dispersion of *χ*^{(2)} and compute it as the purely electronic response of the system at zero frequency. Within these conditions, *χ*^{(2)} is formulated as a third-order derivative of a field-dependent energy functional,

The *χ*^{(2)} tensor delivered by abinit is related to the non-linear optic “d-tensor” as $dijk=(1/2)\chi ijk(2)$ and can be calculated with an additional scissor correction on the electron band energies. Such a second-order optical susceptibility has been recently used for high-throughput second-harmonic generation calculations (see Sec. VI B 2).

The Raman tensor describes the change in linear optical susceptibility, *χ*_{ij}, produced by an atomic displacement *τ*_{κβ} of atom *κ* in direction *β*. It is, therefore, related to the following third energy derivative:^{66}

Using Placzek’s approximation,^{67} the Raman intensity $\alpha ijm$ of a Stokes process associated with a phonon mode *m* depends on the Raman tensors as^{66,68,69}

Depending on the type of phonon mode, i.e., transverse optic (TO) or longitudinal optic (LO), the derivative of the electric susceptibility with respect to an atomic displacement must be written as^{66,70}

where **q** is the direction along which the Brillouin zone center is approached. Both quantities are provided by abinit.

The electro-optic tensor is related to the change in linear optical susceptibility produced by a static field. As such, it combines *a priori* electronic, lattice, and strain contributions. The so-called clamped (fixed strain) electro-optic tensor combining electronic and lattice contributions is directly provided by abinit from the knowledge of previously discussed quantities,^{66}

The unclamped electro-optic tensor includes an additional piezoelectric contribution (strain response) that can be computed but requires an independent calculation of the elasto-optic coefficients by finite differences.^{71}

The computation of the third-order energy derivatives in Eqs. (7) and (8) (providing access to the different tensors discussed above) was first implemented in abinit relying on the “Berry phase” formalism;^{66} this implementation is available using a NCPSP approach within the LDA and for unpolarized systems. Following the scheme of Gonze,^{33} Miwa^{72} used an alternative approach, where the electric field perturbation is treated analytically, leading to a derivative operator with respect to **k** vectors ($\u2202/\u2202k\alpha $). At the non-linear level, this derivative operator is applied to first-order derivatives of wavefunctions, leading to second-order wavefunction derivatives. This way, even in the context of the 2*n* + 1 theorem, a (non-self-consistent) second-order Sternheimer equation has to be solved. Our second implementation follows that the scheme is valid for both NCPSP and PAW approaches. It is available with LDA functionals and is valid for unpolarized or collinear spin-polarized systems. This second implementation also converges faster than the original one in terms of the k-point sampling. The extension to the PAW+*U* formalism is under development.

abinit delivers all the tensors given above in Cartesian coordinates and with clear units (SI or conventional literature choices). These quantities can then be combined to compute the polarization dependent Raman spectra of single crystals or even the Raman spectra of powders. In either case, the Raman intensities can be obtained as^{74}

where *G*_{0}, *G*_{1}, and *G*_{2} are functions of the components of Eq. (9)^{75} and a Lorentzian function is used for the broadening,

Post-processing scripts are provided with abinit to extract the Raman tensor and print the powder-average and polarization dependent Raman intensity to a file. An example of powder-averaged Raman spectra for polycrystalline quartz is shown in Fig. 5.

Using the non-linear tensors computed by abinit, the theoretical Raman spectra have been successfully compared with the experiment for many materials, including oxides,^{76–78} organic compounds,^{79} and transition metal dichalcogenides^{80,81} (see also the WURM database of Raman spectra^{82}). Electro-optic tensors have also been computed for various compounds, including ferroelectric oxides,^{71,78} multiferroics,^{83} and monolayer transition metal dichalcogenides.^{81} Going further, relying on these quantities obtained in DFT at zero Kelvin, non-linear optical properties at finite temperature were also successfully computed, making use of a first-principles based effective Hamiltonian approach.^{84}

### B. The strain perturbation in density-functional perturbation theory

One of the unique capabilities of abinit is the direct calculation of the elastic and piezoelectric properties of materials by DFPT. This involves calculating the responses to three types of perturbations: atomic displacements, electric fields, and strains and appropriately combining the results.^{85} The general structure of DFPT is based on the systematic expansion of the variational expression for the DFT total energy in powers of a parameter *λ* characterizing some dependence of the energy functional. Such parameters as internal atomic coordinates and the macroscopic electric field could be handled in this framework in a conceptually straightforward manner. Treating macroscopic strain as a parameter within this formalism, however, is less straightforward. The approach developed for abinit was based on its existing overall formulation of the DFT energy expression in reduced coordinates. This introduces real- and reciprocal-space metric tensors into every expression. The underlying reduced-coordinate lattice structure is invariant—all unit cells are unit cubes. Strain affects only the metric tensors, and the *λ* derivatives of every term in the Hamiltonian can be developed from derivatives of these tensors.^{63}

The second derivatives of the energy with respect to a pair of perturbations are evaluated using the so-called non-stationary expressions,

where $\psi \alpha (\lambda 2)$ is first-order wavefunction obtained from the self-consistent Sternheimer equation for perturbation *λ*_{2}, the parenthesized superscripts denote partial derivatives, and the notation for the various energy terms and the ground-state wavefunctions should be self-explanatory.^{33} The derivative with respect to strain and the electric field, necessary for the piezoelectric tensor, has the alternative form

where $\psi km(k\u0303j)$ is the first-order wavefunction for the *∂*/*∂k* perturbation.^{86} The metric tensors involved in this context are simple and constant throughout space. They are expressed in terms of the real- and reciprocal-space primitive lattice vectors as

Dot products of real- and reciprocal-space vectors can be expressed in terms of their strain-invariant reduced-space counterparts using these tensors. For example, reciprocal-space vectors **K** = **k** + **G** represented by their reduced counterpart $K\u0303$ have dot products,

First and second strain derivatives of the metric tensors follow from the derivatives of real and reciprocal space vectors with respect to strains,

The only simplifying aspect of all this is that the dot product of a real and a reciprocal lattice vector, typically appearing in phase factors, is merely 2*π* times its reduced counterparts, and therefore strain-independent.^{63}

All the terms in the total energy expression in reduced coordinates can be expressed in terms of dot products. Therefore, the evaluation of Eq. (15) is, in principle, straightforward, while in practice exceedingly tedious. This formalism was initially developed for multi-projector NCPSP and LDA exchange–correlation functionals. There are special considerations for GGA functionals that were added.^{87} Extension to PAW potentials within LDA was recently added.^{88}

The terms described so far, derivatives with respect to two strain components (the elastic tensor), one strain component and a uniform electric field (the piezoelectric tensor) and one strain component and one atomic displacement (the so-called internal strain) are all based on strain-independent atomic positions in reduced coordinates. In fact, infinitesimal strains do produce infinitesimal reduced-coordinate changes in atomic positions. Taking this into account modifies the clamped-atom elastic and piezoelectric tensors calculated so far.

Combining the clamped-ion derivatives with the atomic position contributions is reminiscent of the chain rule for derivatives.^{85} The necessary intermediate derivatives (interatomic force constants and Born effective charges) were already available in abinit. The DFPT calculations of the clamped-ion tensors and these additional derivatives are performed separately in abinit, and the stored results are combined in the auxiliary program ANADDB to produce the final, physical relaxed-atom quantities. Both the NCPSP and PAW formalism and codes have been verified to a high degree of accuracy, by comparison with finite-difference calculations.

### C. Density-functional perturbation theory including van der Waals interaction: Phonons and strain perturbation

Most DFT functionals lack a proper treatment of the long-range e^{−}–e^{−} correlation. This term is crucial to properly describe the structure of weakly bound compounds and, consequently, their underlying properties (phonons, elastic constants, etc.). Several methods have been developed in the past in order to include some long-range e^{−}–e^{−} correlation in DFT computations. Notably, one could mention Grimme’s DFT-D methods,^{89,90} the vdW-DF functionals,^{91,92} as well as the TS-vdW and MBD methods.^{93,94} All of these approaches add new contributions to the energy and its derivatives that have to be properly taken into account in DFPT computations for self-consistency.

abinit allows one to carry out such computations for Grimme’s DFT-D2, -D3 and -D3(BJ) methods.^{95,96} These methods add to the DFT energy a contribution that depends only on the atomic positions,

where *I* and *J* are atomic indices, *C*_{6,IJ} is the dispersion coefficient that depends on the whole set of atomic positions **R**, and *f*(**R**_{IJ}) is a function that depends on the chosen Grimme flavor ($\u223c\u2009RIJ\u22126$ at long range) and on **R**_{IJ}, the difference in atomic position between the *I* and *J* atoms. This dispersion contribution is not compatible with all *ad hoc* exchange–correlation functionals (see Refs. 89 and 90 for the current list of compatible functionals). The common GGA-based functionals, such as GGA-PBE (Perdew-Burke-Ernzerhof)^{208} or revPBE,^{209} are nonetheless available. Due to the way these DFT-D methods were constructed, they may not be the most suited to deal with metals.

Nearly, all response functions available in abinit can be computed in DFPT with the inclusion of DFT-D dispersion corrections, with the only exception of flexoelectricity. Note that since these corrections do not rely on the knowledge of the electronic density, the (total) energy derivatives with respect to electric fields remain untouched by the dispersion scheme. This is both an advantage and a drawback of the DFT-D method, since it eases considerably its implementation for response functions, but it prevents at the same time the method to capture any direct electronic contribution (dielectric constants, Raman susceptibility, etc.). The lack of literature on this topic does not allow us to rule on the importance of dispersion contributions in the proper description of these properties. The reason why the DFT-D dispersion corrections are not available for flexoelectricity comes from the indirect dependence of this property on the spatial derivative of the dynamical matrix, as mentioned in Sec. V D. Such spatial derivative is not yet implemented.

All response functions related to the energy derivatives with respect to atomic displacements or to strains are directly affected by the dispersion corrections. Thanks to their pair-wise form, the DFT-D strain-related quantities can be related to derivatives with respect to atomic positions,

where *κ* is an atomic index.

The investigation of several weakly bound systems (benzene, MoS_{2}, etc.) using DFT-D methods has already highlighted the importance of the dispersion contributions to interatomic force constants and to phonons.^{80,95} Such methods generally lead to a better agreement with the experiments than ordinary GGA-PBE. Similarly to the general DFPT implementation, not only the zone-center phonons but also those associated with an arbitrary wavevector are accessible, allowing the computations of phonon band structures in a fully-consistent manner. In addition, the present implementation can also be straightforwardly used to study thermal expansion within the quasi-harmonic approximation (see Ref. 97 for an example application) or electron–phonon coupling.

### D. Spatial dispersion: Flexoelectricity and dynamical quadrupoles

The flexoelectric (FxE) effect, where a strain gradient deformation results in a macroscopic polarization, is a challenging electromechanical coupling to simulate computationally. The main reason is the *a priori* incompatibility of a spatially varying strain perturbation with periodic boundary conditions. A recent set of developments in abinit has overcome this difficulty by adapting the long-wave method introduced by Born and Huang^{98} in the early days of quantum mechanics to the modern tools of the DFPT. abinit now has the unique capability to calculate the bulk FxE tensor for any insulating material.

Achieving this result has taken nearly a decade of continuing efforts, in order to settle the remarkable number of formal subtleties. Based on the seminal work of Resta,^{99} early attempts^{100,101} made use of supercell geometries to calculate the relevant ingredients discussed below. Later, a long-wave framework in a reciprocal space^{102} was established, together with a curvilinear-coordinates formulation,^{103} which led to the calculation of the full flexoelectric response of SrTiO_{3} slabs^{104} (see Ref. 105 for a summary of the pre-2015 developments). En route toward a practical implementation, additional technical and formal issues were addressed, regarding the generalization of the uniform strain response to a “metric-wave” perturbation^{106,107} and the proper definition of the current-density operator in the presence of nonlocal pseudopotentials.^{108} These efforts culminated with the present long-wave method,^{109} which paves the way toward the computation of many spatial dispersion properties with a computational burden that is comparable to conventional linear-response calculations.

An *ab initio* calculation of the FxE tensor (** μ**) requires the sum of an electronic (clamped-ion), a lattice-mediated, and a mixed contribution. This intricate structure can be formulated as

^{102,105}

where *κ* labels the atomic sublattice and the II superscript indicates that the strain gradient tensor is defined as the gradient of the symmetric strain, i.e., the so-called type-II definition is assumed. The quantities on the right-hand side are the following: $\mu \xaf\alpha \lambda ,\beta \gamma II$, the clamped-ion FxE tensor; $P\xaf\alpha ,\kappa \rho (1,\lambda )$, the first spatial moment of the polarization field induced by an atomic displacement (also known as **P**^{(1)} tensor); $\Gamma \rho \beta \gamma \kappa $, the piezoelectric internal strain tensor; and $L\rho \lambda ,\beta \gamma \kappa $, the flexoelectric internal strain tensor.

However, **Z**^{*} and **Γ** are well-known quantities, which can be obtained by means of the linear response DFPT with *standard* perturbations such as atomic displacements, electric field and strain;^{33,34,63} Eq. (22) demands the calculation of up to four new tensor quantities. These include the $\mu \xaf\alpha \lambda ,\beta \gamma II$ and $P\xaf\alpha ,\kappa \rho (1,\lambda )$ tensors introduced above, but also the first spatial moments of the dynamical matrix ($\Phi \kappa \alpha ,\kappa \u2032\rho (1,\lambda )$) and of the piezoelectric force-response tensor ($C\xaf\alpha \lambda ,\beta \gamma \kappa $), both required to build $L\rho \lambda ,\beta \gamma \kappa $ following the prescriptions in Ref. 102. These new tensors represent four spatial dispersion properties, i.e., they are only sensitive to a perturbation *gradient*, not to a uniform perturbation. In principle, however, it is possible to obtain them as momentum derivatives of counterpart tensors from uniform perturbations, when the corresponding DFPT second order energies are generalized to finite momentum **q**. Thus, for example, the clamped-ion FxE tensor is the spatial dispersion counterpart of the clamped-ion piezoelectric tensor (*ē*_{α,βγ})^{63} and could be calculated as the momentum derivative of a finite-**q** version of *ē*_{α,βγ}.

This is precisely the procedure adopted in the long-wave DFPT of Ref. 109 that has been implemented in abinit. In order to obtain the finite-**q** generalization of the second order energies, the approach reformulates the electric-field and the strain perturbation problems as the time derivative of a vector potential and as the gradient of a metric wave,^{107} respectively. In this way, both perturbations are generalized to finite **q**, as is already the case for atomic displacements. This enables us to carry out an analytical third order derivative of the energy with respect to two of the standard perturbations and to the momentum **q**, which directly provides the sought-after spatial dispersion tensors. Remarkably, by virtue of the 2*n* + 1 theorem,^{9} the third-order energies are computed in one shot using precalculated first-order response functions to the standard perturbations, without the necessity of self-consistently computing any response function to a perturbation gradient. After execution, the long-wave DFPT routines generate a derivative database that is subsequently used by post-processing tools implemented in ANADDB to compute and print the different contributions to the FxE tensor.

The implementation also provides access to another spatial dispersion property, the dynamical quadrupoles, which can be obtained as the symmetrized sum of the **P**^{(1)} tensor,^{109}

The dynamical quadrupoles are the spatial dispersion counterparts of **Z**^{*} and can be used in lattice dynamics calculations to improve the prevalent dipole–dipole treatment of the long-range interactions. The anaddb routines that carry out the process of interpolating the dynamical matrix following Ref. 34 have been adapted to incorporate the dipole–quadrupole and quadrupole–quadrupole electrostatic interactions derived in Ref. 102. This new functionality results in a faster convergence of the phonon band calculation with respect to the density of **q** points and, in some materials, represents the only route to obtain the correct sound velocities.

The abinit implementation of the long-wave DFPT has been used to compute the clamped-ion FxE tensor of Si and SrTiO_{3}, obtaining an excellent agreement with the existing results in the literature.^{109} The computed spatial dispersion properties accurately reproduce the sum rules predicted in Refs. 110 and 102 that relate them with specific quantities from uniform perturbation theories.

Currently, the implementation is restricted to the use of NCPSP without non-linear core corrections and the LDA functional.

### E. Electron–phonon interaction

With abinit, it is possible to compute many physical properties related to electron–phonon (e–ph) interaction. In metallic systems, for instance, one can study conventional superconducting properties within the isotropic Migdal–Eliashberg formalism^{111} and compute transport properties in the normal state by solving the linearized Boltzmann equation within the LOVA approach.^{112–114} However, in this section, we focus on temperature-dependent band structures and the zero-point renormalization of the bandgap in semiconductors,^{115,116} a subject that in recent years has received increased attention within the electronic structure community.^{117–120} As discussed in the review paper by Giustino,^{111} the renormalization of the electron state due to the e–ph interaction is described by the self-energy Σ^{e-ph}(*ω*) = Σ^{FM}(*ω*) + Σ^{DW}. The diagonal matrix elements of the Fan–Migdal self-energy in the KS basis set are defined by

where *η* is a positive real infinitesimal. The e–ph matrix elements *g*_{mnν}(**k**, **q**) are given by

with Δ_{qν}*V*^{KS} the first-order variation of the self-consistent KS potential that can be computed with DFPT.^{33,121} The static Debye–Waller (DW) term involves the second order derivative of the KS potential with respect to the nuclear displacements. State-of-the-art implementations approximate the DW contribution with

where $gmn\nu 2,DW(k,q)$ is an effective matrix element that, within the rigid-ion approximation, can be expressed in terms of the *g*_{mnν}(**k**, **q**) matrix elements.^{111} In principle, the quasi-particle (QP) excitations are defined by the solution(s) in the complex plane of the following equation: $z=\epsilon nk+\Sigma nke-ph(z)$. In practice, the problem is usually simplified by seeking approximated solutions along the real axis following two different approaches. In the on-the-mass-shell approximation, the QP energy is given by the real part of the self-energy evaluated at the bare KS eigenvalue,

The second approach linearizes the self-energy near the KS eigenvalue and evaluates the QP correction using

with the renormalization factor *Z* given by

Both approaches are implemented in abinit although it should be noted that, according to recent works, the on-the-mass-shell approach provides results that are closer to those obtained with more advanced techniques based on the cumulant expansion^{122} or self-energy calculations employing an eigenvalue–self-consistent cycle.^{123}

Accurate calculations of e–ph renormalization are still challenging even by present standards because the e–ph self-energy is quite sensitive to the **q**-point sampling. Moreover, a large number of empty states *m* are usually required to converge the real part of the self-energy. In order to avoid the explicit computation of the DFPT scattering potentials on dense **q**-grids, abinit generalizes the Fourier-based interpolation scheme proposed by Eiguren *et al.*^{124} to account for the non-analytical behavior associated with the long-range interactions present in semiconductors.^{125} The dipolar potentials generated by the Born effective charges in polar materials are treated using a generalized Fröhlich model.^{126,127} Furthermore, details concerning the implementation are given in Ref. 17. In order to accelerate the convergence with the number of empty states, abinit replaces the contributions given by the high-energy states above a certain band index *M* with the solution of a non-self-consistent Sternheimer equation in which only the first *M* states are required. The methodology, proposed in Ref. 128, is based on a quasi-static approximation in which the phonon frequencies in the denominator of Eq. (24) are neglected, and the frequency dependence of Σ is approximated with the value computed at *ω* = *ε*_{nk}. This approximation is justified when the bands above *M* are sufficiently high in energy with respect to the *n***k** states that must be corrected. Furthermore, this upper-band contribution to Σ converges quickly with respect to **q**-point sampling, and it can be safely computed on a coarse **q**-grid.^{123,129}

The code can compute QP corrections and lifetimes due to e–ph scattering as well as spectral functions. The lifetimes obtained from the imaginary part of Eq. (24) can be used to compute carrier mobilities within the self-energy relaxation time approximation,^{111,130,131}

For the computation of electron lifetimes, abinit implements advanced integration techniques that take the advantage of the linear tetrahedron method^{44} and double-grid techniques for the **q**-space integration. Last but not least, abinit provides a specialized driver to interpolate the DFPT scattering potentials, compute e–ph matrix elements on arbitrarily dense **k**-meshes, and save the results to netcdf files. This driver can be used by external codes such as berkeleygw^{132} to treat e–ph interactions at the *GW* level, as discussed in Ref. 133.

### F. a-tdep: Temperature dependent thermodynamic properties using the TDEP approach

The features and the temperature dependence of the phonon spectrum determine a large number of thermodynamic properties of crystals, such as dynamic and thermodynamic stability, elastic properties, and heat transport within the material. The theory of lattice dynamics at finite temperature has long been a central field of research in condensed matter, and its applications are numerous in geophysics, material science, astrophysics, etc.

The calculation of the phonon spectrum is available in abinit since the 1990s and is based on DFPT. In this framework, the ground state undergoes a small perturbation around its 0 K equilibrium positions in order to describe the shape of the potential energy surface (PES) very near the minimum. Thereafter, in order to capture the temperature effects, the harmonic approximation (HA) or quasi-harmonic approximation (QHA) can be applied. The latter approximation is handled by assuming that the temperature dependency of the phonon frequencies can be taken into account implicitly through a variation of the volume. This treatment generally gives excellent agreement with experimental measurements of the thermal expansion and the elastic constants.^{134}

However, the QHA fails or cannot be applied in various particular cases. For example, when the system is close to a phase transition, when the phase of interest is not stable at 0 K, or when one wants properties at high temperature. In these situations, an implicit treatment of the temperature through a volume variation is no longer sufficient and a dedicated calculation, taking into account the temperature explicitly, is needed. One can formalize and synthesize this statement as follows:

The second term of the right-hand side of Eq. (31) is included in the QHA. However, the first term on the right-hand side of Eq. (31) is only included when the temperature dependence is explicit.

Since the beginning of the 1960s, a large number of theoretical studies have been carried out in order to go beyond the simple harmonic crystal and to deal with anharmonic effects.^{135–141} More recently, there has been a revival of such developments. During the last decade, several groups proposed that theoretical frameworks coupled to computational methods are capable of capturing anharmonic effects in crystals.^{142–148} One way, chosen and implemented in abinit, has been proposed by Esfarjani and Stokes^{149} and then developed by Hellman *et al.*^{150–152} In the framework of temperature-dependent effective potential (TDEP), the anharmonic terms are treated in an effective manner.

Let us define a 3-dimensional crystal and consider that a ground state energy *U*_{0} is obtained when the atoms are in their equilibrium positions. The potential energy of this system can be rewritten using a Taylor expansion around equilibrium positions as

with $\Phi (p)i1\u2026ip\alpha 1\u2026\alpha p$ being the interatomic force constants (IFC) at the *p*th order and $ui\alpha $ being the displacement of atom *i* along Cartesian direction *α*.

By performing *ab initio* molecular dynamics (MD) simulations, the forces $Fi,MD\alpha (t)$ and atomic displacements $ui,MD\alpha (t)$ are acquired at each time step *t*. If we differentiate Eq. (32) with respect to atomic displacements and insert the MD quantities, we obtain the following system of equations:

with $\Theta (p)ij\u2026l\alpha \beta \u2026\delta $ being the effective IFC at the *p*th order and *P* being the maximum order of the expansion. Equation (33) is non-linear with respect to the atomic displacements and can be rewritten as

where *θ*_{pλ} being the *λ*^{th} are now the *effective* IFC at the *p*th order and $fi,p\lambda \alpha (uMD(t))$ is a function gathering all the contributions coming from the atomic displacements. This system of equations is now linear as a function of *θ*_{pλ} (the unknown variables of the system) and can be solved using a least-squares method. The solution is given by **Θ** = **f**^{†}·**F**_{MD}, with $FMD\u2261Fi,MD\alpha (t)$, **Θ** ≡ *θ*_{pλ}, and **f**^{†} being the pseudoinverse of $f\u2261fi,p\lambda \alpha (uMD(t))$.

Note the difference between the true $\Phi (p)$ and effective $\Theta (p)$ IFC. The latter includes the effects of all the IFC not taken into account in the expansion (i.e., beyond the *P*th order). For instance, if *P* = 2, the 2nd order effective IFC incorporates all higher anharmonic terms in an effective way. Consequently, the effective IFC acquire a dependency on the temperature at odds with the true IFC.

The number of time steps needed to solve Eq. (34) has to be large with respect to the number of unknown variables *θ*_{pλ}. By taking into account the properties of the system (translation and rotation invariances, crystal symmetries, etc.) the number of independent and non-zero IFC coefficients of a cubic system with one hundred atoms in the supercell can be reduced to around ten at the second order, tens at the third order, and around a hundred at the fourth order.

Once the IFC coefficients are obtained, many dynamic, elastic, and thermodynamic properties can be evaluated:

As usual, the phonon spectrum is obtained after the diagonalization of the dynamical matrix.

Then, the phonon density of states is built and used to compute various thermodynamic properties:

^{153}the free energy, the specific heat, the entropy, etc.By computing the elastic constants using the second-order effective IFC,

^{135,141}it is possible to obtain the bulk and shear moduli.The Grüneisen parameter is obtained using the third order effective IFC,

^{154,155}which leads to other important quantities: the thermal expansion, the isentropic compressibility, the constant pressure specific heat, sound velocities, etc.

Since the effective IFC include an explicit temperature dependency, all the quantities listed above do as well.

The abinit implementation of the TDEP algorithm, referred to as a-tdep, has produced a number of applications in the last five years.^{156–160} We recommend the user to read the a-tdep documentation distributed in the abinit package before any calculations and also to read the article dedicated to this implementation.^{161}

### G. Density functional perturbation theory within non-collinear magnetism

The necessity to treat non-collinear magnetism (NCM) in DFT codes is more and more important for spintronic applications studies. NCM can arise naturally, from geometrical frustration of anti-ferromagnetic interactions, e.g., in triangular lattices, by the magnetic anisotropy generated in the presence of a preferred direction of magnetization or by the competition between exchange interactions. However, the generalization of DFPT for non-collinear magnetic simulations has only recently been attempted. Within the non-collinear DFT scheme, the density is described by a four-component matrix instead of two (in collinear magnetism), which is required to account for the spin–orbit coupling (SOC), except when time-reversal symmetry induces vanishing of the magnetization everywhere.^{162} This last term, in current DFT code implementations, implies a partial loss of symmetry,^{163–165} which worsens the calculation time and explains the scarcity of such calculations. Additionally, the treatment of the electronic exchange–correlation (XC) term is more complicated, since the functionals are usually built in a collinear framework. We have proposed^{166} three different approaches to perform the energy derivatives in the non-collinear regime treating opportunely the XC term and have implemented these for the atomic displacement (phonons) and electric field perturbations.

The first method exploits the local transformation from non-collinear to collinear magnetization: at each point of space, we align the local system of reference to the local quantization axis of the magnetization, then we perform the derivatives in this locally collinear framework, and finally we return to the lab frame of reference. The unitary transformation diagonalizing the 2 × 2 density matrix $\rho ^(0)$, and its derivative $\rho ^(1)$, has to be determined at the zeroth (*U*^{(0)}) and first order (*U*^{(1)}) of perturbation to fully characterize the problem (see Ref. 166).

The second method consists in writing the *U*^{(0)} and *U*^{(1)} analytically in terms of the Pauli matrices *σ*_{αβ}, i.e., the generators of the Lie group SU(2).

The third method evaluates explicitly the expression of the first order XC potential, specifically for the LSDA expression of the XC energy, which is easy to derive.

The three methods have been tested on two particular systems Cr_{2}O_{3} and RuCl_{3}.^{166} Cr_{2}O_{3} has been used to test the implementations because it is one of the simplest and best known collinear antiferromagnetic and magneto-electric systems. Thanks to these tests, we have demonstrated that the two first methods are closely equivalent in terms of second derivatives: quality and timing performance (calculations of phonons, electronic dielectric constant, and Born effective charges). However, the third one requires twice as many self-consistent iterations to give the same quantities. On RuCl_{3}, we have applied our methods to predict the change in phonon frequencies induced by the non-collinearity of the magnetic moments. Comparing the frequencies estimated in collinear and non-collinear frameworks, we have shown the signature of non-negligible NCM: certain phonon modes of specific irreducible representations have their frequency modified by up to 10 cm^{−1} depending on the chosen direction (see Ref. 166 for more details).

The different methods of magnetization rotation can be set through the **ixcrot** abinit input flag. The default value is method 1 (ixcrot = 1). Methods 2 or 3 can be used by setting **ixcrot** to 2 or 3.^{167} Although in the corresponding abinit version of Ref. 166, we implemented the NCM DFPT formalism at the Γ point only, further developments have been completed in order to account for the full **q** dependence in the phonon spectra (available in the next production version, only with **ixcrot** = 3). The implementation is valid for NCPSP only, but work is ongoing to extend it to the PAW formalism.

### H. Spin-magnetic-field perturbation

The linear response of the magnetization to static and dynamical magnetic fields (spin magnetic susceptibility *χ*_{m}) is one of the key characteristics of magnetic materials. This quantity is intrinsically linked to the spin fluctuation spectrum and hence not only allows one to identify and characterize the relevant magnetic excitations but also to understand their character and estimate the strength of spin interactions. The latter case appears to be particularly important in the view of the ever-growing need for accurate lattice spin models to predict the finite-temperature properties of spin systems and understand the physics of magnetic and topological materials. While the theoretical framework for linear spin response calculations within the DFPT formalism has been established several decades ago,^{168} practical implementations of this method are still rare in DFT codes.^{169}

We have implemented and tested the variational DFPT method for computing the linear spin response to a small, q-dependent, magnetic field of the form **B** = **b***e*^{iq⋅r} + **b**^{*}*e*^{−iq⋅r}. The first derivatives of the density and magnetization are obtained by considering the first order external potential $vext(1)=12\sigma \alpha $, where the Cartesian component *α* of the Pauli matrix corresponds to the direction of the external field. The periodic part of the first order (spin) magnetization is represented as $mq(1)=\u2211nkunk\u2020\sigma unk+q(1)\u2009+\u2009unk\u2212q\u2020\u2009(1)\sigma unk$, with the sum running over the occupied Kohn–Sham orbitals *n*. For an arbitrary wavevector **q** ≠ 0, it is, thus, necessary to find the Bloch-periodic parts of the first-order wavefunctions $unk\u2212q(1)$ and $unk+q(1)$ at both **k** + **q** and **k** − **q** points. This is done by performing two state-by-state minimizations of the second order energy functional at each spectral correction factor iteration, using the Self-Consistent Field (SCF) potential $vscf,q(1)$ and its Hermitian conjugate $vscf,\u2212q(1)$, respectively.

The implementation of the algorithm has been tested by comparing the spin susceptibility obtained from DFPT calculations with the finite-difference derivatives estimated from ground state DFT calculations with a finite magnetic field. Specifically, for **q** = 0, we have performed calculations of the longitudinal magnetization response of ferromagnetic bcc Fe as well as longitudinal and transverse spin susceptibilities of antiferromagnetic Cr_{2}O_{3}. These two examples were chosen to test both metallic and insulating occupations. To verify the validity of **q** ≠ 0 implementation, we have further performed ground-state DFT calculations with harmonically varying Zeeman field *B* = **b** cos(**qr**) with **q** = (1/2, 0, 0) and **q** = (1/4, 0, 0) for bcc Fe in the corresponding supercells. The periodic parts of the first order density and magnetization obtained using the first-order finite-difference scheme were then compared to the corresponding results obtained from a single unit cell DFPT run. All considered test cases for **q** = 0 and **q** ≠ 0 showed good agreement between DFPT and finite field results, even in the case of coarse Brillouin zone sampling, e.g., 2 × 2 × 2 and 4 × 4 × 4 k-point grids. Furthermore, to verify the 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 the corresponding Fourier components of the super-exchange parameters *J*(**q**) = *χ*(**q**)^{(−1)}, yielding couplings of 18.9 meV (9.8 meV) for the first (second) shell of nearest neighbors, in good agreement with previous literature.^{170}

All the test simulations described above were performed within the LSDA approximation, which is currently the only exchange–correlation functional implemented for linear spin response calculations in abinit. Furthermore, for the **q** ≠ 0, case, it is imperative to use the explicit evaluation of the first-order XC functional (ixcrot = 3 abinit input flag). Furthermore, at present and similarly to atomic displacement perturbations with non-collinear magnetism, the implementation of Zeeman magnetic field perturbation is limited to NCPSP formalism.

### I. Optical response with the Bethe–Salpeter equation, including temperature dependence

Within MBPT, one can address the computation of the energies and related characteristics of charged excitations as well as those of neutral excitations.^{11} For the former, the *GW* approximation is the state-of-the-art, while for the latter, one relies on the Bethe–Salpeter equation (BSE), which includes excitonic effects. Such approaches produce optical absorption spectra that are much more accurate than independent-particle approaches, such as the simple sum-over-states.

The implementation of the BSE in abinit has been previously described in Refs. 5 and 171. An efficient interpolation allows one to cut down significantly the Central Processing Unit (CPU) time needed to obtain optical spectra converged with respect to the Brillouin zone sampling.

As a salient example, the abinit implementation has been used to obtain frequency-dependent Raman intensities,^{172–174} by finite difference of optical spectra, even for the second-order Raman spectrum.

The MBPT calculations can be combined with the temperature dependent electronic structure presented above, Sec. V E, to compute temperature-dependent optical spectra. The effect of the electron–phonon interaction is included in the diagonal part of the BSE Hamiltonian according to $Hvck,v\u2032c\u2032k\u2032(T)=Hvck,v\u2032c\u2032k\u2032FA+\Delta \epsilon ck(T)\u2212\Delta \epsilon vk(T)\delta vv\u2032\delta cc\u2032\delta kk\u2032,$ where *H*^{FA} is the “frozen-atom” expression for the BSE Hamiltonian introduced in Ref. 5. Figure 6 illustrates the agreement observed for silicon, when compared with experimental data, as well as with Ref. 175. Note that the Bethe–Salpeter Hamiltonian is not Hermitian but can be tackled by an iterative bi-Lanczos algorithm.^{5,174}

## VI. COMMUNITY PROJECTS: PSEUDOPOTENTIALS, HIGH-THROUGHPUT, AND THE libpaw LIBRARY

### A. Pseudopotentials and PAW datasets

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 set.^{7} On the one hand, this approach allows one to reduce significantly the computational cost and enjoy the advantages of the planewave basis set: orthogonality, systematic convergence, ease of implementation, and efficient fast-Fourier-transforms. On the other hand, as these pseudopotential operators are supposed to mimic the KS potential felt by the valence electrons in an all-electron DFT calculation, they imply a trade-off between accuracy and computational efficiency. Generating reliable and accurate tables of pseudopotentials therefore represents a highly nontrivial task, especially for end-users who are not familiar with the pseudopotential formalism and all its intricacies. For this reason, in the last few years, a significant effort has been made by the abinit community in order to provide users with a *recommended* set of pseudopotential tables that have been carefully crafted and validated against ground-state reference results obtained with all-electron codes.^{177,178} In what follows, we briefly describe the main features of the two official PAW and norm-conserving (NC) tables provided by the abinit group, and their design principles.

#### 1. jth table

To perform calculations in the frame of the PAW^{44} methodology, atomic data are needed, which correspond to the pseudopotentials in the norm-conserving approach. We deliver a Mendeleev table of PAW data—the jth table^{179}—on the abinit website. This table has been generated thanks to the atomic code atompaw^{180} for 86 elements from H to Rn. It is available within the scalar relativistic approximation both for LDA and PBE XC functionals. The files are provided in a standard XML format for use in any PAW electronic structure code, following the specifications given in Ref. 181. The XML format allows one to introduce new tags for special developments: for instance, pre-calculated atomic matrix elements are given, to be used for calculations with hybrid functionals. In the v2.0 version of the jth table, data enabling calculations within the LDA-1/2^{182} approximation are provided for some elements. Each file contains the suggested planewave cutoffs, according to a requested accuracy. The input file of the ATOMPAW code is also provided for each element to allow the user to modify the parameters used for the generation of the data, if needed.

#### 2. PseudoDojo

The pseudodojo project^{183} provides two predefined sets of NC pseudopotentials: *standard* and *stringent*. The *standard* set is designed for conventional DFT and DFPT applications, while the *stringent* version contains pseudopotentials with more electrons in valence, smaller pseudization radii, and improved scattering properties at high energies. The *stringent* set is recommended for performing DFT/DFPT calculations in which high precision is needed or for many-body applications such as *GW* that are quite sensitive to the inclusion of semi-core states and to the quality of the logarithmic derivative in the empty region. The majority of the pseudopotentials include non-linear core corrections^{184} with model core charges generated following Teter’s approach,^{185} which produces smooth core charges in real space with reasonably fast decay in reciprocal space. The inclusion of the non-linear core correction improves the transferability of the pseudopotential and turned out to be crucial to avoid unphysical oscillations in the local part of the potential in the region around the atom. These oscillations worsen the convergence rate, as well as the fulfillment of the acoustic sum rule in phonon calculations, in particular when GGA functionals are employed.^{183}

At present, the pseudodojo provides NC tables generated with three different XC functionals: LDA, PBE, and PBE-sol. For each XC flavor, one can opt for the scalar-relativistic or for the fully relativistic version, which includes the additional projectors required to treat spin–orbit coupling.^{186} Finally, a *specialized* set of pseudopotentials for the lanthanides is also available, with *f*-electrons frozen in the core. These pseudopotentials are supposed to be used for calculations in which the lanthanide is in the 3^{+} oxidation state. For each pseudopotential, three different hints (low, normal, and high) for the planewave cutoff energy are suggested on the basis of the convergence studies performed during the validation tests.^{183} These cutoff hints are employed in high-throughput studies, to implement machine-friendly workflows, but can also be used as a starting point in more conventional convergence studies. The pseudopotentials can be downloaded from the official website or alternatively from a GitHub repository. The data are available in three different formats: the original *psp*8 file format implemented by Hamann, the UPF format, and the *psml* format.^{187} The graphical interface of the official website allows users to select an element of the periodic table, the XC functional as well as the accuracy level, and the relativistic version. Jupyter notebooks with pre-generated figures showing the convergence of physical properties as a function of the cutoff energy and the results produced by the ONCVPSP pseudopotential generator are also provided. The same interface can be used to download the PAW atomic data of the jth table.

### B. High-throughput: Phonons, second-harmonic generation, and *GW*

Using a first-principles package for high-throughput applications places more stringent demands on the implementation details than the traditional approach based on single-step operations. abinit implements two important features to accommodate these needs: portable machine-readable output files in the binary netcdf format^{188,189} to communicate results to external software, as well as automatic algorithms to select optimal parallelization settings at runtime. Besides these features directly implemented in the code itself, efficient high-throughput calculations require a high-level library to programmatically interact with the Fortran code; the abinit group provides this framework in the abipy project.^{190} These features have enabled the development of the following high-throughput projects.

#### 1. DFPT phonons

The vibrational properties of a material represent a fundamental ingredient for understanding a variety of physical phenomena, and abinit provides accurate and efficient algorithms based on the DFPT formalism (Sec. V). In order to fully exploit the features of the code and handle a large number of calculations, a new python project has been developed, abiflows. The package contains all the functionalities required to execute high-throughput workflows with the fireworks framework^{191} and to store the results in a MongoDB database in a standardized format for later queries and analysis. abiflows heavily relies on the API implemented in abipy not only for the aforementioned properties but also for the automatic generation of abinit inputs for different kinds of calculations and the post-processing of results. In order to provide sensible input configurations for phonon calculations and avoid expensive convergence studies for each material, a global study has been performed and a set of heuristic optimal parameters has been determined.^{192} These parameters have been used to compute in an automated way the dynamical matrix, the phonon dispersion and the density of states for more than 2000 materials.^{193} Taking advantage of the post-processing tools implemented in anaddb, a set of derived quantities has been obtained for these materials. These include the dielectric permittivity tensors, Born effective charges, temperature dependent thermodynamic properties (entropy, heat capacity, Helmholtz free energy, and vibrational internal energy), frequency-dependent dielectric tensor, speed of sound along high symmetry directions, and the Debye–Waller tensor.^{153} All these quantities are freely available on the materials project website.^{194}

#### 2. DFPT SHG

The workflow for the calculation of vibrational properties discussed in the previous paragraph has been extended to compute the optical dielectric tensor *ε*^{∞} and the second-order susceptibility *χ*^{(2)} (see Sec. V). This gives access to optical properties beyond the linear regime with minimum effort. For instance, the second-order susceptibility leads to the nonlinear Second Harmonic Generation (SHG) coefficient *d*_{eff}.^{195} SHG processes play an important role in modern optics, especially in laser-related science and technology.^{196} Only non-centrosymmetric materials are SHG-active, and the calculation of the *d*_{eff} is possible in abinit through the 2*n* + 1 theorem.^{9,66}

By selecting the non-centrosymmetric materials from databases^{197,198} in which the optical dielectric tensors were already available, we built up a new set of candidates including more than 800 compounds aiming at calculating their second-order susceptibility. At this stage, we have computed the nonlinear properties for more than 400 materials employing the above-mentioned abiflows workflow. The results of this study will be presented in more detail in a future paper, and many more materials are being added to the list.

#### 3. GW

Computing the excited-state properties of materials (such as bandgaps) implies going beyond standard DFT and its fundamental shortcomings. MBPT provides a rigorous and well-established formalism for accurate band-structure calculations.^{11} The main ingredients of MBPT are the Green’s function *G*, the screened Coulomb interaction *W*, and the electronic self-energy Σ that is usually approximated with the *GW* method.^{199} Automating *GW* calculations at the high-throughput level represents, however, a significant challenge. Compared to the DFT electron density in which only occupied states are needed, the Green’s function requires the knowledge of all unoccupied states of the systems. In practical implementations, the sum over an infinite number of bands is avoided by introducing a cutoff (**nband**) beyond which the contribution of the high-energy states is considered negligible and the Fourier expansion of the two-point function *W* is restricted to **G**-vectors lying inside a sphere of kinetic energy **ecuteps**. The required computational resources increase steeply when these two cutoffs are increased. An additional problem is that these two convergence parameters are linked, i.e., the convergence with respect to one cutoff is determined by the value of the other parameter with a rate that is system-dependent. The level of coupling, moreover, varies strongly from one material to another. Determining reliable input parameters is, hence, of crucial importance to achieve precise and efficient MBPT calculations.

To enable high-throughput *GW* calculations, an automated convergence framework has been developed using the abipy package.^{200,201} A typical workflow is illustrated in Fig. 7. At a low density **k**-grid, the two interdependent cutoffs are treated in the same loop ensuring convergence of both. The convergence in **k**-space is treated separately since it was found to be unconnected. In the last step, the convergence rates of the screening cutoff and the unoccupied states cutoff at high **k**-space sampling are compared to those at low sampling density to finally ensure converged results. This method was used to compute the *GW* correction for about 90 crystalline compounds.^{200}

### C. The libpaw library

The abinit implementation of PAW^{44}—although not the only one—is one of the most complete ones present in DFT codes. Ground-state properties (DFT), excited states (MBPT), and response functions (DFPT) can be computed within PAW. However, the PAW method should theoretically be implemented in the same way in all DFT codes, regardless of the basis used to represent the electronic wavefunctions.

From our original PAW implementation,^{42} we have developed a portable PAW library called libpaw^{202} in order to facilitate the porting of PAW in other codes. This library provides a formalism for computational physics—the pseudopotential/PAW approach—at a relatively high implementation level (not only low-level methods, such as I/O, error-handling, or system resolution). It is packaged in such a way that it gives access to implementations of core PAW procedures with basis-independent data interfaces and thus can potentially be combined with any other basis used in different programs.

At present, the libpaw library—initially used in abinit—has been integrated in at least two other codes: the wavelet-based BIGDFT software^{202,203} as well as a recently developed Gaussian-type atomic orbital based DFT code.^{204} The use of libpaw in three codes, using three different types of basis functions (planewaves, wavelets, and Gaussians) highlights its portability. Other projects based on the libpaw are in progress.

The library is currently built “on the fly” from abinit source files. A standalone and abinit independent package is generated and can be directly inserted in a “host” code. At present, this is the only way to get the libpaw library; in the near future, the libpaw package will be directly downloadable as a separate package.

As soon as a host code uses the libpaw library, it can access at low cost a full PAW implementation and take the advantage of regular library updates, new features, and debugging. Some features are automatically available because they are only implemented in “on-site” PAW contributions such as, for example, Hubbard Hamiltonian calculations (DFT+*U*). The library provides a module to access PAW atomic datasets in the XML format and in the abinit proprietary format (legacy). The XML format^{181} then gives access to several PAW atomic dataset tables: jth,^{179} gbrv,^{177} atompaw,^{205} and gpaw.^{206}

The library fits very easily into the host code. A simple libpaw.h header file has to be adapted to integrate the LIBPAW. In the libpaw.h file, one has to specify which routines are used for the I/O and the error handling, as well as the dependencies shared with the LIBPAW (e.g., the libxc electronic exchange and correlation library^{207} or netcdf interfaces). Most of the low-level dependencies can be shared with the host code. The library is made of a complete collection of Fortran 2003 modules. Therefore, it can most easily be incorporated in codes written in Fortran. Each Fortran module is related to a specific “object” of the PAW formalism and contains associated methods.

Among the provided objects, the following are worth mentioning (a more complete description is available in Ref. 202):

m_pawrad: contains all functions related to the PAW radial meshes and associated derivation/integration routines for different kinds of meshes (linear or logarithmic).

m_pawpsp: used to read PAW atomic datasets (XML or abinit legacy formats),

m_pawtab: used to define non-self-consistent tabulated PAW data either read from the atomic dataset or directly deduced from it,

m_pawxc: computes the exchange and correlation potential/energy in the PAW augmentation regions using developments over spherical harmonics,

m_pawcprj: calculates, stores, and manipulates projections of pseudo-wavefunctions $\Psi \u0303nk$ on PAW non-local projectors $p\u0303i$: $cnki=\u27e8p\u0303i|\Psi \u0303nk\u27e9$. Note that the routine dedicated to the computation of these projections has to be provided by the host code,

m_pawrhoij: computes and manipulates the PAW occupancy matrix $\rho ij=\u2211nkcnki*cnk\u2009j$, especially ensuring the use of symmetries, and

m_pawdij: computes and manipulates the PAW non-local potential intensities

*D*_{ij}(see Ref. 42 for a complete expression). In this routine, all the physical ingredients are included (e.g., electronic correlations, spin–orbit, and hybrid functionals).

All available methods are fully compatible with the *Message Passing Interface* (MPI) distributed parallelism, including parallelization over atomic sites, as well as over the real space grid sampling the augmentation regions. A MPI communicator has to be provided.

The libpaw library is a perfect example of “code sharing” in the perspective of outsourcing and sharing components common to DFT pseudopotential-based codes. The advantages of code sharing are numerous:

It enforces clean and structured object-oriented programming (not obvious in the case of Fortran). Strong links to a specific code must be removed, and strong optimization can be implemented once and used widely.

It allows one to implement a cross-validation process, thanks to the use of the code in different contexts (here, different electronic wavefunction bases).

It allows one to substantially expand the user base and thus to reinforce user feedback.

It avoids developers re-coding what already exists and concentrates human work in implementations of novel physical quantities.

In the case of abinit, this code sharing effort is only in its early stages; other features could be easily packaged and shared: the handling of crystal symmetries, sampling of the Brillouin zone, self-consistent cycle preconditioning/mixing algorithms, etc.

## VII. CONCLUSION

With the increasing use of DFT for materials science applications, the need to have user-friendly packages that calculate and predict structural, electronic, vibrational, and elastic properties is essential. In this paper, we have presented the case of abinit, an electronic structure package for materials and nanosystem simulations. Although abinit originally focused on the “simple” solution of the ground state Kohn–Sham equations of DFT, it has broadened to include many different theories with far reaching applications, including DFT, DFPT, TD-DFT, MBPT, and DMFT. Similarly, two different treatments of core electrons are provided: norm-conserving pseudopotentials and projector-augmented waves. Each of these methodologies implements a different set of properties, which are summarized in Secs. II–VI. A list of the capabilities offered by abinit has been presented, with particular attention given to the graphical interfaces created by the abipy workflow manager and the post-processing tool agate.

In this paper, we have also concentrated on the properties that make abinit rather unique. We started with DMFT, the application of finite electric fields, probes of nuclei properties (such as the Mössbauer and nuclear spectroscopy), and positron annihilation. We then moved to response functions that are the most developed features of abinit: sections on Raman spectroscopy, electro-optic effects and strain perturbation, van der Waals interactions, flexoelectricity, the electron–phonon interaction, anharmonic interatomic force constants, non-collinear magnetism, spin magnetic field, and temperature dependence of the optical response with the Bethe–Salpeter equation.

After examining individual properties, we presented some of our efforts in community projects and library generation, which can be used by abinit and many other electronic structure codes. For example, we mention the distribution of reliable and well tested PAW datasets, libpaw, the pseudodojo project, and high-throughput developments to calculate advanced properties in the MBPT and DFPT formalisms.

In conclusion, we have presented here the evolution of abinit from its history to the most recent developments. We have also included links to the code documentation and tutorials. The basic theory and frameworks were discussed, which allow users to calculate a broad set of properties. Beyond the features with unique implementations in abinit, we encourage interested users to visit the abinit web page for more general documentation, tutorials, and information.

## ACKNOWLEDGMENTS

This work (A.H.R.) was supported by the DMREF-NSF Grant No. 1434897, National Science Foundation OAC-1740111, and U.S. Department of Energy DE-SC0016176 and DE-SC0019491 projects.

N.A.P. and M.J.V. gratefully acknowledge funding from the Belgian Fonds National de la Recherche Scientifique (FNRS) under Grant No. PDR T.1077.15-1/7. M.J.V. also acknowledges a sabbatical “OUT” grant at ICN2 Barcelona as well as ULiège and the Communauté Française de Belgique (Grant No. ARC AIMED G.A. 15/19-09).

X.G. and M.J.V. acknowledge funding from the FNRS under Grant No. T.0103.19-ALPS.

X.G. and G.-M. R. acknowledge support from the Communauté française de Belgique through the SURFASCOPE Project (No. ARC 19/24-057).

X.G. acknowledges the hospitality of the CEA DAM-DIF during the year 2017.

G.H. acknowledges support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DE-AC02-05-CH11231 (Materials Project Program No. KC23MP).

The Belgian authors acknowledge computational resources from supercomputing facilities of the University of Liège, the Consortium des Equipements de Calcul Intensif (Grant No. FRS-FNRS G.A. 2.5020.11), and Zenobe/CENAERO funded by the Walloon Region under Grant No. G.A. 1117545.

M.C. and O.G. acknowledge support from the Fonds de Recherche du Québec Nature et Technologie (FRQ-NT), Canada, and the Natural Sciences and Engineering Research Council of Canada (NSERC) under Grant No. RGPIN-2016-06666.

The implementation of the libpaw library (M.T., T.R., and D.C.) was supported by the ANR NEWCASTLE project (Grant No. ANR-2010-COSI-005-01) of the French National Research Agency.

M.R. and M.S. acknowledge funding from Ministerio de Economia, Industria y Competitividad (MINECO-Spain) (Grants Nos. MAT2016-77100-C2-2-P and SEV-2015-0496) and Generalitat de Catalunya (Grant No. 2017 SGR1506). This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation program (Grant Agreement No. 724529).

P.G. acknowledges support from FNRS Belgium through PDR (Grant No. HiT4FiT), ULiège and the Communauté française de Belgique through the ARC project AIMED, the EU and FNRS through M.ERA.NET project SIOX, and the European Funds for Regional Developments (FEDER) and the Walloon Region in the framework of the operational program “Wallonie-2020.EU” through the project Multifunctional thin films/LoCoTED.

The Flatiron Institute is a division of the Simons Foundation.

A large part of the data presented in this paper is available directly from the Abinit Web page www.abinit.org. Any other data not appearing in this web page can be provided by the corresponding author upon reasonable request.

## REFERENCES

_{2}O

_{3}and Pu

_{2}O

_{3}with the Hubbard I solver and comparison to DFT + U

_{2}RuO

_{4}from first-principles electronic structure

_{3}

^{17}O NMR parameters in SiO

_{2}: Assignment of the zeolite ferrierite spectrum

_{2}

_{2}

^{0}ferromagnetism in gallium nitride: A two-component density functional theory simulation

_{3}double perovskites

_{3}perovskite and post-perovskite at high pressure

_{3}and PbTiO

_{3}prototypical ferroelectrics from density functional theory

_{1}ab ferroelectric phase of aurivillius Bi

_{2}WO

_{6}

_{3}thin films

_{3}from first principles

_{2}PO

_{4}including dispersion forces: Stability and related vibrational, dielectric, and elastic properties

_{3}O

_{6}

_{2}, W, and the (1 × 1) H-covered W(110) surface

_{3}from first principles

_{1−x}K

_{x}BiO

_{3}

_{2}-type SiO

_{2}at high pressure

_{2}α-quartz and stishovite

Note that the meaning of ixcrot was unfortunately inverted in Ref. 166 with respect to the actual implementation.