QuantumATK: an integrated platform of electronic and atomic-scale modelling tools

QuantumATK is an integrated set of atomic-scale modelling tools developed since 2003 by professional software engineers in collaboration with academic researchers. While different aspects and individual modules of the platform have been previously presented, the purpose of this paper is to give a general overview of the platform. The QuantumATK simulation engines enable electronic-structure calculations using density functional theory or tight-binding model Hamiltonians, and also offers bonded or reactive empirical force fields in many different parametrizations. Density functional theory is implemented using either a plane-wave basis or expansion of electronic states in a linear combination of atomic orbitals. The platform includes a long list of advanced modules, including Green’s-function methods for electron transport simulations and surface calculations, first-principles electron-phonon and electron-photon couplings, simulation of atomic-scale heat transport, ion dynamics, spintronics, optical properties of materials, static polarization, and more. Seamless integration of the different simulation engines into a common platform allows for easy combination of different simulation methods into complex workflows. Besides giving a general overview and presenting a number of implementation details not previously published, we also present four different application examples. These are calculations of the phonon-limited mobility of Cu, Ag and Au, electron transport in a gated 2D device, multi-model simulation of lithium ion drift through a battery cathode in an external electric field, and electronic-structure calculations of the composition-dependent band gap of SiGe alloys.

