**MAE715-Atomistic Modeling of Materials**

MAE715-Atomistic Modeling of Materials

(Parte **3** de 7)

4.1. PWscf

PWscf implements an iterative approach to reach self-consistency, using at each step iterative diagonalization techniques, in the framework of the plane-wave pseudopotential method. An early version of PWscf is described in Ref [78].

Both separable NC-PPs and US-PPs are implemented; recently, also the projector augmented-wave method [37] has been added, largely following the lines of Ref [79] for its implementation. In the case of US-PPs, the electronic wave functions can be made smoother at the price of having to augment their square modulus with additional contributions to recover the actual physical charge densities. For this reason, the charge density is more structured than the square of the wavefunctions, and requires a larger energy cutoff for its

QUANTUM ESPRESSO 12 plane wave expansion (typically, 6 to 12 times larger; for a NC-P, a factor of 4 would be mathematically sufficient). Hence, different real-space Fourier grids are introduced - a "soft" one that represents the square of electronic wave functions, and a "hard" one that represents the charge density [80, 81]. The augmentation terms can be added either in reciprocal space (using an exact but expensive algorithm) or directly in real space (using an approximate but faster algorithm that exploits the local character of the augmentation charges).

PWscf can use the well established LDA and GGA exchange-correlation functionals, including spin-polarization within the scheme proposed in Ref [82] and can treat non-collinear magnetism[48, 49] as e.g. induced by relativistic effects (spin-orbit interactions) [83, 84] or by complex magnetic interactions (e.g. in the presence of frustration). DFT + Hubbard U calculations [85] are implemented for a simplified (“no-J”) rotationally invariant form [86] of the Hubbard term. Other advanced functionals include TPSS meta-GGA [39], functionals with finite-size corrections [87], and the PBE0 [40] and B3LYP [41, 42] hybrids.

Self-consistency is achieved via the modified Broyden method of Ref [8], with some further refinements that are detailed in Appendix A.1. The sampling of the Brillouin Zone (BZ) can be performed using either special [89, 90] k-points provided in input or those automatically calculated starting from a uniform grid. Crystal symmetries are automatically detected and exploited to reduce computational costs, by restricting the sampling of the BZ to the irreducible wedge alone (See Appendix A.4). When only the Γ point (k = 0) is used, advantage is taken of the real character of the KS orbitals, allowing one to store just half of the Fourier components. BZ integrations in metallic systems can be performed using a variety of smearing/broadening techniques, such as Fermi-Dirac, Gaussian, Methfessel-Paxton [91], and Marzari-Vanderbilt cold smearing [92]. The tetrahedron method [93] is also implemented. Finite-temperature effects on the electronic properties can be easily accounted for by using the Fermi-Dirac smearing as a practical way of implementing the Mermin finite-temperature density-functional approach [94].

Structural optimizations are performed using the Broyden-Fletcher-Goldfarb-Shanno

(BFGS) algorithm [95, 96, 97] or damped dynamics; these can involve both the internal, microscopic degrees of freedom (i.e. the atomic coordinates) and/or the macroscopic ones (shape and size of the unit cell). The calculation of minimum-energy paths, activation energies, and transition states uses the Nudged Elastic Band (NEB) method [57]. Potential energy surfaces as a function of suitably chosen collective variables can be studied using Laio-Parrinello metadynamics [98].

Microcanonical (NVE) MD is performed on the BO surface, i.e. achieving electron selfconsistency at each time step, using the Verlet algorithm[9]. Canonical (NVT) dynamics can be performed using velocity rescaling, or Anderson’s or Berendsen’s thermostats [100]. Constant-pressure (NPT) MD is performed by adding additional degrees of freedom for the cell size and volume, using either the Parrinello-Rahman Lagrangian [101] or the so-called invariant Lagrangian of Wentzcovitch [53].

The effects of finite macroscopic electric fields on the electronic structure of the ground state can be accounted for either through the method of Ref [102, 103] based on the Berry phase, or (for slab geometries only) through a sawtooth external potential [104, 105]. A

QUANTUM ESPRESSO 13 quantum fragment can be embedded in a complex electrostatic environment that includes a model solvent [106] and a counterion distribution [107], as is typical of electrochemical systems.

4.2. CP

The CP code is the specialized module performing Car-Parrinello ab initio MD. CP can use both NC PPs [108] and US PPs [109, 80]. In the latter case, the electron density is augmented through a Fourier interpolation scheme in real space (“box grid”) [80, 81] that is particularly efficient for large scale calculations. CP implements the same functionals as PWscf, with the exception of hybrid functionals; a simplified one-electron self-interaction correction (SIC)[110] is also available. The Car-Parrinello Lagrangian can be augmented with Hubbard U corrections [1], or Hubbard-based penalty functionals to impose arbitrary oxidation states [112].