Atomic-scale modelling is increasingly important for industrial and academic research and development in a wide range of technology areas, including semiconductors, 1,2 batteries, 3 catalysis, 4 renewable energy, 5 advanced materials, 6 next-generation pharmaceuticals, 7 and many others. Surveys indicate that the return on investment of atomic-scale modelling is typically around 5:1. 8 With development of increasingly advanced simulation algorithms and more powerful computers, we expect that the economic benefits of atomic-scale modelling will only increase.
The current main application of atomic-scale modelling is in early-stage research into new materials and technology designs, see Refs. 9 and 10 for examples. The early research stage often has a very large design space, and experimental trial and error is a linear process that will explore only a small part of this space. Atomic-scale simulations make it possible to guide experimental investigations towards the most promising part of the technology design space. Such insights are typically achieved by simulating the underlying atomic-scale processes behind failed or successful experiments, to understand the physical or (bio-)chemical origins. Such insight can often rule out or focus research to certain designs or material systems. 8 Recently, materials screening has also shown great promise. In this approach, atomic-scale calculations are used to obtain important properties of a large pool of materials, and the most promising candidates are then selected for experimental verification and/or further theoretical refinement. 6,11,12 The scientific field of atomic-scale modelling covers everything from near-exact quantum chemical calculations to approximate simulations using empirical force fields. Quantum chemical methods (based on wave-function theory) attempt to fully solve the many-body Schrödinger equation for all electrons in the system, and can provide remarkably accurate descriptions of molecules. 13 However, the computational cost is high: in practice, one is usually limited to calculations involving far below 100 atoms in total. Such methods are currently not generally useful for industrial research into advanced materials and next-generation electronic devices.
On the contrary, force-field (FF) methods are empirical but computationally efficient: all inter-atomic interactions are described by analytic functions with preadjusted parameters. It is thereby possible in practice to simulate systems with millions of atoms. Unfortunately, this often also hampers the applicability of a force field for system types not included when fitting the FF parameters.
As an attractive intermediate methodology, density functional theory (DFT) [14][15][16][17] provides an approximate but computationally tractable solution to the electronic many-body problem. This allows for good predictive power with respect to experiments with minimal use of empirical parameters at a reduced computational cost. Standard DFT simulations may routinely be applied to systems containing more than one thousand atoms, and DFT is today the preferred framework for industrial applications of ab initio electronic-structure theory.
Semi-empirical (SE) electronic-structure methods based on tight-binding (TB) model Hamiltonians are more approximate, but have a long tradition in semiconductor research. 18 Whereas DFT ultimately aims to approximate the true many-body electronic Hamiltonian in an efficient but parameter-free fashion, a TB model relies on parameters that are adjusted to very accurately describe the properties of a number of reference systems. This leads to highly specialized electronic-structure models that typically reduce the computational expense by an order of magnitude compared to DFT methods. Such SE methods may be convenient for large-scale electronicstructure calculations, for example in simulations of electron transport in semiconductor devices.
The QuantumATK platform offers simulation engines covering the entire range of atomic-scale simulation methods relevant to the semiconductor industry and materials science in general. This includes force fields, SE methods, and several flavors of DFT. These are summarized in Table I, including examples of other platforms that offer similar methodology.
To give a bird's-eye view of the computational cost of the different atomic-scale simulation methods mentioned above, we compare in Fig. 1 the computational speed of the methods when simulating increasingly larger structures of amorphous Al 2 O 3 . The measure of speed is here the number of molecular dynamics steps that are feasible within 24 hours when run in parallel on 16 computing cores. Although the parallel computing techniques used may differ between some of the methods, we find that Fig. 1 gives a good overview of the scaling between the different methods.
It is important to realize that the simulation methods showing the total number of molecular dynamics steps performed in 24 hours (# MD steps) against system size (# atoms) for amorphous Al2O3 with constant density. Each step includes evaluation of the total energy and atomic forces. The simulations were run on a 16-core central processing unit (CPU) of the type Intel R Xeon R E5-2670. The FF simulations (Section VI) were performed using threading only, whereas full MPI parallelization was used for the TB (Section V) and DFT (Section IV) simulations. For the latter, we have considered either semi-local exchange-correlation functionals using linear-combination-ofatomic-orbitals (DFT-LCAO) and plane-wave (DFT-PW) basis sets, or a hybrid exchange-correlation functional using a PW basis set (Hybrid DFT). Further details of the calculations are given in Appendix A.
listed in Table I should ideally complement each other: for successful use of atomic-scale modelling, it is essential to have easy access to all the methods, in order to use them in combination. The vast majority of atomicscale simulation tools are developed by academic groups, and most of them focus on a single method. Using the tool typically requires a large effort for compilation, installation, learning the input/output syntax, etc. The tool is often not fully compatible with any other tool, so learning an additional tool within a new modelling class requires yet another large effort. Even within one modelling class, for example DFT, a single simulation tool may not have all the required functionality for a given application, so several different tools within each modelling class may be needed to solve a given problem, and a significant effort must be invested to master each of them. As a commercially developed platform, Quantu-mATK aims to circumvent these issues.
Academic development of atomic-scale simulation platforms, often made available through open-source licenses, is essential for further technical progress of the field. However, the importance of commercial platforms in progressing the industrial uptake of the technology is often underestimated. Commercial software relies on payment from end users. This results in a strong focus on satisfying end-user requirements in terms of usability, functionality, efficiency, reliability, and support. The revenue enables the commercial software provider to establish a stable team of developers and thereby provide a software solution that will be maintained, extended, and supported for decades.
The ambition of the QuantumATK platform is to provide a state-of-the-art and easy-to-use integrated toolbox with all important atomic-scale modelling methodologies for a growing number of application areas. The methods are made available through a modern graphical user interface (GUI) and a Python scripting-based frontend for expert users. Our current focus is semiconductor devices, polymers, glasses, catalysis, batteries, and materials science in general. In this context, semiconductor devices is a broad area, ranging from siliconbased electronic logics and memory elements, 31,32 to solar cells composed of novel materials 33 and next-generation electronic devices based on spintronic phenomona. 34 One key strength of a unified framework for a large selection of simulation engines and modelling tools is within multiphysics and multiscale problems. Such problems often arise in physical modelling of semiconductor devices, and the QuantumATK platform is widely used for coupling technology computer-aided design (TCAD) tools with atomic-scale detail, for instance to provide first-principles simulations of defect migration paths and subsequently the temperature-dependent diffusion constant for continuum-level simulation of semiconductor processes. 2 Furthermore, QuantumATK provides a highly flexible and efficient framework for coupling advanced electrostatic setups with state-of-the-art transport simulations including electron-phonon coupling and light-matter interaction. This has enabled predictions of gate-induced phonon scattering in graphene gate stacks, 35 atomistic description of ferroelectricity driven by edge-absorbed polar molecules in gated graphene, 36 and new 2D material science such as prediction of the room-temperature photocurrent in emerging layered Janus materials with a large dipole across the plane. 37 The flexibility of the QuantumATK framework supports the imagination of researchers, and at the same time enables solutions to both real-world and cutting-edge semiconductor device and material science problems.
The purpose of this paper is to give a general overview of the QuantumATK platform with appropriate references to more thorough descriptions of several aspects of the platform. We also provide application examples that illustrate how the different simulation engines can complement each other. The paper is organized as follows: In Section II we give a general overview of the QuantumATK platform, while Section III introduces the types of system geometries handled by the platform. The next three sections (IV-VI) describe the DFT, SE, and FF simulation engines, respectively. We then introduce a number of simulation modules that work with the different engines. These modules include ion dynamics (VII A), phonon properties (VIII), polarization (IX), magnetic anisotropy energy (X), and quantum transport (XI). We next describe the parallel computing strategies of the different engines, and present parallel scaling plots in Section XII. We then in Section XIII describe the scripting and GUI simulation environment in the Quan-tumATK platform. This is followed by four application examples in Section XIV, and the paper is summarized in Section XV.

II. OVERVIEW
The core of QuantumATK is implemented in C++ modules with Python bindings, such that all C++ modules are accessible from ATK-Python, a customized version of Python built into the software. The combination of a C++ backend and a Python-based frontend offers both high computational performance and a powerful but user-friendly scripting platform for setting up, running, and analyzing atomic-scale simulations. All simulation engines listed in Table I are invoked using ATK-Python scripting. More details are given in Section XIII A. Quan-tumATK also relies on a number of open-source packages, including high-performance numerical solvers.
All computationally demanding simulation modules may be run in parallel on many processors at once, using message passing between processes and/or sharedmemory threading, and often in a multi-level approach. More details are given in Section XII.
The full QuantumATK package is installed on Windows or Linux/Unix operating systems using a binary installer obtained from the Synopsys SolvNet website, https://solvnet.synopsys.com. All required external software libraries are pre-compiled and shipped with the installer. Licensing is handled using the Synopsys Common Licensing (SCL) system.

III. ATOMISTIC CONFIGURATIONS
The real-space physical system to be simulated is defined as an ATK-Python configuration object, including lattice vectors, element types and positions, etc. Quan-tumATK currently offers four main types of such configurations: molecule, bulk, device, and surface. Examples of these are given in Fig. 2.
The simplest configuration is the molecule configuration shown in Fig. 2(a). It is used for isolated (nonperiodic) systems, and is defined by a list of elements and their positions in Cartesian coordinates.
The bulk configuration, shown in Fig. 2(b), defines an atomic-scale system that repeats itself in one or more directions, for example a fully periodic crystal (periodic in 3D), a 2D nanosheet (or a slab), or a 1D nanowire. The bulk system is defined by the Bravais lattice and the position of the atomic elements inside the primitive cell.
The two-probe device configuration is used for quantum transport simulations. As shown in Fig. 2(c), the device consists of a central region connected to two semiinfinite bulk electrodes. The central region, where scattering of electrons travelling from one electrode to the other may take place, can be periodic in zero (1D wire), one (2D sheet), or two (3D bulk) directions, but is bounded by the electrodes along the third dimension. The device configuration is used to simulate electron and/or phonon transport via the non-equilibrium Green's function (NEGF) method. 38 Finally, for physically correct simulations of a surface, QuantumATK provides the one-probe surface configuration. This is basically a device configuration with only one electrode, as illustrated in Fig. 2(d). By construction, the surface configuration realistically describes the electronic structure of a semi-infinite crystal beyond the approximate slab model. 19 The remainder of this paper is devoted to describing the computational methods available for calculating the properties of such configurations using QuantumATK.

IV. DFT SIMULATION ENGINES
Density functional theory is implemented in the Kohn-Sham (KS) formulation [14][15][16][17] within the framework of the linear combination of atomic orbitals (LCAO) and planewave (PW) basis set approaches, combined with the pseudopotential method. The electronic system is seen as a non-interacting electron gas of density n in the effective potential V eff [n], where V H is the Hartree potential describing the classical electrostatic interaction between the electrons, V xc is the exchange-correlation (XC) potential, which in practise needs to be approximated, and V ext is the sum of the electrostatic potential energy of the electrons in the , and a right electrode (transparent yellow). Both electrodes are semi-infinite in the left-right transport direction. The device is in this case periodic in the two directions perpendicular to the transport direction, but would be nonperiodic in one or both perpendicular directions in case of a nanosheet or nanotube device, respectively. (d) Surface configuration of a gold surface. The structure consists of a left electrode (transparant yellow) and a central region (orange box). We note that an electric field can be applied to the surface by choice of boundary condition on the right-hand face of the central region.
external potential of ions and other electrostatic field sources. The total external potential is in QuantumATK given by where V pseudo a includes the local (V loc a ) and nonlocal (V nl a ) contributions to the pseudopotential of the a-th atom. The term V gate is a potential that may originate from other external sources of electrostatic fields, for example metallic gates.
The KS Hamiltonian consists of the single-electron kinetic energy and the effective potential, and the single-electron energies ( α ) and wave functions (ψ α ) are solutions to eigenvalue problem The electronic ground state is found by iteratively minimizing the KS total-energy density functional, E[n], with respect to the electron density, where T is the kinetic energy. The forces (acting on the atoms) and stress tensor of the electronic system may then be computed as derivatives of the ground-state total energy with respect to the atomic coordinates and the strain tensor, respectively.

A. LCAO Representation
The DFT-LCAO method uses a LCAO numerical representation of the KS equations, closely resembling the SIESTA formalism. 20 This allows for a localized matrix representation of the KS Hamiltonian in (3), and therefore an efficient implementation of KS-DFT for molecules, bulk materials, interface structures, and nanoscaled devices.
In the DFT-LCAO method, the single-electron KS eigenfunctions, ψ α , are expanded in a set of finite-range atomic-like basis functions φ i , The KS equation can then be represented as a matrix equation for determining the expansion coefficients c αi , where the Hamiltonian matrix H KS ij = φ i |Ĥ KS |φ j and overlap matrix S ij = φ i |φ j are given by integrals with respect to the electron coordinates. Two-center integrals are computed using 1D radial integration schemes employing a Fourier transform technique, while multiplecenter integrals are computed on a real-space grid. 20 For molecules and bulk systems, diagonalization of the Hamiltonian matrix yields the density matrix D ij , where f is the Fermi-Dirac distribution of electrons over energy states, ε F the Fermi energy, T the electron temperature, and k B the Boltzmann constant. For device and surface configurations, the density matrix is calculated using the NEGF method, as described in Section XI. The electron density is computed from the density matrix, and is represented on a regular real-space grid, which is the same grid as used for the effective potential in (1).

B. PW Representation
A PW representation of the KS equations was recently implemented in QuantumATK. It is complimentary to the LCAO representation discussed above. The ATK-PlaneWave engine is intended mainly for simulating bulk configuratins with periodic boundary conditions. The KS eigenfunctions are expanded in terms of PW basis functions, where α denotes both the wave vector k and the band index n, and g are reciprocal lattice vectors. The upper threshold for the reciprocal lattice-vector lengths included in the PW expansion (g max ) is determined by a kinetic-energy (wave-function) cutoff energy E cut , The DFT-PW method has its distinct advantages and disadvantages compared to the DFT-LCAO approach. In particular, the PW expansion is computationally efficient for relatively small bulk systems, and the obtained physical quantities can be systematically converged with respect to the PW basis-set size by increasing E cut . However, the PW representation is computationally inefficient for low-dimensional systems with large vacuum regions. It is also incompatible with the DFT-NEGF methodology for electron transport calculations in nanoscaled devices, unlike the LCAO representation, which is ideally suited for dealing with open boundary conditions, and is also more efficient for large systems.
The ATK-PlaneWave engine was implemented on the same infrastructure as used by the ATK-LCAO engine, though a number of routines were modified to reach stateof-the-art PW efficiency. For example, we have adopted iterative algorithms for solving the KS equations, 39 and fast Fourier transform (FFT) techniques for applying the Hamiltonian operator and evaluating the electron density. 40,41 In Fig. 3 we compare the CPU times of DFT-PW vs. DFT-LCAO calculations for different LCAO basis sets. The figure shows the CPU time for the different methods as function of the system size. The PW approach is computationally efficient for smaller systems, while the LCAO approach can be more than an order of magnitude faster for systems with more than 100 atoms.

C. Pseudopotentials and LCAO Basis Sets
QuantumATK uses pseudopotentials (PPs) to avoid explicit DFT calculations of core electrons, and currently supports both scalar-relativistic and fully relativistic normconserving PPs. 42 Projector augmented-wave (PAW) potentials 43 is currently available for the ATK-PlaneWave simulation engine only.
The QuantumATK platform is shipped with built-in databases of well-tested PPs, covering all elements up to Z = 83 (Bi), excluding lanthanides. The current default PPs are those of the published SG15 44 and PseudoDojo 45 sets. These are two modern normconserving PP types with multiple projectors for each angular momentum, to ensure high accuracy. Both sets contain scalar-relativistic and fully relativistic PPs for each element. The fully relativistic PPs are generated by solving the Dirac equation for the atom, which naturally includes spin-orbit coupling, and then mapping the solution onto the scalarrelativistic formalism. 42,46 For each PP, we have generated an optimized LCAO basis set, consisting of orbitals φ nlm , where Y lm are spherical harmonics, and χ nl are radial functions with compact support, being exactly zero outside a confinement radius. The basis orbitals are obtained by solving the radial Schrödinger equation for the atom in a confinement potential. 20 For the shape of the confinement potential, we follow Ref. 47.
To construct high-accuracy LCAO basis sets for the SG15 and PseudoDojo PPs, we have adopted a large set of pseudo-atomic orbitals that are similar to the "tight tier 2" basis sets used in the FHI-aims package. 47 These basis sets typically have 5 orbitals per PP valence electron, and a range of 5Å for all orbitals, and include angular momenta up to l = 5. From this large set, we have constructed three different series of reduced DFT-LCAO basis sets implemented in QuantumATK: 1. Ultra: Generated by reducing the range of the original pseudo-atomic orbitals, requiring that the overlap of each contracted orbital with the corresponding original orbital can change by no more than 0.1%. Also denoted "LCAO-U".
2. High: Generated by reducing the number of basis orbitals in the Ultra basis set, requiring that the DFT-obtained total energy of suitably chosen test systems change by no more than 1 meV/atom. Also denoted "LCAO-H".
3. Medium: Generated by further reduction of the number of orbitals in the High basis set, requiring that the subsequent change of the DFT-obtained total energies do not exceed 4 meV/atom. Also denoted "LCAO-M". The number of pseudo-atomic orbitals in a Medium basis set is typically comparable to that of a double-zeta polarized (DZP) basis set.
To validate the PPs and basis-sets, we have used the ∆-test 48,50 to check the accuracy of the equation of state for elemental, rock-salt, and perovskite solids against allelectron reference calculations, as shown in Table II. For each bulk crystal, the equation of state was calculated at fixed internal ion coordinates, and the equilibrium lattice constant and bulk modulus were computed. In Table II, the ∆-value is defined as the root-mean-square (RMS) energy difference between the equations of state obtained with QuantumATK and the all-electron reference, averaged over all crystals in a purely elemental benchmark set.
Table II suggests a general trend that the PseudoDojo PPs are slightly more accurate than the SG15 ones. Since the PseudoDojo PPs are in general also softer, requiring a lower real-space density mesh cutoff energy, these are the default PPs in QuantumATK.
Table II also shows that the accuracy of the DFT-LCAO calculations done with High or Ultra basis sets is rather close to that of the PW calculations. The Medium basis sets give on average a larger deviation from the PW results. However, we also find that LCAO-M provides sufficient accuracy for many applications, and it is therefore the default ATK-LCAO basis set in QuantumATK. We note that in typical applications, using Medium instead of the High (Ultra) basis sets decreases the computational cost by a factor of 2-4 (10-20), as seen in Fig. 3.
More details on the construction and validation of the LCAO basis sets can be found in Ref. 19.

D. Exchange-Correlation Methods
The XC functional in (5) is the only formal approximation in KS-DFT, since the exact functional is unknown. [15][16][17] QuantumATK supports a large range of approximate XC functionals, including the local density approximation (LDA), generalized gradient approximations (GGAs), and meta-GGA functionals, all supplied through the Libxc library. 51 The ATK-PlaneWave engine also allows for calculations using the HSE06 screened hybrid functional. [52][53][54] The ATK-LCAO and ATK-PlaneWave engines both support van der Waals dispersion methods using the two-body and three-body dispersion corrections by Grimme. 55 Both DFT engines support different spin variants for each XC functional: spinunpolarized and spin-polarized (both collinear and noncollinear). Spin-polarized noncollinear calculations may include spin-orbit interaction through the use of fully relativistic PPs.

Semilocal functionals
During the past 20 years, the semilocal (GGA) XC approximations have been widely used, owing to a good balance between accuracy and efficiency for DFT calculations. QuantumATK implements many of the popular GGAs, including the general-purpose PBE, 56 the PBEsol (designed for solids), 57 and the revPBE/RPBE functionals (designed for chemistry applications). 58 Recently, the meta-GGA SCAN functional 59 was also included in QuantumATK, often providing improved accuracy of DFT calculations as compared to PBE.

Hybrid functionals
Hybrid XC approximations mix local and/or semilocal functionals with some amount of exact exchange in order to provide higher accuracy for electronic-structure calculations. 52,60 However, the computational cost is usually much higher than for semilocal approximations. New methodological developments based on the adaptively compressed exchange operator (ACE) method 61 allow reducing the computational burden of hybrid functionals. The ACE algorithm was recently implemented in Quan-tumATK for HSE06 calculations, which gives a systematically good description of the band gap of most semiconductors and insulators, see Table III.

Semiempirical methods
Using hybrid functionals is computationally demanding for simulating large systems, often even prohibitive. QuantumATK offers a number of semiempirical XC methods that allow for computationally efficient simulations while giving rather accurate semiconductor band gaps. These include the DFT-1/2 method, 62,63 the TB09 XC potential, 64 and the pseudopotential projector-shift approach of Ref. 19.
The selfconsistent DFT-1/2 methods, including LDA-1/2 and GGA-1/2, do contain empirical parameters. In QuantumATK, these parameters are chosen by fitting the calculated band gaps to measured ones for bulk crystals. Table III suggests that the DFT-1/2 method, as implemented in QuantumATK, allows for significantly improved band gaps at almost no extra computational cost. We note that a recent study has shown certain limitations of the DFT-1/2 method, in particular for antiferromagnetic transition metal oxides. 65 Furthermore, this method does not provide reliable force and stress calculations. It is also important to note that not all species in the system necessarily require the DFT-1/2 correction. In general, it is advisable to apply this correction to the anionic species only, keeping the cationic species as normal. 62, 63 The Tran-Blaha meta-GGA XC functional (TB09) 64 introduces a parameter, c, which can be calculated self-consistently according to an empirical formula given in Ref. 64. Table III includes band gaps computed using this approach. The c-parameter may also be adjusted to obtain a particular band gap for a given material, and QuantumATK allows for setting different TB09 cparameters on different regions in the simulation cell. This may be useful for studying electronic effects at interfaces between dissimilar materials, for example in oxide-semiconductor junctions, where the appropriate (and material-dependent) c-parameter may be significantly different in the oxide and in the semiconductor.
QuantumATK also offers a pseudopotential projectorshift (PPS) method, that introduces empirical shifts of the nonlocal projectors in the PPs, in spirit of the empirical PPs proposed by Zunger and co-workers. 66 The PPS method is usually combined with ordinary PBE calculations. 19 The two main advantages of this PPS-PBE approach are that (1) for each semiconductor, the projector shifts can be fitted such that the DFTpredicted fundamental band gap and lattice parameters are both fairly accurate compared to measured ones, and (2) the PPS method does yield first-principles forces and stress, and therefore can be used for geometry optimization, unlike the DFT-1/2 and TB09 methods. Table IV shows that the PPS-PBE predicted equilibrium lattice parameters are only slightly overestimated, and the PPS-PBE band gaps are fairly close to experiments. We note that the PPS-PBE parameters are currently available in QuantumATK for the elements silicon and germanium only.

DFT+U methods
QuantumATK supports the mean-field Hubbard-U correction by Dudarev et al. 71 and Cococcioni et al., 72 denoted DFT+U, LDA+U, GGA+U, or XC+U. This method aims to include the strong on-site Coulomb interaction of localized electrons (often localized d and f electrons), which are not correctly described by LDA or GGA. A Hubbard-like term is added to the XC functional, where n l is the projection onto an atomic shell l, and U l is the Hubbard U for that shell. The energy term E U is zero for a fully occupied or unoccupied shell, but positive for a fractionally occupied shell. This favors localization of electrons in the shell l, typically increasing the band gap of semiconductors.

E. Boundary Conditions and Poisson Solvers
As already mentioned in Section IV A, the electron density, n(r) in (9), and the effective potential, V eff (r) in (3), are in QuantumATK represented on a real-space regular grid. The corresponding Hartree potential V H (r) is then calculated by solving the Poisson equation on this grid with appropriate boundary conditions (BCs) imposed on the six facets of the simulation cell, where e is the elementary charge, and 0 is the vacuum permittivity.
In QuantumATK, one may also specify metallic or dielectric continuum regions in combination with a micro-  (14) in the following way. For a metallic region denoted Ω, the electrostatic potential is fixed to a constant potential value (V 0 ) within this region, i.e., the Poisson equation is solved with the constraint For a dielectric region denoted Υ, the right-hand side of the Poisson equation will be modified as follows: where r is the relative dielectric constant, which can be specified as an external parameter in QuantumATK calculations.

Boundary conditions
QuantumATK implements four basic types of BCs; multipole, periodic, Dirichlet and Neumann BCs. It is also possible to impose mixed BCs on the six facets of the simulation cell to simulate a large variety of physical systems at different levels of approximation.
A multipole BC means that the Hartree potential at the boundary is determined by calculating the monopole, dipole and quadrupole moments of the charge distribution inside the simulation cell, and that these moments are used to extrapolate the value of the potential at the boundary. A Dirichlet BC means that the potential has been fixed to a certain potential V 0 (r) at the boundary, such that, for a facet S of the simulation cell, A Neumann BC means that the normal derivative of the potential on a facet has been fixed to a given function V 0 (r), where n denotes the normal vector of the facet. Next, we briefly describe applications of the different BCs.
• Multipole BCs are used for molecule configurations, ensuring the correct asymptotic behavior of the Hartree potential, even for charged systems (ions or charged molecules), as shown in Fig. 4(a).
• Periodic BCs is the natural choice along all directions for fully periodic bulk materials, as shown in Fig. 4(b). Periodic BCs are also often used to model heterostructures or interfaces, as well as surfaces using a slab model.
• Dirichlet-Neumann BCs for a slab model. In slab calculations, it can be more advantageous to impose mixed BCs, such as Neumann (fixed potential gradient) and Dirichlet (fixed potential) on the leftand right-hand side of the slab, respectively, combined with periodic BCs in the in-plane-directions, as shown in Fig. 4(c). These mixed BCs provide a physically sound alternative to the often-used dipole correction for slab calculations. 73 • Dirichlet-Neumann BCs for a surface configuration. For accurate surface simulations, the surface configuration may be used, in combination with mixed BCs: Neumann in the right-handside vacuum region, Dirichlet at the left electrode, and periodic BCs in the in-plane directions, see Fig. 4(d). In this case, one can account, e.g., for the charge transfer from the near-surface region to the semi-infinite electrode, which acts as an electron reservoir. 19 • Dirichlet BCs for a device configuration. Twoprobe device simulations are in QuantumATK done using Dirichlet BCs at the left and right boundaries to the electrodes. Periodic BCs may then be applied in the directions perpendicular to the electron transport direction, as shown in Fig. 4(e). For complex devices, one may need to apply a more mixed set of BCs, as discussed in the following.
• General mixed BCs. QuantumATK also allows for combining Neumann, Dirichlet and periodic BCs. This can be used to, e.g., model a 2D device in a field-effect transistor setup, such as that in Fig. 14.
We note that for systems with periodic or Neumann BCs in all directions, the Hartree potential can only be determined up to an additive constant. In this case, in order to obtain a uniquely defined solution, we require the average of the Hartree potential to be zero when solving the Poisson equation.

Poisson solvers
To handle such different BCs, the QuantumATK simulation engines use Poisson solvers based on either FFT methods or real-space finite-difference (FD) methods. The FD methods are implemented using a multigrid TABLE V. Classes of TB models currently supported by ATK-SE. The model types are either two-center Slater-Koster (SK) or based on environment-dependent parameters (Env). The model may be orthogonal (H) or non-orthogonal (H, S). Short-ranged models include nearest-neighbour interactions only (range up to a fewÅ), while the long-ranged Hückel models have a typical range of 5-10Å. As indicated in the right-hand column, not all models support calculation of total energies, forces, and stress, but are used mainly for simulating the electronic structure of materials.

Model
Ref solver, 74 a parallel conjugate-gradient-based solver, 75 and the MUMPS direct solver. 76 The real-space methods also allow for specifying spatial regions with specific dielectric constants or values of the electrostatic potential, as mentioned above. For systems with 2D or 3D periodic BCs, and no dielectric regions or metallic gates, the Poisson equation (14) is most efficiently solved using the FFT solvers. For a bulk configuration with 3D periodic directions, we use a 3D-FFT method, see Fig. 4(b). In the case of only 2 periodic directions, for example in slab models, surface configurations, and device configurations, we use a 2D-FFT method combined with a 1D finite-difference method, see Figs. 4(c)-(d). 77

V. SEMI-EMPIRICAL MODELS
As a computationally fast alternative to DFT, the ATK-SE engine allows for semi-empirical TBtype simulations. 24 The TB models consist of a nonselfconsistent Hamiltonian that can be extended with a selfconsistent correction for charge fluctuations and spin polarization. These corrections closely follow the density functional tight-binding (DFTB) approach. 78 The main aspects of these TB models have been described in Ref. 24 and below we give only a brief description of the models. Table V summarizes the available models for the nonselfconsistent part of the SE Hamiltonian, H 0 ij . Most of the models are non-orthogonal, that is, include a parametrization of the overlap matrix S ij . In most of the models, the Hamiltonian matrix elements depend only on two centers, parameterized in terms of Slater-Koster parameters. These models include Hückel models, 79,82 Slater-Koster orthogonal TB models, 18,83 and DFTB models. 78 The ATK-SE engine also supports models that take into account the position of atoms around the two centers. These currently include the environmentdependent TB models from Purdue 80 and those from the U.S. Naval Research Laboratory. 81 It is possible to add a selfconsistent correction to the non-selfconsistent TB models. 24 The selfconsistent correction use the change in the onsite Mulliken population of each orbital, relative to a reference system, to assign an orbital-dependent charge to each atom. The charge on the orbital is represented by a Gaussian orbital, and the width of the Gaussian, σ l , can be related to an onsite repulsion, U l , where l is the angular momentum of the orbital. The relation is given by 24 This onsite repulsion can be calculated from the charge-dependent onsite energies, 78 where ε l is the orbital energy of the atom and n l the charge in orbital l. QuantumATK comes with a database of U l calculated using DFT all-electron simulations of the atom. In practice, it is more reliable for each element to use a single averaged value, 78 where the average is determined by the number of valence electrons of each orbital, n l ; N = l n l . The ATK-SE default is to use such a single value.
In the ATK-SE selfconsistent loop, the Mulliken population is calculated for each orbital. Based on the change in charge relative to the reference system, a Gaussian charge is added at the orbital position. We note that in the default case, where an atom-averaged U is used on each orbital, only changes in the atomic charge will have an affect. From the atom-centered charge we set up a real-space charge density from which the Hartree potential V (r) is calculated using the same methods as used for DFT, see Section IV E. It is added to the TB Hamiltonian through where r i is the position of orbital i. The ATK-SE engine also supports spin polarization through the term 84 where the sign depends on the spin. The spin splitting of shell l, dE li , is calculated from the spin-dependent Mulliken populations m l↑ , m l↓ of each shell at the local site µ l : The shell-dependent spin-splitting strength W ll is calculated from a spin-polarized atomic calculation, 84 and ATK-SE provides a database with the parameters. The main advantage of the SE models compared to DFT methods are their computational efficiency. For large systems, the main computational cost of both DFT and TB simulations is related to diagonalization of the Hamiltonian, the speed of which depends strongly on the number of orbitals on each site and their range. This makes TB Hamiltonians very attractive for large systems, provided the SE parametrization is appropriate for the particular simulation. Furthermore, orthogonal Hamiltonians have inherent performance advantages. The Empirical and Purdue environment-dependent models are the most popular TB models for electron transport calculations. We also note that for many two-probe device systems, it is mainly the band structure and quantum confinement that determine electrical characteristics such as current-voltage curves. TB model Hamiltonians can provide good results for such simple devices. Finally, DFTB models are popular for total-energy calculations, although we find in general that the accuracy should be cross-checked against DFT.

VI. EMPIRICAL FORCE FIELDS
ATK-ForceField is a state-of-the-art FF simulation engine that is fully integrated into the Python framework. This has already been described in detail in Ref. 28, and we therefore only summarize some of the main features. Table VI lists the empirical potential models supported by ATK-ForceField, which includes all major FF types. The simulation engine also allows for combining models, such that different FFs can be assigned to different subsystems. The empirical potential for each sub-system, and the interactions between them, can be customized as desired, again using Python scripting. ATK-ForceField currently includes more than 300 predefined literature parameter sets, which can conveniently be invoked from the NanoLab GUI. Additionally, it is also possible to specify custom FF parameters via the Potential Editor tool in NanoLab or in a Python script, or even use built-in Python optimization modules to optimize the parameters against reference data.
Table VII compares the computational speed of ATK-ForceField molecular dynamics simulations to that of the popular LAMMPS package. 97 For most of the FF potential types, the two codes have similar performance.

VII. ION DYNAMICS
One very powerful feature of QuantumATK is that ion dynamics is executed using common modules that are not specific to the chosen simulation engine. This means that modules for calculating energy, forces, and stress may be used with any of the supported engines, including DFT, SE methods, and classical FFs. Options for ion-dynamics simulations are defined using Python scripting, which allows for easy customization, extension, and combination of different simulation methods, without loss of performance. In Section XIV C we illustrate this by combining the DFT and FF engines in a single molecular dynamics simulation. Several methods related to ion dynamics in QuantumATK have been described in detail in Ref. 28, so here we only summarize the main features.

A. Local Structural Optimization
The atomic positions in molecules and clusters are optimized by minimizing the forces, while for periodic crystals, the unit-cell vectors can also be included in the optimization, possibly under an external pressure that may be anisotropic. The simultaneous optimization of positions and cell vectors is based on Ref. 98, where the changes to the system are described as a combined vector of atomic and strain coordinates.

B. Global Structural Optimization
The previous section considered methods for local geometry optimization, which locate the closest local minimum-energy configuration. However, often the goal is to find the globally most stable configuration, for example, the minimum-energy crystal structure. Quantu-mATK therefore implements a genetic algorithm for crystal structure prediction. It works by generating an initial set of random configurations and then evolving them using genetic operators, as described in Ref. 101. An alternative approach is to perform simulated annealing using molecular dynamics. 102

C. Reaction Pathways and Transition States
The minimum-energy path (MEP) for changes to the atomic positions from one stable configuration to another may be found using the nudged elastic band (NEB) method. 103 The QuantumATK platform implements state-of-the-art NEB, 104 including the climbingimage method. 105 The initial set of images are obtained from linear interpolation between the NEB end points, or by using the image-dependent pair potential (IDPP) method. 106 The IDPP method aims to avoid unphysical starting guesses, and leads in general to an initial NEB path that is closer to the (unknown) MEP. This typically reduces the number of required NEB optimization steps by a factor of 2. In some implementations, the projected NEB forces for each image are optimized independently. However, the L-BFGS algorithm is in that case known to behave poorly. 107 In QuantumATK, the NEB forces for each image are combined into a single vector, F NEB ∈ R 3mn , where m is the number of images and n the number of atoms. This combined approach is more efficient when used with L-BFGS, and has been referred to as the global L-BFGS method. 107

D. Molecular Dynamics
Molecular dynamics (MD) simulations provide insights into dynamic atomic-scale processes or sample microscopic ensembles. The essential functional blocks in a typical QuantumATK MD loop are depicted in Fig. 5. Different thermodynamic ensembles can be simulated. The basic NVE ensemble uses the well-known velocity-Verlet algorithm. 108 Additionally, thermostats or barostats can be applied to different parts of the system to simulate NVT or NPT ensembles, using for example the chained Nosé-Hoover thermostat, 109   vectors, velocities, etc. This is often used to implement custom constraints on atoms or to apply a non-trivial strain to the simulation cell. The post-step hook is typically used to modify the forces and/or stress. It may, for example, be used to add external forces and stress contributions, such as a bias potential, to the regular interaction forces.
QuantumATK is shipped with a number of predefined hook functions, implementing thermal transport via reverse non-equilibrium molecular dynamics (RNEMD), 112 metadynamics, and other methods. For metadynamics, QuantumATK integrates with the PLUMED package, 113 so that all methods implemented in PLUMED are available in QuantumATK as well. Figure 6 illustrates the free-energy map of surface vacancy diffusion on Cu(111) using the ATK-ForceField engine with an EAM potential. 114

E. Adaptive Kinetic Monte Carlo
Adaptive kinetic Monte Carlo (AKMC) is an algorithm for modelling the long-timescale kinetics of solid-state materials. [115][116][117] For a given configuration, AKMC involves 3 steps: (1) locate all kinetically relevant product states; (2) determine the saddle point between the reactant and product states; (3) select a reaction using kinetic Monte Carlo (KMC).
Step 1 is in QuantumATK performed using hightemperature MD. At regular intervals, the MD simulation is stopped and a geometry optimization is performed to check if the system has left the initial basin. This procedure is repeated until all relevant reactions are found within a user-specified confidence. 117,118 In step 2, the saddle-point geometry for each reaction is determined by performing a NEB optimization for each reaction, and the reaction rates k are determined via harmonic transition-state theory (HTST), 119 (25) where N is the number of atoms, ν min i and ν ‡ i are the positive (stable) normal-mode frequencies at the minimum and saddle points, E min and E ‡ the corresponding energies, k B is the Boltzmann constant, and T the temperature. The ratio of the products of the vibrational frequencies in (25) is often called the attempt frequency or the prefactor, and can be computationally expensive to obtain. Instead of calculating the prefactor for each reaction mechanism, a user-given value may therefore be used.
Finally, in step 3, a reaction is selected using KMC, the system evolves to the corresponding product configuration, and the entire procedure is repeated. More details of the QuantumATK implementation of AKMC may be found in Ref. 117.

VIII. PHONONS
The ground-state vibrational motion of atoms is of paramount interest in modern materials science. Within the harmonic approximation, which is valid for small thermal displacements of atoms around their equilibrium position, the vibrational frequencies of a configuration are eigenvalues of the dynamical matrix D, where m a (m b ) is the atomic mass of atom a (b) and dF b,β /dr a,α is the force constant. Computing and diagonalizing D yields the vibrational modes of the system (molecular or periodic), and is also used to obtain the phonon density of states for a periodic crystal.

A. Calculating the Dynamical Matrix
QuantumATK calculates the dynamical matrix using a FD method, where each matrix element in (26) is computed by displacing atom a along Cartesian direction α, and then calculating the resulting forces on atom b along directions β. This approach is sometimes referred to as the frozen-phonon or supercell method, and applies equally well to isolated (molecular) systems. The method lends itself to heavy computational parallelization over many computing cores, since all displacements may be calculated independently. Crystal symmetries are taken into account in that only symmetrically unique atoms in the unit cell are displaced, and the forces resulting from displacement of the equivalent atoms are obtained using the corresponding symmetry operations. 120

B. Wigner-Seitz Method
For crystals with small unit cells, periodic repetition of the cell is usually needed to accurately account for long-range interactions in D. For larger simulation cells, including cells with defects and amorphous structures, this is not always necessary, since the cell may already include the entire interaction range. In order to recover the correct phonon dispersion across periodic boundaries, the Wigner-Seitz method can be employed. Here, a Wigner-Seitz cell is centered around the displaced atom and the forces on each atom in the simulation cell is assigned to its periodic image that is located within the Wigner-Seitz cell. 121

C. Phonon Band Structure and Density of States
The phonon band structure (or phonon dispersion) consists of bands with index λ of vibrational frequencies ω = ω λq throughout the Brillouin zone (BZ) of phonon wave vectors q. The phonon density of states (phonon DOS) per unit cell, g(ω), is defined as where N is the number of q-points in the sum. In practice, the phonon DOS is calculated using the tetrahedron method. 122 Additionally, quantities such as vibrational free energy, entropy, and zero-point energy can easily be calculated from the vibrational modes and energies. Figure 7 gives an example of phonon simulations for different metals using the ATK-ForceField and ATK-LCAO engines. The ATK-LCAO supercell calculation yields accurate vibrational properties, as exemplified by the excellent agreement between the two methods. The dispersions follow the same trends, which is expected, since the three metals have the same FCC crystal symmetry. We note that the higher phonon frequencies in Cu can be understood from the similar bond strength as in Ag and Au, but significantly lower Cu atomic mass.

D. Electron-Phonon Coupling
The electron-phonon coupling (EPC) is an imporant quantity in modern electronic-structure theory. It is, for example, used to calculate the transport coefficients in bulk crystals (see Section VIII E) and inelastic scattering of electrons in two-probe devices (see Section XI H).
To obtain the EPC, we calculate the derivative of the Hamiltonian matrix with respect to the position of atom a, r a , where ∂Ĥ/∂r a is calculated using finite differences, similar to the calculation of the dynamical matrix described above. A unit cell is repeated to form a supercell (for a device configuration, only the atoms in the central region are displaced). The terms that contribute to the Hamiltonian derivative is the local and non-local PP terms. The real-space Hamiltonian matrix is expanded in electron eigenstates, nk, and Fourier transformed using the phonon polarization vectors, to finally obtain the electron-phonon couplings g, g λnn kk q = n k |δĤ λq |nk , where q is the phonon momentum and λ the phonon branch index. Further details of how QuantumATK calculates the EPC are given in Ref. 123.

E. Transport Coefficients
The electron/hole mobility in a semiconductor material is an important quantity in device engineering, and also determines the conductivity of metals. Electronic transport coefficients for bulk materials, including the conductivity, Hall conductivity, and thermoelectric response, may be calculated from the Boltzmann transport equation (BTE) as linear-response coefficients related to the application of an electric field, magnetic field, or temperature gradient. In QuantumATK, this is done by expanding the current density j to lowest order in the electric field E, magnetic field B, and temperature gradient ∇T , where the indices label Cartesian directions and σ αβ , σ αβγ and ν αβ are the electronic conductivity, Hall conductivity, and thermoelectric response, respectively. Following Ref. 124, the band-dependent thermoelectric transport coefficients and Hall coefficients are obtained as where µ is the chemical potential and γuv the Levi-Civita symbol. The band group velocities v(nk) and effective mass tensors M(nk) are obtained from perturbation theory. Importantly, we may in (31) include the full scattering rate τ nk , and thereby go beyond the constant scattering-rate approximation used in Ref. 124. As we will see in Section XIV B, this may not only be important in order to obtain quantitatively correct results; it is also required to reproduce experimental trends in the conductivity of different materials.
The scattering rate is given by where B is a temperature-dependent scattering weight, where f is the Fermi function, and the scattering angle is defined by Furthermore, P (P ) are transition rates due to phonon absorption (emission). They are obtained from Fermi's golden rule, where n λq is the phonon occupation operator, and g λnn kk q the EPC constant from (29). QuantumATK offers two different methods for performing the q-integral in (32). In the first method, the delta functions in (35) are represented by Gaussians with a certain width, and we perform the discrete sum over q.
In the second method, we realize that the integral closely resembles the numerical problem of obtaining a density of states, and use the tetrahedron method 122 for the integration. In particular for metals, we find the tetrahedron method to be most efficient. Figure 8(c) shows the convergence of the Au resistivity as the number of q-points increases, using both Gaussian and tetrahedron integration. The tetrahedron calculation seems converged for In semiconductors, it is possible to make a clever selection of kand q-points to minimize the computational load while including all relevant scattering processes. Typically, a sparse k-point sampling is used for the mobility integral, while a denser q-point sampling is needed to secure a correct scattering rate at each k-point.
(b) Fermi-surface of bulk Au. In metals, k-points contributing to the neighborhood of the Fermi surface are not located in a small subset of the Brillouin zone. Therefore, kand q-points are sampled in the full Brillouin zone, and q-space is integrated using the tetrahedron method to minimize the sampling density. (c) Resistivity convergence with respect to the number of q-points for bulk Au with either a direct or tetrahedron integration over the full Brillouin zone. Resistivities were calculated with N k × N k × N k k-points for a sequence of Nq × Nq × Nq q-points.
N q = 20, that is, a 20 × 20 × 20 q-point sampling. The result with a finite Gaussian broadening may converge fast if using a rather large broadening, but the resistivity then appears to converge to a wrong result. In general, we therefore recommend the tetrahedron integration method for calculation of metallic resistivity.
To further improve the computational performance when calculating transport coefficients, it is possible to use the energy-dependent isotropic-scattering-rate approximation, introduced in Ref. 125. A two-step procedure is used for the k-point sampling, which significantly reduces simulation time without affecting the resulting mobilities for many materials (those that have a fairly isotropic scattering rate in momentum space). In step one, an initial k-space with a low sampling density and a well-converged q-point sampling are used. The initial kpoint grid is automatically reduced further by including only k-points where the band structure has energies in a specific range around the Fermi level. This limits the simulations to the relevant range of initial states (and relevant carrier densities), which significantly increases simulation speed and reduces memory usage. Typically, the variation of the scattering rates from the different directions in momentum space will be small. Fom the obtained data, we may therefore generate an isotropic scattering rate that only depends on energy, where we have integrated over bands n and wave vectors k, and n(E) is the density of states. In the second step, we then perform a calculation on a fine k-point grid, but using the energy-dependent isotropic scattering rate τ (E). Since the scattering rate often varies slowly on the Fermi surface (for metals), this is a good approximation. The second step therefore requires only an evaluation of band velocities and effective masses on the dense k-point grid, while the scattering rate is reused. This two-step procedure, combined with either direct integration for semiconductors and semimetals, or tetrahedron integration for metals, makes QuantumATK an efficient platform for simulating phonon-limited mobilities of materials. In addition, it is possible to input a predefined scattering rate as a function of energy. This is relevant for adding extra scattering mechanisms, for example impurity scattering, on top of the electron-phonon scattering, or in the case where a scattering-rate expression is known analytically. One special case of the last situation is the limit of a constant relaxation time, which is the basis of the popular Boltztrap code. 124 We note that such constant-relaxation-time calculations are easily performed within the more general QuantumATK framework outlined above. Moreover, since electron velocities are calculated from perturbation theory, accuracy is not lost due to band crossings, which is the case when velocities are obtained from FD methods, as is done in Ref. 124. In some cases, the constant relaxation time approximation can give a good first estimate of thermoelectric parameters for a rough screening of materials, but for quantitative predictions, the more accurate models of the relaxation time outlined above must be used.

IX. POLARIZATION AND BERRY PHASE
Electronic polarization in materials has significant interest, for example in ferroelectrics, where the electric polarization P can be controlled by application of an external electric field, or in piezoelectrics, where charge accumulates in response to an applied mechanical stress or strain. 126 It is common to divide the polarization into ionic and electronic parts, P = P i + P e . The ionic part can be treated as a classical electrostatic sum of point charges, where Z ion a and r a are the valence charge and position vector of atom a, Ω is the unit-cell volume, and the sum runs over all ions in the unit cell.
The electronic contribution to the polarization in direction α is obtained as 126 where R α is the lattice vector in direction α, and the Berry phase Φ α is obtained as where the sum runs over N ⊥ k ⊥ -points in the BZ plane perpendicular to R α , and with the overlap integrals and with the J k-points given by k j = k ⊥ + k ,j lying on a line along the R α direction. The polarization depends on the coordinate system chosen since it is related to the real-space charge position, and is determined by the Berry phase, which is only defined modulo 2π. Consequently, the polarization is a periodic function and constitutes a polarization lattice itself. The polarization lattice in direction α is written as where n is an integer labeling a polarization branch, and the polarization quantum in direction α is P Q,α = |e| Ω R α . All measurable quantities are related to changes in the polarization, which is a uniquely defined variable, provided that the different polarization values are calculated for the same branch in the polarization lattice.
QuantumATK supports calculation of the polarization itself, as well as the derived quantities piezoelectric tensor, where Voigt notation is used for the strain component, that is, i ∈ (xx, yy, zz, yz, xz, xy), and the Born effective charge tensor where the derivative is with respect to the position of atom a in direction β. Table VIII shows calculated values of the Born effective charges (only the negative components for each structure) and elements of the piezoelectric tensor for III-V wurtzite nitrides and zincblende GaAs. The calculated Born effective charges and piezoelectric tensor components agree well with the reference calculations.

X. MAGNETIC ANISOTROPY ENERGY
The magnetic anisotropy energy (MAE) is an important quantity in spintronic magnetic devices. The MAE is defined as the energy difference between two spin orientations, often referred to as in-plane ( ) and out-of-plane (⊥) with respect to a crystal plane of atoms, a surface, or an interface between two materials: The MAE can be split into two contributions: A classical dipole-dipole interaction resulting in the so-called shape anisotropy, and a quantum mechanical contribution often refered to as the magnetocrystalline anisotropy, which arises as a consequence of spin-orbit coupling (SOC). In this section we will focus on the magnetocrystalline anisotropy and refer to this as the MAE. There are at least three different ways of calculating the MAE: (i) Selfconsistent total-energy calculations including SOC with the noncollinear spins constrained in the in-plane and out-of-plane directions, respectively, (ii) using the force theorem (FT) to perform non-selfconsistent calculations (including SOC) of the band-energy difference induced by rotating the noncollinear spin from the in-plane to the out-of-plane direction, and (iii) second-order perturbation theory (2PT) using constant values for the SOC. While it has been demonstrated that methods (i) and (ii) give very similar results, 128,129 the 2PT method can lead to significantly different results. 129 In QuantumATK we have implemented an easy-to-use workflow implementing the FT method (ii). Using the FT gives the advantage over method (i) that the calculated MAE can be decomposed into contributions from individual atoms or orbitals, which may give valuable physical and chemical insight.
The QuantumATK workflow for calculating the MAE using the FT method is the following: 1. Perform a selfconsistent spin-polarized calculation.
2. For each of the considered spin orientations (a) Perform a non-selfconsistent calculation, in a noncollinear spin representation including SOC, using the effective potential and electron density from the polarized calculation but rotated to the specified spin direction.
(b) Calculate the band energies n and projection weights w n,p .

Calculate the total MAE as
where f n is the occupation factor for band n (including both band and k-point index) for the spin orientation and n is the corresponding band energy, and likewise for the ⊥ spin orientation.
The contribution to the total MAE for a particular projection p (atom or orbital projection) is where the projection weight is w n,p = ψ n |(SP + PS)/2|ψ n , with |ψ n being the eigenstate, S the overlap matrix, and P the projection matrix. P is a diagonal, singular matrix with ones in the indices corresponding to the orbitals we wish to project onto and zeros elsewhere.  129. We first note that the calculated MAEs agree rather well among the four codes, the only exception being FeAu, where the LCAO representations give somewhat smaller values than obtained with PW expansions. In this case it seems that the LCAO basis set has insufficient accuracy, which could be related to the fact that the LCAO basis functions are generated for a scalar-relativistic PP derived from a fully relativistic pseudopotential. Figure 9 shows the atom-and orbital-projected MAE for a Fe/MgO interface. The structure is similar to the one reported in Ref. 130. We use periodic BCs in the transverse directions. The calculated interfacial anisotropy constant K 1 = MAE/(2A), where A is  > 0), while the atoms in the center of the Fe slab contribute with much smaller values. From the orbital projections it is evident that the MAE peak at the interface is caused primarily by a transition from negative to positive MAE contributaions from the Fe d xy and d x 2 −y 2 orbitals, which hybridize with the nearby oxygen atom.

XI. QUANTUM TRANSPORT
The signature feature of QuantumATK is simulation of device systems. While most DFT device simulation codes are constructed on top of an electronic structure code designed for simulating bulk systems, QuantumATK is designed from scratch to achieve the highest accuracy and performance for both bulk and device systems. Figure 10 shows a device (two-probe) geometry. It consists of a left electrode, a central region, and a right electrode. The three regions have the same BCs in the two lateral directions perpendicular to the left-right electron transport direction, as defined in Fig. 10. The left and right electrodes are assumed to have bulk properties, and the first step of the device simulation is to perform a bulk calculation of each electrode with periodic BCs in the transport direction. Using Bloch's theorem, we describe the wave functions in terms of transverse k-points, and to seamlessly connect the three regions, the same kpoint sampling is used in the transverse directions for all three regions. In the transport direction, the centralregion wave functions are described by using scattering BCs, while the electrode wave functions are described by using periodic BCs. To have a seamless connection, it is important that the electrode wave functions very accurately reproduce the infinite-crystal limit in the transport direction. A very dense electrode k-point grid is therefore needed in the transport direction.
The left and right electrodes are modelled in their ground states with chemical potentials µ L and µ R , respectively. This is only a correct model if the electrodes are not affected by the contact with the central region. The central-region electrostatic potential should therefore be sufficiently screened in the regions interfacing with the electrodes (denoted "electrode extensions"), such that the potential in each electrode extension virtually coincides with that in the electrode. Furthermore, the approximation is not valid if the finite-bias current density is high; in this case a non-equilibrium electron occupation is needed to accurately model the electrodes. A device with no electron scattering in the central region can therefore not be modelled reliably at finite bias.
The electronic structures of the isolated electrodes are defined with respect to an arbitrary energy reference. When used in a device simulation, they must be properly aligned to a common reference. This is achieved by applying a potential shift to the electronic structure of the right electrode, chosen to fulfill the condition where V bias is the bias applied on the electrodes. It is clear that µ R = µ L at zero bias. The electrode electrostatic potentials, including the right-electrode potential shift, sets up the BCs for the central-region electrostatic potential. Thus, the whole system is aligned to a common reference, and device built-in potentials, if any, are properly included.  ) have an equilibrium electron distribution with chemical potentials µL and µR, related through the applied sample bias, µR − µL = eV bias . At T = 0 K, the electrons with energies in the bias window, µL ≤ ε ≤ µR, give rise to a steady-state electrical current from the right to left electrode. Note that the electron transport direction is from the left to right electrode. For higher temperatures, the electrons above (below) µL (µL) will also contribute to the current because of the corresponding broadening of the Fermi-Dirac distribution at T > 0 K. The system is modelled selfconsistently at the DFT or TB level using the NEGF method. It is possible to include the effect of gate potentials in the selfconsistent solution. Inelastic effects due to phonon or photon scattering can be included through perturbation theory.
The electrostatic potential enters the KS equation from which the electron density in the central region is determined. We assume the system is in a steady state, that is, the central-region electron density does not change with time. The density can then be described in terms of extended electronic states from the left and right electrodes, as well as bound states in the central region, n(r) = n L (r) + n R (r) + n B (r).
We now focus on the contribution from the extended states of the left (n L ) and right (n R ) electrodes, and delay the discussion of bound states (n B ) for later. The former may be obtained by calculating the scattering states incoming from the left (ψ L α ) and right (ψ R α ) electrodes, which can be obtained by first calculating the Bloch states in the electrodes, and subsequently solving the KS equation for the central region using those Bloch states as matching BCs.
The left and right electron densities can then be cal-culated by summing up the occupied scattering states, where f (x) = (1 + e x ) −1 is the Fermi-Dirac distribution.