Since the main applications of CP are for large systems without translational symmetry (e.g. liquids, amorphous materials), Brillouin zone sampling is restricted to the Γ point of the supercell, allowing for real instead of complex wavefunctions. Metallic systems can be treated in the framework of “ensemble DFT” [113].

In the Car-Parrinello algorithm, microcanonical (NVE) MD is performed on both electronic and nuclear degrees of freedom, treated on the same footing, using the Verlet algorithm. The electronic equations of motion are accelerated through a preconditioning scheme [114]. Constant-pressure (NPT) MD is performed using the Parrinello-Rahman Lagrangian [101] and additional degrees of freedom for the cell. Nosé-Hoover thermostats [115] and Nosé-Hoover chains [116] allow to perform simulations in the different canonical ensembles.

CP can also be used to directly minimize the energy functional to self-consistency while keeping the nuclei fixed, or to perform structural minimizations of nuclear positions, using the “global minimization” approaches of Refs. [117, 118], and damped dynamics or conjugate-gradients on the electronic or ionic degrees of freedom. It can also perform NEB and metadynamics calculations.

Finite homogeneous electric fields can be accounted for using the Berry phase method, adapted to systems with the Γ point only [102]. This advanced feature can be used in combination with MD to obtain the infrared spectra of liquids [102, 119], the low- and highfrequency dielectric constants [102, 120] and the coupling factors required for the calculation of vibrational properties, including infrared, Raman [121, 122, 123], and hyper-Raman [124] spectra.

4.3. PHonon

The PHonon package implements density-functional perturbation theory (DFPT) [54, 5, 56] for the calculation of second- and third-order derivatives of the energy with respect to atomic displacements and to electric fields. The global minimization approach [125, 126] is used for the special case of normal modes in finite (molecular) systems, where no BZ sampling is

QUANTUM ESPRESSO 14 required (Gamma code). In the general case a self-consistent procedure [5] is used, with the distinctive advantage that the response to a perturbation of any arbitrary wavelength can be calculated with a computational cost that is proportional to that of the unperturbed system. Thus, the response at any wavevector, including very small (long-wavelength) ones, can be inexpensively calculated. This latter approach, and the technicalities involved in the calculation of effective charges and interatomic force constants, are described in detail in Refs. [5, 127] and implemented in the PH code.

Symmetry is fully exploited in order to reduce the amount of computation. Lattice distortions transforming according to irreducible representations of small dimensions are generated first. The charge-density response to these lattice distortions is then sampled at a number of discrete k-points in the BZ, which is reduced according to the symmetry of the small group of the phonon wavevector q. The grid of the q points needed for the calculation of interatomic force constants reduces to one wavevector per star: the dynamical matrices at the other q vectors in the star are generated using the symmetry operations of the crystal. This approach allows us to speed up the calculation without the need to store too much data for symmetrization.

The calculation of second-order derivatives of the energy works also for US P [128, 129] and for all GGA flavors [130, 131] used in PWscf and in CP. The extension of PHonon to PAW [132], to noncollinear magnetism and to fully relativistic US PPs which include spinorbit coupling [133] will be available by the time this paper will be printed.

Advanced features of the PHonon package include the calculation of third order derivatives of the energy, such as electron-phonon or phonon-phonon interaction coefficients. Electron-phonon interactions are straightforwardly calculated from the response of the selfconsistent potential to a lattice distortion. This involves a numerically-sensitive “doubledelta” integration at the Fermi energy, that is performed using interpolations on a dense k-point grid. Interpolation techniques based on Wannier functions [134] will speed up considerably these calculations. The calculation of the anharmonic force constants from third-order derivatives of the electronic ground-state energy is described in Ref. [135] and is performed by a separate code called d3. Static Raman coefficients are calculated using the second-order response approach of Refs. [136, 137]. Both third-order derivatives and Raman-coefficients calculations are currently implemented only for NC PPs.

4.4. atomic

The atomic code performs three different tasks: i) solution of the self-consistent all-electron radial KS equations (with a Coulomb nuclear potential and spherically symmetric charge density); i) generation of NC PPs, of US PPs, or of PAW data-sets; ii) test of the above PPs and data-sets. These three tasks can be either separately executed or performed in a single run. Three different all-electron equations are available: i) the non relativistic radial KS equations, i) the scalar relativistic approximation to the radial Dirac equations [138], ii) the radial Diraclike equations derived within relativistic density functional theory [139, 140]. For i) and i) atomic magnetism is dealt with within the local spin density approximation, i.e. assuming