A. NEGF Method
Instead of using the scattering states to calculate the non-equilibrium electron density, QuantumATK uses the NEGF method; the two approaches are formally equivalent and give identical results. 38 The electron density is given in terms of the electron density matrix. We split the density matrix into left and right contributions, The left contribution is calculated using the NEGF method as 38 where is the spectral density matrix, expressed in terms of the retarded Green's function G and the broadening function Γ L of the left electrode, which is given by the left electrode self-energy Σ L . Note that while there is a non-equilibrium electron distribution in the central region, the electron distribution in the left electrode is described by a Fermi-Dirac distribution f with an electron temperature T L . Similar equations exist for the right density matrix contribution. The next section describes the calculation of G and Σ in more detail.
We note that the implemented NEGF method supports spintronic device simulations, using a noncollinear electronic spin representation, and possibly including spinorbit coupling. This enables, for example, studies of spintransfer torque driven device physics. 131

B. Retarded Green's Function
The NEGF key quantity to calculate is the retarded Green's function matrix for the central region. It is calculated from the central-region Hamiltonian matrix H and overlap matrix S by adding the electrode self-energies, where δ + is an infinitesimal positive number. Calculation of G at a specific energy ε requires inversion of the central-region Hamiltonian matrix. The latter is stored in a sparse format, and we only need the density matrix for the same sparsity pattern. This is done by block diagonal inversion, 132 which is O(N ) in the number of blocks along the diagonal.
The self-energies describe the effect of the electrode states on the electronic structure in the central region, and are calculated from the electrode Hamiltonians. QuantumATK provides a number of different methods, [133][134][135][136] where our preferred algorithm use the recursion method of Ref. 134, which in our implementation exploits the sparsity pattern of the electrode. This can greatly speed up the NEGF calculation as compared to using dense matrices.

C. Complex Contour Integration
The integral in (54) requires a dense set of energy points due to the rapid variation of the spectral density along the real axis. We therefore follow Ref. 38 and divide the integral into an equilibrium part, which can be integrated on a complex contour, and a non-equilibrium part, which needs to be integrated along the real axis, but only for energies within the bias window. We have where where ρ B is the density of states of any bound states in the central region. Equivalently, we could write the density matrix as where L and R are exchanged in (59) and (60). Due to the finite accuracy of the integration along the real axis, (58) and (61) are numerically different. We therefore use a double contour, 38 where (58) and (61) are weighted such that the main fraction of the integral is obtained from the equilibrium parts, D L eq and D R eq , which are usually much more accurate than the non-equilibrium parts, due to the use of high-precision contour integration. We have where W L and W R are chosen according to Ref. 38, i.e. such that at each site, the equilibrium part of the density matrix gives the largest contribution and W L + W R = 1.

D. Bound States
The non-equlibrium integrals, ∆ L neq and ∆ R neq , do not include any density from bound states in the central region. However, the equilibrium part of the density matrix is calculated from a complex contour integral of the retarded Green's function, and this calculation includes bound states with energies below the chemical potential of the contour.
Assume µ L < µ R , then a bound state with energy ε B < µ L will be included in both D L eq and D R eq , but a bound state in the bias window, µ L < ε B < µ R , will only be included in D R eq . Thus, from (62) we see that the state will be included with weight 1 if ε B < µ L and only with a fractional weight if µ L < ε B < µ R . The weight will depend on the position of the bound state along the transport direction, that is, if the bound state is in a region that is well connected with the right electrode, the occupation will follow the right electrode and thus be close to 1. If it is in a region that is not well connected with the right electrode, the occupation will follow the left electrode, and thus for the current example the occupation be close to 0.
The true occupation of a bound state in the bias window will depend on the physical mechanism responsible for the occupation and de-occupation, for example electron-phonon scattering, defects, etc. However, the matrix element will typically be higher with the electrode that is well connected with the region around the bound state, so we believe that the use of a double contour gives a qualitatively correct description of the occupation of the bound states in the bias window. Furthermore, we find that if we do not use such weighting schemes, bound states in the bias window can cause instabilities in the selfconsistent finite-bias NEGF calculation.

E. Spill-in Terms
Given the density matrix D, the electron density is obtained from the LCAO basis functions φ: The Green's function of the central region gives the density matrix of the central region, D CC . However, to calculate the density correctly close to the central-region boundaries towards the electrodes, the terms involving D LL , D LC , D CR , and D RR are also needed. These are denoted spill-in terms. 137 QuantumATK implements an accurate scheme for including all the spill-in terms, both for the electron density and for the Hamiltonian integrals. 137 This gives additional stability and well-behaved convergence in device simulations.

F. Device Total Energy and Forces
A two-probe device is an open system where charge can flow in and out of the central region through the left and right electrode reservoirs. Since the two reservoirs may have different chemical potentials, and the particle number from a reservoir is not conserved, it is necessary to use a grand canonical potential to describe the energetics of the system, 138 where N L/R is the number of electrons contributed to the central region from the left/right electrode, and E KS [n] is the KS total energy. Due to the screening approximation, the central region will be charge neutral, and therefore N L +N R = N , where N is the ionic charge in the central region. At zero bias (µ L = µ R ), the particle term is constant, so that N µ L = N µ R , and is thus independent of atom displacements in the central region. However, at finite bias (µ L = µ R ), the particle terms in Ω will affect the forces.
If one neglects current-induced forces, 139,140 as done in QuantumATK simulations, the force acting on atom a at position r a in the device central region is given by It can be shown that the calculation of this force is identical to the calculation of the equilibrium (zero-bias) force, but in the non-equilibrium (finite-bias) case the density and energy density matrix must be calculated within the NEGF framework. 38,138,141

G. Transmission Coefficient and Current
When the selfconsistent non-equilibrium density matrix has been obtained, it is possible to calculate various transport properties of the system. One of the most notable is the transmission spectrum from which the current and differential conductance are obtained. The transmission coefficient T at electron energy ε is obtained from the retarded Green's function, 142 and the electrical current is given by the Landauer formula,

H. Inelastic Transmission and Inelastic Current
QuantumATK implements the lowest-order expansion (LOE) method 143 for calculating the inelastic current due to electron-phonon scattering, which is not included in (66) and (67). The LOE method is based on perturbation theory in the first Born approximation, and requires calculation of the dynamical matrix and the Hamiltonian derivative with respect to atomic positions in the central region, ∇H(r). Calculation of these derivatives are described in Section VIII.
First-principles calculation of ∇H(r) can be prohibitive for large device systems. However, if the atomic configuration of the central region can be generated by repeating the left electrode along the transport direction, then ∇H(r) can be obtained to a good approximation by using the ∇H(r) of the left electrode only. 144 From ∇H(r) of the central region we get the electronphonon matrix elements in reciprocal space, 123 where the (mn)-sum runs over repeated unit cells in the supercell calculation of the Hamiltonian derivatives, 123 and the subscript 0 indicates that the derivatives are only calculated for atoms in the unit cell with index 0. Moreover, |φ i R n (|φ j R m ) denotes the i(j)'th LCAO basis orbital in the unit cell displaced from the reference cell by the lattice vector R n (R m ), while q is the phonon momentum, and v λ,q is the mass-scaled mode vector of phonon mode λ with frequency ω λ,q . Following Ref. 143, we obtain the inelastic transmission functions for a finite transfer of momentum. From these we calculate the total electrical current, including inelastic effects. 144, 145 The complete formulas for the QuantumATK implementation can be found in Ref. 144.