QUANTUM ESPRESSO 15 an axis of magnetization. The atomic code uses the same exchange and correlation energy routines of PWscf and can deal with the same functionals.

The code is able to generate NC PPs directly in separable form (also with multiple projectors per angular momentum channel) via the Troullier-Martins [141] or the Rappe- Rabe-Kaxiras-Joannopoulos [142] pseudization. US PPs can be generated by a two-step pseudization process, starting from a NC PPs, as described in Ref. [143], or using the solutions of the all-electron equation and pseudizing the augmentation functions [80]. The latter method is used also for the PAW data-set generation. The generation of fully relativistic NC and US PPs including spin-orbit coupling effects is also available. Converters are available to translate pseudopotentials encoded in different formats (e.g. according to the Fritz-Haber [75] or Vanderbilt [76] conventions) into the UPF format adopted by QUANTUM ESPRESSO.

Transferability tests can be made simultaneously for several atomic configurations with or without spin-polarization, by solving the non relativistic radial KS equations generalized for separable nonlocal PPs and for the presence of an overlap matrix.

4.5. PWcond

The PWcond code implements the scattering approach proposed by Choi and Ihm [60] for the study of coherent electron transport in atomic-sized nanocontacts within the Landauer- Büttiker theory. Within this scheme the linear response ballistic conductance is proportional to the quantum-mechanical electron transmission at the Fermi energy for an open quantum system consisting of a scattering region (e.g., an atomic chain or a molecule with some portions of left and right leads) connected ideally from both sides to semi-infinite metallic leads. The transmission is evaluated by solving the KS equations with the boundary conditions that an electron coming from the left lead and propagating rightwards gets partially reflected and partially transmitted by the scattering region. The total transmission is obtained by summing all transmission probabilities for all the propagating channels in the left lead. As a byproduct of the method, the PWcond code provides the complex band structures of the leads, that is the Bloch states with complex kz in the direction of transport, describing wave functions exponentially growing or decaying in the z direction. The original method formulated with

NC PPs has been generalized to US PPs both in the scalar relativistic [144] and in the fully relativistic forms [145].

4.6. GIPAW

The GIPAW code allows for the calculation of physical parameters measured by i) nuclear magnetic resonance (NMR) in insulators (the electric field gradient (EFG) tensors and the chemical shift tensors), and by i) electronic paramagnetic resonance (EPR) spectroscopy for paramagnetic defects in solids or in radicals (the hyperfine tensors and the g-tensor). The code also computes the magnetic susceptibility of nonmagnetic insulators. GIPAW is based on the PW-P method, and uses many subroutines of PWscf and of PHonon. The code is currently restricted to NC PPs. All the NMR and EPR parameters depend on the detailed shape of the electronic wave-functions near the nuclei and thus require the reconstruction of

QUANTUM ESPRESSO 16 the all-electron wave-functions from the P wave-functions. For the properties defined at zero external magnetic field, namely the EFG and the hyperfine tensors, such reconstruction is performed as a post-processing step of a self-consistent calculation using the PAW reconstruction, as described for the EFG in Ref. [146] and for the hyperfine tensor in Ref. [147]. The g-tensor, the NMR chemical shifts and the magnetic susceptibility are obtained from the orbital linear response to an external uniform magnetic field. In the presence of a magnetic field the PAW method is no more gauge- and translationally invariant. Gauge and translational invariances are restored by using the gauge including projector augmented wave (GIPAW) method [63, 64] both i) to describe in the P Hamiltonian the coupling of orbital degrees of freedom with the external magnetic field, and i) to reconstruct the all-electron wave-functions, in presence of the external magnetic field. In addition, the description of a uniform magnetic field within periodic boundary conditions is achieved by considering the long wave-length limit of a sinusoidally modulated field in real space [148, 149]. The NMR chemical shifts are computed following the method described in Ref. [63], the g-tensor following Ref. [150] and the magnetic susceptibility following Refs. [148, 63]. Recently, a “converse” approach to calculate chemical shifts has also been introduced [151], based on recent developments on the Berry-phase theory of orbital magnetization; since it does not require a linear-response calculation, it can be straightforwardly applied to arbitrarily complex exchange-correlation functionals, and to very large systems, albeit at a computational cost that is proportional to the number of chemical shifts that need to be calculated.

(Parte **3** de 7)