Special Thermal Displacement Method
In Ref. 146 we showed that the average transmission from a thermal distribution of configurations accurately describes the inelastic electron transmission spectrum due to electron-phonon scattering at this temperature. In the special thermal displacement (STD) method, the average is replaced with a single representative configuration, which may drastically reduce the computational cost of inelastic transport simulations. 147 To obtain the STD configuration, we first calculate the phonon eigenspectrum using the dynamical matrix of the central region. We consider only q = 0, since only relative displacements between atoms in the cell will be important, and to account for finite q-vectors we will have to increase the cell size. The phonon modes are labeled by λ with frequency ω λ , eigenmode vector e λ , and characteristic length l λ .
The STD vector of atomic displacements is given by 147 where s λ denotes the sign of the first non-zero element in e λ , enforcing the same choice of "gauge" for the modes. The Gaussian width σ is related to the mean square dis- An essential feature of the STD method is the use of opposite phases for phonons with similar frequencies; in this way phonon-phonon correlation functions average to zero and the transmission spectrum of the STD configuration becomes similar to a thermal average of single phonon excitations.
The final step in the STD method is to calculate the selfconsistent Hamiltonian of the system displaced by u STD , and use that to calculate the transmission spectrum. Thus, the computational cost of the inelastic transmission calculation is for the STD method similar to that of an ordinary elastic transmission calculation.
Formally, this method becomes accurate for systems where the central region is a large unit cell generated by the repetition of a basic unit cell.

I. Thermoelectric Transport
The thermoelectric figure of merit, ZT, quantifies how efficiently a temperature difference (heat) can be converted into a voltage difference in a thermoelectric material, where G e is the electronic conductance, S the Seebeck coefficient, T the temperature, and κ = κ e + κ ph the summed electron and phonon heat transport coefficients. Following Ref. 148, and given a set of electron and phonon transmission spectra for a device configuration, QuantumATK uses linear-response theory to compute the above-mentioned thermoelectric coefficients and the Peltier coefficient, Π, where I Q = dQ/dT is the electronic contribution to the heat current. It is calculated in a similar way as the electronic current, 149 where µ = (µ L + µ R )/2 is the average chemical potential, and the difference to (67) is the inclusion of the factor (ε − µ) in the integral. Note that one may use DFT or a TB model for obtaining the electron transmission and a force field to calculate the phonon transmission, constituting a computaionally efficient workflow for investigating thermoelectric materials.

J. Photocurrent
QuantumATK allows for calculating photocurrent using first-order perturbation theory within the first Born approximation. [150][151][152] In brief, the electron-light interaction is added to the Hamiltonian, whereĤ 0 is the Hamiltonian without the electron-light interaction, e the electron charge, m 0 the free-electron mass,p the momentum operator, and A ω the electromagnetic vector potential from a single-mode monocromatic light source with frequency ω. The first-order coupling matrix is where |j is an LCAO basis function. The first Born electron-photon self-energies are (79) where ε ± = ε± ω, and N is the number of photons. The Green's function including electron-photon interactions to first order is then where G r,>,< 0 denote the non-interacting Green's functions, and Σ >/< L,R are the lesser and greater self-energies due to coupling to the electrodes. The current in electrode α (left or right) with spin σ is calculated as where the effective transmission coefficients are given by 152 (82) We note that it is possible to include also the effect of phonons through the STD method, which is important for a good description of photocurrent in indirect-bandgap materials such as silicon. 153

XII. QUANTUMATK PARALLELIZATION
Atomic-scale simulations for small configurations (systems with only a few atoms) may often be executed in serial on a single CPU core, but most production simulations require execution in parallel on several cores (often many) to increase computational speed and/or to reduce the per-core memory footprint. The QuantumATK platform offers several parallelization techniques depending on the type of computational task. . Scaling performance of QuantumATK DFT simulations for a 64-atom Si0.5Ge0.5 random-alloy supercell when executed in parallel (using MPI) on 1, 2, 4, and 8 computing nodes (16 cores per node). a) Total wall-clock times for LCAO and PW selfconsistent total-energy calculations, and b) the corresponding peak memory requirements per core. Grey lines indicate ideal scaling of the wall-clock time. Pseu-doDojo PPs with LCAO-High basis sets were used. Note that the Ge PP contains semicore states. The supercell has 32 irreducible k-points, corresponding to two computing nodes for full MPI parallelization over k-points. With 4 (8) full nodes, 2 (4) MPI processes are assigned to eack k-point.
of, CPU cores, and also allows for assigning multiple processes to each work unit. Moreover, each MPI process may be further distributed in a hybrid parallelization scheme by employing shared-memory threading of each process. Figure 11 shows an example of how the total wall-clock time and peak memory requirement for DFT-LCAO and DFT-PW calculations scale with the number of 16-core computing nodes used with MPI parallelization. We 2)+*[EPPGPSGOXMQIQMRYXIW E 0MRIEVWGEPMRK 14-14-XLVIEHW 14-XLVIEHW 14-XLVIEHW GSQTYXMRKGSVIW 2)+*TIEOQIQSV]TIVGSVI+& 14-

14-XLVIEHW
14-XLVIEHW 14-XLVIEHW F FIG. 12. Scaling performance of equilibrium DFT-NEGF simulations for a 10 nm long silicon p-n junction with doping levels of 5 · 10 20 cm −3 . The junction cross section is 1.33 nm 2 , corresponding to 684 Si atoms in the device central region. The NEGF calculations were done using a PseudoDojo PP with the LCAO-Low basis set, and 2 irreducible k-points in the central-region 2D Brillouin zone, resulting in 96 generalized contour points. The simulations were run on up to eight 24-core Intel Xeon nodes, using both MPI (purple) and hybrid parallelization schemes. Hybrid parallelization was done using 2 (orange), 4 (green), and 24 (blue) threads per MPI process, with processes distributed evenly over the nodes. Gray dashed line indicates ideal scaling of the wall-clock time.
considered a 64-atom SiGe random-alloy supercell with N k = 32 k-points. In this case, 2 full nodes, N n = 2, with 32 cores in total (N c = N n × 16 = 32), yields full MPI parallelization over k-points. The PW calculations were done using a blocked generalized Davidson algorithm 39,155 to iteratively diagonalize the Hamiltonian matrix, which in the QuantumATK implementation parallelizes the computational work over both k- points and plane waves. The LCAO calculations use the LAPACK 156 (when N c /N k ≤ 1) or ELPA 157 (when N c /N k > 1) libraries to distribute Hamiltonian diagonalization over MPI processes. It is clear from Fig. 11 that the LCAO engine is both fast and requires less memory than the PW representation for the 64-atom supercell, although communication overhead causes the LCAO computational speed to start breaking off from ideal scaling when the number of processes (cores) exceeds the number of k-points in the DFT calculation (when N c /N k > 1). On the contrary, MPI parallelization over both k-points and plane waves enables approximately ideal scaling of the PW wall-clock time up to at least 8 nodes (128 cores), corresponding to 4 MPI processes per k-point.

B. DFT-NEGF Device Simulations
As discussed in Section XI C, the NEGF equilibrium density matrix at a single k-point is obtained from integrating the spectral density matrix over M ε energy points on a complex contour. This integral must be performed at all transverse k-points in the 2D Brillouin zone of the device central region, yielding N k × M ε generalized contour points. Each of these constitute a unit of computational work in equilibrium NEGF calculations, equivalent to k-point parallelization in DFT calculations for periodic bulks.
Since we typically have M ε = 48 contour energies, an equilibrium NEGF simulation may easily require evaluation of hundreds of generalized contour points. MPI parallelization over contour points is therefore a highly efficient strategy. For devices with relatively large trans-verse cross sections, and therefore relatively few contour points (because of small N k ), assignment of several processes to each contour point enables scaling of NEGF computational speed to numbers of computing cores well beyond the number of contour points. This can also be combined with more than one thread per process in a hybrid parallelization scheme, for a smaller speedup, but with a reduced per-core memory footprint. Figure 12 shows an example of how the total wall-clock time and peak memory usage for a DFT-NEGF calculation scale with the number of computing nodes used with both MPI and hybrid parallelization schemes. Calculations for this 10 nm long silicon p-n junction require evaluation of 96 generalized contour points, in this case corresponding to 4 nodes for full MPI distribution of computational work. As expected, we find that using only MPI parallelization requires most memory per core, but also results in the smallest wall-clock time for the NEGF calculation, although communication overhead causes a deviation from ideal scaling for more than 1 node, see Fig. 12(a). We also note that the per-core memory consumption is in this case almost constant in Fig. 12(b), except for a modest decrease for 8 nodes, where 2 processes (cores) are assigned to each contour point. It is furthermore clear from Fig. 12 that hybrid parallelization enables significant memory reduction, although at the cost of decreased computational speed. Taking simulation on 4 nodes as an example, hybrid parallelization with 4 threads per process (green lines) requires in this case 50% more wall-clock time as compared to the MPIonly simulation (purple lines), but at a 70% smaller memory footprint.
Although NEGF computational efficiency and memory consumption depend significantly on the device length and transverse dimensions, the general trend is that MPI parallelization over contour points yields computational speedup, while threading of processes reduce the NEGF memory footprint at a comparatively smaller computational speedup.

C. FF Simulations
The ATK-ForceField engine uses shared-memory threading for parallelization of relatively small systems, while additional parallelization by domain decomposition over MPI processes is available for large systems. As explained in detail in Ref. 28, the MPI distribution of ATK-ForceField workload is implemented via functionality from the Tremolo-X MD package, 158 which is developed by the Fraunhofer Institute for Algorithms and Scientific Calculations (SCAI).
In Fig. 13, we show the wall-clock time per MD step for a simulation of SiO 2 with 1 million atoms, using a force field from Pedone et al. 154 This illustrates how the use of domain decomposition over MPI processes results in a significant speedup when parallelizing over a large number of nodes and cores.

A. Python Scripting
The QuantumATK software is programmed in the C++ and Python languages. Around 80% of the code lines are in Python, and only low-level numerically demanding parts are written in C++. The use of Python allows for using a large number of high-level physics and mathematics libraries, and this has greatly helped building the rich functionality of QuantumATK in a relatively short time.
The user input file is a Python script and the user has through the script access to the same functionality as a QuantumATK developer. This enables the user to transform input files into advanced simulation scripts, which do not only set up advanced workflows and analysis, but may also alter the functionality of the simulation engines, for example by adding new total-energy terms. Quan-tumATK supplies a public application programming interface (API) with currently more than 350 classes and functions. These all take a number of arguments with detailed checks of the input parameters to ensure correct usage. For example, if the input argument is a physical quantity, the physical units must be supplied. A wide range of units are supported, e.g., for energy, the user may select units of joule, calorie, electron volt, kilojoule per mole, kilocalories per mole, Hartree, or Rydberg. All physical units are automatically converted to the internal units used by QuantumATK. The user also has access to internal quantities such as the Hamiltonian, Green's function, self-energies, etc., through the API. Through Python scripting it is possible to build advanced workflows that automate complex simulations and analysis. However, some simulations may require a large number of time consuming calculation tasks that are combined into a final result, and scripting such workflows can be impractical. For instance, if the computer crashes during a loop in the script, how to restart the script at the right step in a loop in the middle of the script? Or perhaps some additional tasks are needed after a custom simulation has finished; how to combine the already calculated data with the new data?
To simplify such simulations, QuantumATK has introduced a framework called a study object. The study object keeps track of complex simulations that rely on execution and combination of a number of basic tasks. It allows for running the basic tasks in parallel and will be able to resume if the calculation is terminated before completion. A study object also allows for subsequently extending the number of tasks, and will only perform tasks that have not already completed. This framework is currently used for a number of complex simulations, for instance for coupling atomic-scale simulations with continuum-level TCAD tools. Examples include simulation of the formation energy and diffusion paths of charged point defects, scans over source-drain and gate bias for two-terminal devices, relaxation of devices, and calculation of the dynamical matrix and Hamiltonian derivatives by finite differences.
To store data we use the cross-platform HDF5 binary format, 159 which allows for writing and reading data in parallel to/from a single file. This file can also hold many different objects, so the entire output from a Quantu-mATK simulation can be stored efficently in a single file.

B. NanoLab Graphical User Interface
While scripting is very efficient for production runs, it requires knowledge of the scripting language, and it takes time to manually build up scripts for setting up the configuration, simulation, and analysis of interest. The NanoLab GUI eliminates this barrier to productivity by enabling the user to fully set up the Python input script in a professional GUI environment. NanoLab is itself programmed in Python, and each tool in NanoLab can interpret and generate Python scripts, thus, it is possible to seamlessly shift from using the GUI tools in NanoLab to manually editing the Python scripts. It is the ambition that all NanoLab functions are also available as Python commands, such that any GUI workflow can be documented and reproduced in a Python script.
NanoLab is developed around a plugin concept, which makes it easy to extend it and add new functionality. Plugins can be downloaded and installed from an add-on server, and the majority of the plugins are available as source code, making it easy to modify or extend them with new user-defined functionality.
NanoLab also provides GUI tools for communicating with online databases ("Databases"), setting up the atomic-scale geometry of configurations ("Builder"), writing the Python script ("Scripter"), submitting the script to a remote or local computing unit ("Job Manager"), and visualizing and analyzing the results ("Viewer"). It is possible to connect third-party simulation cods with NanoLab by writing plugins that translate the input/output files into the internal NanoLab format. Such plugins are currently available for the VASP, 22 Quantum ESPRESSO, 23 ORCA, 160 GPAW, 161 and CASTEP 162 codes.
The plugin concept also allows for many specialized functions, for example specialized Builder tools like surface builders, interface builders, 163 NEB setups, 106 etc. The Job Manager has plugins that provide support for a wide range of job schedulers on remote computing clusters. Moreover, NanoLab has a large selection of graphical analysis tools, which can be used to visualize and analyze simulations with respect to a wide range of properties, all implemented as plugins. For instance, with the "MD analyser" plugin, a MD trajectory can be analyzed with respect to angular and radial distribution functions, or different spatial and time correlation functions. Other examples are interactive band structure analysis with extraction of effective masses, and analysis of transmission in device simulations with on-the-fly inspection of transmission eigenstates at specified points in the transmission spectrum. NanoLab currently ships with more than 100 preinstalled plugins, and additional plugins are available through the add-on server.

C. Documentation
Keeping an updated documentation system for the large set of QuantumATK classes and functions pose a challenge. To synchronize the documentation with the source code, we have developed an automated documentation system where the information for the Quantu-mATK reference manual is extracted directly from the Python source code using the Sphinx documentation generator. 164 The reference manual is available from an online platform 165 together with tutorials, whitepapers, webinars, etc. Through a search engine it is thus easy to find all available information for a given problem.

XIV. QUANTUMATK APPLICATIONS
A. Large-scale Simulations of 2D Field-effect Transistors As already described in Section XI, the combination of DFT-LCAO with the NEGF method makes it possible to use QuantumATK to simulate the electronic structure and electrical characteristics of devices at the atomistic level. Field-effect transistor (FET) device configurations 166,167 are simulated by including dielectric regions and electrostatic gates, see Section IV E.
Here, we show how this framework can be used to study the electrical characteristics of a tunnel FET (TFET) device, where the channel is formed by a heterojunction based on two-dimensional semiconductors. 168,169 We demonstrate how the characteristics of the device can be tuned by using an asymmetric contact scheme. The latter is similar to that proposed for graphene-based photodetectors, 170 where two different metals are used to contact the graphene channel. Figure 14 shows the 2D-TFET device considered here. The device comprises a semiconducting channel formed by a MoTe 2 /SnS 2 heterojunction. 171   Understanding these trends requires considering that the asymmetric contact scheme has a two-fold effect on the electronic structure of the device. On the one hand, the use of two metals with different work functions leads to an additional built-in electric field in the channel region, when the chemical potentials of the drain and source electrodes, µ D and µ S , are aligned on a common energy scale. On the other hand, the interaction between the metallic contact and MoTe 2 is expected to depend also on the chemical nature of the metal.
The presence of an additional built-in electric field, and its effect on the device electrostatics, are evident by comparing the Hartree difference potential (∆V H ) in the two devices at V GS = 0 V along the channel, as shown in Fig. 16(a,b). While in the SC device the potential changes smoothly along the channel region, a sudden increase in the potential is observed in the ASC device around the E(MoTe 2 ) region. Here, the potential lines run parallel to the transport direction, indicating the presence of a left-pointing local electric field. The sign of this field is consistent with that generated by an asymmetric contact scheme with Φ MS > Φ MD , that is, the same as that of the ASC device.
The projected local density of states (PLDOS) along the devices reveal that the different electrostatics also affect their electronic structure. For both contact schemes, the DOS within the bias window, [µ D − µ S ] ± k B T = ∆µ ± k B T , is strongly inhomogeneous along the channel, as the conduction bands (CBs) of MoTe 2 and SnS 2 are pinned to µ D and µ S , respectively (see Fig. 16(c,d)). This results in a vanishing DOS in the E(MoTe 2 ) region within the bias window. Here, the DOS is even smaller in the ASC device, due to (i) the weaker pinning of the CBs to µ D , and (ii) the effect of the local electric field, which bends and depletes even more the CBs, moving Temperature-dependent phonon-limited resistivity of the three metals Au, Ag and Cu evaluated from firstprinciples simulations using the ATK-LCAO engine. them further away from the bias window. In the SC device, the field is much weaker, and the CBs are bent only in the proximity of the overlap region.
The transconductance behavior can be understood from the combined analysis of ∆V H and of the PLDOS. The DOS within ∆µ ± k B T in the E(MoTe 2 ) region, described in terms of an effective barrier φ MoTe2 , ultimately determines the reverse-bias current in the channel. In the SC device, φ MoTe2 is lower for the case V GS = 0 V, and depends only weakly on V GS , as shown in Fig. 16(e). This results in a higher absolute value of I DS , and in a lower variation of I DS with V GS . Conversely, in the ASC device, φ MoTe2 is higher at comparable values of V GS , and varies appreciably when V GS is increased, see Fig. 16(f). This explains the lower values of the drain-source current, and its higher variation with the gate-source voltage. These trends are consistent with those of the transconductance curves shown in Fig. 15.
In summary, DFT-NEGF simulations for M D /MoTe 2 /SnS 2 /M S ultra-scaled 2D-TFET devices show that the transconductance can be engineered by an appropriate choice of the metallic electrodes, and highlight the importance of atomistic device simulations for optimization of the electrical characteristics of devices based on non-conventional semiconductors.

B. Phonon-limited Mobility of Metals
The continued downscaling of nanoelectronics makes the metal interconnects an increasingly critical part of transistor designs. 175 Present-day transistors use Cu as an interconnect material, and a good understanding of the origin of resistance increase with downscaling of interconnects will be important for the design and performance of future nanoscale devices.
We here present first-principles calculations of the phonon-limited resistivity of three FCC metals; Cu, Ag, and Au. We solve the Boltzmann transport equation for the mobility, using first-principles EPC constants, as described in section VIII E. Such DFT calculation of the resistivity of metals is computationally demanding, as one needs to integrate the EPC over both electron and phonon wave vectors (k-and q-space), and we know of only few studies of the EPC in metals that includes a full integration. 146,176,177 We here show that the tetrahedron integration method enables computationally efficient mobility calculations. The method may therefore be used for computational screening of materials, and first-principles simulations become accessible for identifying promising replacement materials for future interconnects.
To calculate the scattering rate related to EPC, the phonon modes and derivatives of the Hamiltonian with displacements are needed. The supercell method for calculation of phonons and EPC from first principles was described in section VIII, and Fig. 7 showed the phonon band structures of Cu, Ag and Au, calculated using the ATK-LCAO simulation engine. For the integration of the scattering rate in (32) we use a sampling of 20 × 20 × 20 q-points and tetrahedron integration. In addition, we apply the two-step procedure, where a k-space isotropic but energy-dependent scattering rate is used to efficiently evaluate the resistivity. Figure 17 shows the DFT results for the temperaturedependent phonon-limited resistivity of bulk Cu, Ag, and Au (Debye temperatures of 347, 227 and 162 K, 178 respectively). The resistivity increases with temperature as the phonon occupation increases, and becomes linearly dependent on temperature above the Debye temperature. Table X presents the calculated room-temperature bulk resistivities, and compares them to experiments and to calculated values for metal nanowires (NWs) with diameters d = 1 nm. In agreement with experiments, we find that Au has the largest resistivity, and that Ag is more conductive than Cu. In addition, the resistivity increases significantly when forming nanowires of the elements. Despite the fact that the phonon dispersions of bulk Au and Ag are very similar, the resistivity is quite different. In the minimal free-electron model of metals, the conductivity is given by 1/ρ(T ) = 1 3 e 2 v 2 F τ (T )n(ε F ). In the three FCC metals considered here, the Fermi velocity, v F , and the DOS, n(ε F ), (and resulting carrier density) are almost identical, and the difference in the resistivity is traced back to the variation in the scattering rate. This shows how full first-principles Boltzmann transport simulations of the scattering rate is needed to capture the origin of the resistivities of different metals. While the resistivity of bulk Ag is slightly underestimated by the simulations, we find good agreement with experiments for bulk Au and Cu, as well as the correct ranking of the individual metals. This illustrates the predictive power of the method. In general, we find that the resistivity of d = 1 nm nanowires is increased by a factor of three for Au and even more for Ag and Cu, as compared to bulk, due to the increased EPC in nanowires.

C. Multi-model Dynamics with an Applied Electric Field
The tight integration of different atomic-scale simulation engines within the same software framework allows for straight-forward combination of multiple atomistic models into one single simulation workflow. This enables elaborate computational workflows and extend the functionality of QuantumATK beyond that of methods based on a single atomistic model. We here show how such a multi-model approach can be used to implement a hybrid method that combines classical FF MD simulations with a DFT description of time-dependent fluctuations of the atomic charges as the MD simulation progresses.
We study here LiFePO 4 , a promising cathode material of the olivine family for Li-ion batteries. 141,179 In this class of materials, the olivine scaffold provides natural diffusion channels for the Li + ions, which have been shown to diffuse via a hopping mechanism, preferentially through [010]-oriented channels. 180,181 MD simulations aimed at understanding the diffusion process have focused mainly on its temperature (T ) dependence. In this case, relatively high temperatures, usually in the range 500 to 2000 K, are required to reach a sufficiently high hopping probability within a reasonable MD simulation time, and allow for calculation of the associated diffusion constants. These simulations have demonstrated that the diffusion increases with T , as a natural consequence of the increased hopping probability favored by Brownian motion.
However, in an electrochemical cell under operating conditions, the motion of the Li + ions may also have a non-negligible drift component, due to the displacement field resulting from the voltage difference applied between the anode and the cathode. This potentially rather important effect is rarely taken into account in atomistic simulations. 182,183 Another significant issue in the simulation of Li-ion batteries is related to the inclusion of electronic effects. In order to reach reasonably long simulation times, to describe atom diffusion at temperatures close to 300 K, most low-T MD simulations are based on FFs, which by construction neglect any time-dependent fluctuations of the electronic density during the MD run. A number of models have tried to address this issue by either including approximate models to account for the charge fluctuation, 184 or by running semi-classical dynamics on precalculated potential-energy surfaces based on DFT. 185 A QuantumATK multi-model approach can be used to address these issues by including first-principles charge fluctuations in the FF MD. The applied displacement field should add a force term F a = Q a D on the a-th ion with formal charge Q a and D being the field vector. However, in a FF, Q a is a time-independent parameter, so the field-induced force will also be time-independent. In a multi-model approach, we instead use DFT simulations to determine the instantaneous charge Q a at regular intervals during the MD. Time-dependency in the field-induced force term is then included by use of a MD hook function (see Section VII D) by defining the timedependent formal charge Q a (t) as where ∆Q a describes the time-dependent fluctuation. In principle, ∆Q a can be defined arbitrarily, provided that charge neutrality is maintained in the system. In the present case, we chose a simple definition, where Q DFT a (t) and Q DFT,ref a are the time-dependent charge of the i-th atom obtained from a DFT calculation for the MD configuration at time t, and a timeindependent charge obtained for a reference configuration at T = 0 K, respectively. We note that, in the present case, the lack of consistency between the methods used to calculate Q FF a and ∆Q a (t) does not constitute an issue, since the charge fluctuations during the dynamics are of the order ∆Q a ∼ 0.1 e − .
We have applied this multi-model approach to investigate the interplay between Brownian and drift components of the diffusion of Li + ions along the [010] channels in LiFePO 4 in the presence of an applied displacement field. The system was described by a 1 × 2 × 1 LiFePO 4 112-atom supercell, that is, 2 times the conventional unit cell (16 formula units). For the classical part of the multi-model simulations, we used a FF potential by Pedone et al., 154 which has been shown to describe qualitatively correctly the geometry and transport properties of olivine materials. 186 The ATK-LCAO engine was used for the DFT part. MD simulations were performed at temperatures 300 K and 1000 K for a displacement field D = [0, D y , 0], with 0.0 V/Å ≤ D y ≤ 0.3 V/Å. For each temperature, a 5 ps equilibration run using a NPT ensemble was performed, starting from the structure optimized at 0 K, using a Maxwell-Boltzmann distribution of initial velocities, followed by a 45 ps production run using a NVT ensemble. The MD time step was 1.0 fs, and ∆Q a (t) was recalculated every 100 MD steps, see (84), with Q DFT a (t) and Q DFT,ref a obtained from Mulliken population analysis. Further computational details are given in Appendix A. Figure 18(a,c) shows the average displacement y Li + of the Li + ions along the y Cartesian direction, that is, along the [010] channels of the FePO 4 scaffold, calculated for temperatures 300 and 1000 K and for increasingly higher values of the applied field, using either FFs only or the FF+DFT multi-model approach. In the absence of an applied field and at 300 K, the average Li + -ion displacement remains constant at y Li + = 4.67 ± 0.11Å during the entire simulation, indicating the absence of hopping events. At 1000 K, the situation is rather similar, as y Li + increases only slightly from an initial value of 7.12 ± 0.19Å (obtained from an average of the snapshots collected during the first picosecond of the FF-only MD) to a final value of 9.16 ± 0.13Å (obtained from an average of the snapshots collected during the last picosecond). For the multi-model simulation, we observe instead a small decrease of y Li + over time. This indicates that, at both temperatures, Li + hopping due to Brownian motion is a rare event.
Applying an increasingly stronger displacement field leads to a progressive increase in the Li + hopping probability. At 300 K, the average Li + -ion displacement increases steadily from the beginning of the MD run for D y ≥ 0.20 V/Å, indicating that, for these values of D y , Li + hopping is primarily due to field-induced drift. The Li + ions accelerate until they reach a constant velocity, as shown by the tendency of the y Li + vs. time curves to continually decrease their slope, corresponding to a straight line on a linear plot.
In the absence of an applied field, increasing the temperature should increase the probability of Li + ion diffusion due to increased Brownian motion. 181,186 However, in the present case we find that the Li + ions move less at 1000 K than at 300 K. For comparable values of D y , the y Li + vs. time curve has a smaller slope at 1000 K than those calculated at 300 K. The reason is that collision events of the Li + ions with the LiFePO 4 lattice, where phonons are considerably more excited at higher temperatures than at room temperature, limits the effective velocity of the Li + ions. This is evident by comparing the LiFePO 4 structures at the two temperatures. Figure 18(c,d) shows two snapshots extracted at the end of the MD runs at D y = 0.3 V/Å and at temperatures 300 and 1000 K, respectively. At 300 K, the LiFePO 4 structure is relatively unperturbed. Consequently, the Li + ions are able to travel through the [010] channels with relatively few scattering events with the LiFePO 4 lattice. Conversely, at 1000 K, the LiFePO 4 structure is significantlty perturbed, leading to a high probability of collisions between the Li + ions and the olivine lattice.
In summary, we have studied the diffusion of Li + in olivine LiFePO 4 , using a multi-model computational approach that combines a classical FF with DFT, the latter to include the effect of the field and of time-dependent charge fluctuations. Our analysis highlights the importance of considering the combined effect of both Brownian and drift contributions to the Li + hopping to describe the overall process, which strongly depends on not only the temperature itself, but also on the probability of collision events between the diffusing ions and the FePO 4 lattice.

0%
20% 40% 60% 80% 100% Ge content in SiGe alloy  19. Conduction band energies (EX and EL) of Si1−xGex alloy as a function of Ge content, x, calculated using PPS-PBE and HSE06 functionals in combination with the LCAO basis set and the PW basis set, respectively. The EX and EL energies are defined with respect to the top of the valence band (E val ) at the Γ-point. Details on the definition of band energies at special k-points for disordered alloys can be found elsewhere. 187 Reference experimental data (open markers) on the band-gap compositional dependence, Egap(x), are given for low (4.2 K) and room (296 K) temperatures. 188 The dashed (solid) lines correspond to linear (quadratic) interpolation of the DFT-calculated band energies, EL (EX), given with filled markers; the interpolation formulas are given in Table XI.

D. Electronic Structure of Binary Alloys
Understanding the physical properties of semiconductor alloys, such as silicon-germanium binary compounds, is highly relevant, since such alloys are commonly used in microelectronics as a semiconductor material for, e.g., heterojunction bipolar transistors or as a strained semiconductor layer in CMOS transistors. 189 Moreover, device-level TCAD simulations, frequently used in industrial semiconductor research and development, usually require material-dependent input parameters such as band gap, effective masses, deformation potentials, and many others. 190 Atomic-scale simulations may be used to calculate such parameters from first principles if ex-perimental values are not available, including composition dependence. 187 However, simulating randomly disordered alloys may be computationally challenging since the traditional approach to random-alloy (RA) simulations use stastical sampling of multiple large supercells with random atomic arrangements (configurational averaging) to take into account the effect of disorder on the physical properties of alloys.
We here adopt the special quasi-random structure (SQS) approach 191 for DFT modelling of SiGe random alloys, which significantly reduces the computational cost. Unlike in the RA approach, in the SQS method the configurational averaging of band energies is captured by a single supercell structure. We study 64-atom Si 1−x Ge x supercells in the full range of compositions, 0 ≤ x ≤ 1, by calculating composition-dependent lattice constants and band energies. The PPS-PBE method, 19 discussed in Section IV D, was used with the ATK-LCAO simulation engine, and we compare the band energies to those obtained with the HSE06 hybrid functional using the ATK-PlaneWave engine. We also compare SQS band energies to those calculated using the traditional RA approach, obtained by averaging over 5 randomly generated RA configurations. In both the SQS and RA cases, the band energies were computed by averaging the energies of the conduction (valence) band states split by alloy disorder. 187 We used a NanoLab SQS module to generate the SQS configurations. The module uses a genetic algorithm to optimize finite-size alloy configurations to reproduce selected correlation functions of the infinite alloy system. The genetic optimization algorithm is very efficient, and systems with many hundred atoms are easily handled. In this case, the SQS structures were generated by fitting all pair, triplet, and quadruplet correlation functions with figure sizes up to 7.0, 5.0 and 4.0Å, respectively, such as to match those correlation functions for the truly random alloy, as detailed in Refs. 191 and 192. Generation of a single 64-atom SiGe SQS alloy takes about 4 minutes on a modern 4-core processor. The alloy configurations were then relaxed using PPS-PBE followed by band structure analysis. HSE06-level band structures were calculated without further relaxation. More computational details are given in Appendix A. Figure 19 shows the Si 1−x Ge x composition dependent conduction band minima (CBM), referenced to the valence band maximum, for both the X-and L-valley in the SiGe BZ. We first note that SQS band energies are very similar to the those calculated using the more expensive RA approach. It is well known that the Si 1−x Ge x fundamental band gap changes character at x ∼ 0.85. The PPS-PBE and HSE06 predictions of the transition point are x ∼ 0.88 and 0.82, respectively. As expected, the calculated X-valley conduction-band energies exhibit bowing, i.e., nonlinear behavior of these quantities with respect to Ge content, x. The best-fit interpolation formulas, shown as lines in Fig. 19, are listed in Table XI ters. The PPS-PBE band gaps are in good agreement with room-temperature experiments (within ∼50 meV for the entire range of Ge content), while the HSE06 band gaps are in better agreement with low-temperature experiments. Moreover, the HSE06-based approach appears to more accurately describe the band-gap bowing parameter, while PPS-PBE tends to overestimate it. Finally, the calculated SiGe lattice constant also exhibits compositional bowing, as indicated by the interpolation formula in Table XI. The bowing parameter of 0.034Å is overestimated by ∼26% as compared to experiments (0.027Å).
To benchmark the empirical PPS-PBE method against the parameter-free HSE06 approach, we also calculated the band structure of bulk Si and Ge using both methods, as shown in Fig. 20. The PPS-PBE conduction and valence bands around the Fermi energy are in good agreement with the HSE06 band structure. This is consistent with the fact that the PPS-PBE method was fitted to ex-perimental data, and that the HSE06 hybrid functional accurately simulates the band structure of bulk semiconductors.
In summary, we find that the SQS approach is well suited to describe the compositional bowing of the band energies in Si 1−x Ge x random alloys, suggesting that SQS provides an accurate and efficient approach to randomalloy simulations. The HSE06 hybrid functional accurately describes the conduction-band energies of SiGe alloys and their compositional bowing, while the PPS-PBE method offers a computationally efficient alternative if only bands around the Fermi level are important.

XV. SUMMARY
In this paper we have presented the QuantumATK platform and details of its atomic-scale simulation engines, which are ATK-LCAO, ATK-PlaneWave, ATK-SE, and ATK-ForceField. We have compared the accuracy and performance of the different engines, and illustrated the application range of each. The platform includes a wide range of modules for application of the different simulation engines in solid-state and device physics, including electron transport, phonon scattering, photocurrent, phonon-limited mobility, optical properties, static polarizations, molecular dynamics, etc.
The simulation engines are complimentary and through the seamless Python integration in the Quantu-mATK platform, it is easy to shift between different levels of theory or integrate different engines into complex computational workflows. This has been illustrated in several application examples, where we for example showed how ATK-LCAO and ATK-ForceField can be combined to study Li + -ion drift in a battery cathode material. We also presented applications of QuantumATK for simulat-ing electron transport in 2D materials, phonon-limited resistivity of metals, and electronic-structure simulations of SiGe random alloys.
While several of the simulation engines and methods have been described independently before, 19,24,28,106,123,132,137,147,153,163 we have here provided an overview of the entire platform, including implementation details not previously published. We expect that this paper can become a general reference for documenting the QuantumATK platform, and is a reference to its applications for atomic-scale modelling in semiconductor physics, catalysis, polymer materials, battery materials, and other fields.

Appendix A: Computational Details
In the simulations presented in Fig. 1, we have considered noncrystalline a-Al 2 O 3 structures with a constant density of 2.81 g/cm 3 . The system sizes considered were formed by 5,30,60,120,240,480,960,1920,3840,7680,15360, and 30720 atoms, respectively, each system with the appropriate unit-cell volume. The amorphous phases were generated by randomizing the structure at 5000 K and then quenching to 0 K to avoid extremely small bond distances. The MD simulations were then performed at 300 K using a random Boltzmann distribution of initial velocities and a Langevin thermostat. The FF simulations used a Pedone potential, 154 while the TB simulations used a Slater-Koster parametrization. For the DFT-LCAO and DFT-PW simulations, we used normconserving PseudoDojo pseudopotentials with a Medium basis set and a kinetic-energy cutoff energy of 1360 eV (50 Ha), respectively. For TB, DFT-LCAO, and DFT-PW simulations, the Brillouin zone was sampled using a Monkhorst-Pack 193 (MP) k-point density of 3-4Å. For systems with sizes between 240 and 960 atoms, 2 processes/k-point was used, whereas for the 1920-atom system, 16 processes/k-point was used.
For the DFT-NEGF device simulations presented in Section XIV A, we used the PBE density functional with SG15-Medium (FHI-DZP) combinations of pseudopotentials and basis sets for MoTe 2 and SnS 2 (Au and Al). The real-space cutoff energy was 2721 eV (100 Ha), and MP k-point grids of 12×1×100 and 12×1 were used to sample the BZ of the electrode and of the device, respectively.
In the study of multi-model dynamics presented in Section. XIV C, we used the ATK-LCAO engine with a DZP basis set and a real-space cutoff energy of 2180 eV (80 Ha). Exchange-correlation effects were described by the PBE functional, and the FePO 4 BZ was sampled using a 3 × 3 × 2 MP k-point grid.
For the electronic-structure calculations for SiGe random alloys presented in Section XIV D, we used a 3×3×3 MP k-point grid and an electron temperature of 0.025 eV for the Fermi-Dirac occupation function. SG15 (FHI) pseudopotentials were used for the PSS-PBE (HSE06) simulations. The LCAO mesh density cutoff was 2721 eV (100 Ha), and the PW kinetic-energy cutoff was 544 eV (20 Ha). The LCAO simulations used Medium (High) bais sets for silicon (germanium). Relaxation of unit-cell volume and ion positions was done using the PPS-PBE method with total energy, forces and stress converged to 10 −5 eV, 0.01 eV/Å, and 0.05 GPa, respectively.