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

TL;DR

QuantumATK is a comprehensive atomic-scale modeling platform integrating various simulation methods like DFT, tight-binding, and force fields, enabling complex multi-physics simulations for materials and electronic devices.

Contribution

This paper provides a general overview of QuantumATK, highlighting its integrated modules, implementation details, and application examples, showcasing its versatility and capabilities.

Findings

Demonstrated phonon-limited mobility calculations for Cu, Ag, Au

Simulated electron transport in a gated 2D device

Modeled lithium ion drift in battery cathodes

Abstract

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.

Tables

Table 1. Table 1: Simulation engines in the QuantumATK platform, with examples of other simulation platforms using the same underlying methodology. LCAO and PW means linear combination of atomic orbitals and plane wave, respectively.

Engine	Description	First release	Related platforms
ATK-LCAO	Pseudopotential DFT using LCAO basisSmidstrup et al. (2017)	2003	SIESTA,Soler et al. (2002) OpenMXOzaki (2003)
ATK-PlaneWave	Pseudopotential DFT using PW basis	2016	VASP,Kresse and Hafner (1993) Quantum ESPRESSOGiannozzi et al. (2009)
ATK-SE	Semi-empirical TB methodsStokbro et al. (2010)	2010	DFTB+,Aradi et al. (2007) NEMO,Klimeck et al. (2002) OMENKlimeck and Luisier (2010)
ATK-ForceField	All types of empirical force fieldsSchneider et al. (2017)	2014	LAMMPS,LAM GULPGale and Rohl (2003)

Table 2. Table 2: Summary of QuantumATK  Δ Δ \Delta -tests for elemental solids, Lejaeghere  et al. ( 2016 ) and RMS errors of the lattice constant ( a 𝑎 a ) and bulk modulus ( B 𝐵 B ) of rock-salt and perovskite test sets, Garrity  et al. ( 2014 ) using SG15 and PseudoDojo PPs, and the ATK-PlaneWave and ATK-LCAO engines, the latter with different basis sets. All errors are relative to all-electron calculations.

	Medium	High	Ultra	PW
Elemental solids: Delta tests
SG15 $\Delta$ (meV)	3.45	1.88	2.03	1.32
PseudoDojo $\Delta$ (meV)	4.53	1.52	1.40	1.04
Rock salts: RMS of $a$ and $B$
SG15 (%)	0.40	0.24	0.23	0.16
PseudoDojo (%)	0.50	0.18	0.15	0.09
Perovskites: RMS of $a$ and $B$
SG15 (%)	0.36	0.24	0.18	0.13
PseudoDojo (%)	0.35	0.21	0.13	0.06

Table 3. Table 3: Fundamental band gaps (in units of eV) for a range of semiconductors and simple oxides, calculated using different XC methods, and compared to experimental values. The ATK-LCAO simulation engine was used for PBE, TB09, and PBE-1/2 calculations, while the ATK-PlaneWave engine was used for simulations using HSE06. PseudoDojo PPs were used, combined with Ultra basis sets for DFT-LCAO, except for TB09 calculations, which were done using FHI-DZP. Default cutoff energies were used, and a 𝐤 𝐤 \mathbf{k} -point grid density of 7 Å. For bulk silicon, this corresponds to a 15 × 15 × 15 15 15 15 15\times 15\times 15 𝐤 𝐤 \mathbf{k} -point grid. Experimental band gaps are from Ref.  Heyd  et al. , 2005 unless otherwise noted. The bottom row lists the RMS deviation between theory and experiments.

Material	Experiment	PBE	TB09	PBE-1/2	HSE06
C	$5.48$	$4.19$	$5.11$	$5.59$	$5.33$
Si	$1.17$	$0.57$	$1.20$	$1.16$	$1.17$
Ge	$0.74$	$0.00$	$1.11$	$0.81$	$0.55$111Direct band gap ($\Gamma \to \Gamma$), different in size from the 0.72 eV reported in Ref. Schimka et al., 2011, but similar to the 0.56 eV reported in Ref. Heyd et al., 2005, both using theoretical lattice constants rather than experimental ones.
SiC	$2.42$	$1.36$	$2.31$	$2.66$	$2.27$
BP	$2.40$	$1.24$	$1.79$	$1.63$	$2.01$
BAs	$1.46$	$1.25$	$1.94$	$1.58$	$2.05$
AlN	$6.13$	$4.16$	$6.97$	$5.83$	$5.54$
AlP	$2.51$	$1.55$	$2.36$	$2.46$	$2.30$
AlAs	$2.23$	$1.45$	$2.45$	$2.38$	$2.27$
AlSb	$1.68$	$1.22$	$1.82$	$1.92$	$1.76$
GaN	$3.50$	$1.89$	$4.10$	$3.27$	$2.87$
GaP	$2.35$	$1.59$	$2.38$	$2.22$	$2.26$
GaAs	$1.52$	$0.63$	$1.81$	$1.23$	$1.11$
GaSb	$0.73$	$0.11$	$0.76$	$0.52$	$0.64$
InN	$0.69$	$0.00$	$1.74$	$1.20$	$0.49$
InP	$1.42$	$0.69$	$2.17$	$1.30$	$1.26$
InAs	$0.41$	$0.00$	$1.08$	$0.51$	$0.23$
InSb	$0.23$	$0.00$	$0.49$	$0.32$	$0.27$
TiO2	$3.0$222Ref. Landmann et al., 2012.	$1.91$	$3.11$	$3.00$	$3.37$
SiO2	$8.9$333Ref. Bersch et al., 2008.	$6.07$	$11.31$	$8.16$	$7.83$
ZrO2	$5.5$33footnotemark: 3	$3.65$	$4.96$	$5.26$	$5.16$
HfO2	$5.7$33footnotemark: 3	$4.17$	$5.54$	$5.87$	$5.76$
ZnO	$3.44$444Ref. Berger, 2017.	$0.95$	$3.24$	$2.78$	$2.47$
MgO	$7.22$	$4.79$	$8.51$	$6.75$	$6.49$
RMS error		$1.34$	$0.71$	$0.33$	$0.43$

Table 4. Table 4: Silicon and germanium equilibrium lattice constants and fundamental band gaps, both calculated using the PPS-PBE XC method, and compared to experiments at 300 K. The SG15-High combination of PPs and LCAO basis sets was used, and a 15 × 15 × 15 15 15 15 15\times 15\times 15 𝐤 𝐤 \mathbf{k} -point grid. The lattice constants were determined by minimizing the first-principles stress on the primitive unit cells, using a maximum stress criterion of 0.1 GPa (0.6 meV/Å 3 ).

Material	Property	PPS-PBE	Experiment
Silicon	Lattice constant	5.439 Å	5.431 Å
	Band gap	1.14 eV	1.12 eV
Germanium	Lattice constant	5.736 Å	5.658 Å
	Band gap	0.65 eV	0.67 eV

Table 5. Table 5: 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 𝐻 H ) or non-orthogonal ( H , S 𝐻 𝑆 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.	Type	Range	$E,F,\sigma$
Hückel	Ammeter et al.,1978	SK, ($H,S$)	long	no
Empirical TB	Vogl et al.,1983	SK, ($H$)	short	no
DFTB	Elstner et al.,1998	SK, ($H,S$)	medium	yes
Purdue	Boykin et al.,2002	Env, ($H$)	short	no
NRL	Bernstein et al.,2000	Env, ($H,S$)	long	yes

Table 6. Table 6: Selected potential models included in ATK-ForceField.

Potential model	Special properties	References
Stillinger–Weber (SW)	Three-body	Stillinger and Weber,1985
Embedded atom model (EAM)	Many-body	Mishin et al.,2001
Modified embedded atom model (MEAM)	Many-body	Baskes,1997
	Directional bonding
Tersoff/Brenner	Bond-order	Tersoff,1988; Brenner et al.,2002
ReaxFF	Bond-order	Chenoweth et al.,2008
	Dynamical charges
COMB/COMB3	Bond-order	Yu et al.,2007
	Dynamical charges
	Induced dipoles
Core-shell	Dynamical charge fluctuations	Mitchell and Fincham,1993
Tangney–Scandolo (TS)	Induced dipoles	Tangney and Scandolo,2002
Aspherical ion model	Induced dipoles and quadrupoles	Rowley et al.,1998
	Dynamical ion distortion
Biomolecular and valence force fields	Static bonds	Mackerell,2004; Keating,1966

Table 7. Table 7: ATK-ForceField timings as compared to LAMMPS. Plimpton ( 1995 ) Absolute timings for molecular dynamics (MD) simulations, in units of microseconds per atom per MD step, using one computing core for all potentials. Potential abbreviations are defined in Table  6 . In addition, LJ means Lennard–Jones. More details of the benchmark systems can be found in Ref.  Schneider  et al. , 2017 .

	LJ	Tersoff	SW	EAM	ReaxFF	COMB	TS
QuantumATK	3.8	6.3	5.2	3.8	180	320	360
LAMMPS	1.9	7.8	5.2	2.4	190	240	N/A

Table 8. Table 8: Born effective charges ( Z ∗ superscript 𝑍 Z^{*} ) and piezoelectric tensor components ( ϵ 33 subscript italic-ϵ 33 \epsilon_{33} and ϵ 14 subscript italic-ϵ 14 \epsilon_{14} ) for III-V wurtzite nitrides and zincblende GaAs. Reference vales for the nitrides are from Ref.  Bernardini  et al. , 1997 and from Ref.  King-Smith and Vanderbilt, 1993 for GaAs. QuantumATK calculations were performed using the DFT-LCAO engine with the LDA XC functional and a DZP basis set.

	$Z^{\ast }$		${ϵ}_{33}$
	Reference	QuantumATK	Reference	QuantumATK
AlN	$-$2.70	$-$2.67	1.46	1.65
GaN	$-$2.72	$-$2.75	0.73	0.86
InN	$-$3.02	$-$2.98	0.97	1.21
			${ϵ}_{14}$
GaAs	$-$1.98	$-$2.07	$-$0.28	$-$0.26

Table 9. Table 9: MAE (in units of meV) for various Fe-based L 1 0 subscript 1 0 1_{0} phases. Atomic structures and reference results (SIESTA and VASP) are from Ref.  Blanco-Rey  et al. , 2018 . The QuantumATK selfconsistent and non-selfconsistent calculations were performed with a 17 × × 17 × 14 17\times\times 17\times 14 𝐤 𝐤 \mathbf{k} -point grid, while the band energies were sampled on a 40 × 40 × 34 40 40 34 40\times 40\times 34 𝐤 𝐤 \mathbf{k} -point grid. PseudoDojo pseudopotentials were used for both LCAO and PW calculations. The High basis set was used for LCAO.

Structure	SIESTAa	VASPa	QuantumATK	QuantumATK
	LCAO	PW	LCAO	PW
FeCo	0.45	0.55	0.66	0.66
FeCu	0.42	0.45	0.45	0.45
FePd	0.20	0.13	0.12	0.15
FePt	2.93	2.78	2.43	2.57
FeAu	0.36	0.62	0.22	0.56

Table 10. Table 10: First-principles phonon-limited resistivities at 300 K (in units of n ​ Ω ⋅ m ⋅ n Ω m \textrm{n}\Omega\cdot\textrm{m} ), compared with experimental values from Ref.  Kasap  et al. , 2017 . Au nanowire results from Ref.  Markussen  et al. , 2017 .

	DFT	Experiment
Au, Bulk	15.9	20.5
Au, NW ($d\approx 1$ nm)	56.0	-
Ag, Bulk	4.9	14.7
Ag, NW ($d\approx 1$ nm)	28.7	-
Cu, Bulk	14.2	15.4
Cu, NW ($d\approx 1$ nm)	98.3	-

Table 11. Table 11: First-principles interpolation formulas for the Si 1-x Ge x composition-dependent band gap and lattice constant. The variables b w subscript 𝑏 w b_{\text{w}} , b w ′ superscript subscript 𝑏 w ′ b_{\text{w}}^{\prime} , and b w ′′ superscript subscript 𝑏 w ′′ b_{\text{w}}^{\prime\prime} are bowing parameters.

	Band gap interpolation formula (eV)
PPS-PBE	$E_{\text{X}}=1.116-0.764x+b_{\text{w}}^{\prime }{x}^2$
	$b_{\text{w}}^{\prime }=0.526$ eV
	$E_{\text{L}}=2.104-1.425x$
HSE06	$E_{\text{X}}=1.204-0.444x+b_{\text{w}}^{\prime \prime }{x}^2$
	$b_{\text{w}}^{\prime \prime }=0.228$ eV
	$E_{\text{L}}=2.032-1.267x$
	Lattice constant interpolation formula (Å)
PPS-PBE	$a\left(x\right)=5.431+0.257x+b_{\text{w}}{x}^2$
	$b_{\text{w}}=0.034$ Å

Equations

$$V^{\text{eff}}[n]=V^{\text{H}}[n]+V^{\text{xc}}[n]+V^{\text{ext}}[n],$$

$$V^{\text{ext}}=\sum _a{V}_a^{\text{pseudo}}+V^{\text{gate}},$$

$${\hat{H}}^{\text{KS}}=-\frac{{\hslash }^2}{2m}{\nabla }^2+V^{\text{eff}},$$

$${\hat{H}}^{\text{KS}}{\psi }_{\alpha }={ϵ}_{\alpha }{\psi }_{\alpha }.$$

$$E[n]=T+E^{\text{H}}[n]+E^{\text{xc}}[n]+E^{\text{ext}}[n],$$

$${\psi }_{\alpha }(\mathbf{r})=\sum _i{c}_{\alpha i}{\varphi }_i(\mathbf{r}).$$

$$\sum _j{H}_{ij}^{\text{KS}}{c}_{\alpha j}={\epsilon }_{\alpha }\sum _j{S}_{ij}{c}_{\alpha j},$$

$$D_{ij}=\sum _{\alpha }{c}_{\alpha i}^{\ast }{c}_{\alpha j}f\left(\frac{{\epsilon }_{\alpha }-{\epsilon }_{\text{F}}}{k_{\text{B}}T}\right),$$

$$n(\mathbf{r})=\sum _{ij}{D}_{ij}{\varphi }_i(\mathbf{r}){\varphi }_j(\mathbf{r}),$$

$${\psi }_{\alpha }(\mathbf{r})=\sum _{|\mathbf{g}|<g_{\text{max}}}{c}_{\alpha ,\mathbf{g}}{\mathrm{e}}^{\mathrm{i}\mathbf{g}\cdot \mathbf{r}},$$

$$\frac{{\hslash }^2{g}_{\text{max}}^2}{2m}<E_{\text{cut}}.$$

$${\varphi }_{nlm}(\mathbf{r})={\chi }_{nl}(r)Y_{lm}(\hat{\mathbf{r}}),$$

$$E_U=\frac{1}{2}\sum _l{U}_l(n_l-n_l^2),$$

$${\nabla }^2{V}^{\text{H}}(\mathbf{r})=-\frac{e^2}{4\pi {ϵ}_0}n(\mathbf{r}),$$

$$V^{\text{H}}(\mathbf{r})=V_0,\mathbf{r}\in \Omega .$$

$${\nabla }^2{V}^{\text{H}}(\mathbf{r})$$

$$V^{\text{H}}(\mathbf{r})=V_0(\mathbf{r}),\mathbf{r}\in S.$$

$$\frac{\partial V^{\text{H}}(\mathbf{r})}{\partial \mathbf{n}}=\mathbf{n}\cdot \nabla V^{\text{H}}(\mathbf{r})=V_0^{\prime }(\mathbf{r}),\mathbf{r}\in S,$$

$$U_l=\frac{2e^2}{\sqrt{\pi }{\sigma }_l}.$$

$$U_l=\frac{d{\epsilon }_l}{d{n}_l},$$

$$U=\frac{1}{N}\sum _l{n}_l{U}_l,$$

$$H_{ij}=H_{ij}^0+\frac{1}{2}(V({\mathbf{r}}_i)+V({\mathbf{r}}_j))S_{ij},$$

$$H_{ij}^{\sigma }=\pm \frac{1}{2}{S}_{ij}\left(d{E}_{l_i}+d{E}_{l_j}\right),$$

$$d{E}_l=\sum _{l^{\prime }\in {\mu }_l}{W}_{l{l}^{\prime }}(m_{l^{\prime }\uparrow }-m_{l^{\prime }\downarrow }).$$

$$k_{\mathrm{HTST}}=\frac{{\prod }_i^{3N}{\nu }_i^{\text{min}}}{{\prod }_i^{3N-1}{\nu }_i^{‡}}\exp \left[-\left(E^{‡}-E^{\text{min}}\right)/k_{\text{B}}T\right],$$

$$D_{a,\alpha ;b,\beta }=\frac{1}{\sqrt{m_a{m}_b}}\frac{d{F}_{b,\beta }}{d{r}_{a,\alpha }},$$

$$g(\omega )=\frac{1}{N}\sum _{\mathbf{q}\lambda }\delta (\omega -{\omega }_{\lambda \mathbf{q}}),$$

$$(\delta {\hat{H}}_{{\mathbf{r}}_a}{)}_{ij}=⟨{\varphi }_i|\frac{\partial \hat{H}}{\partial {\mathbf{r}}_a}|{\varphi }_j⟩,$$

$$g_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}^{\lambda n{n}^{\prime }}=⟨n^{\prime }{\mathbf{k}}^{\prime }|\delta {\hat{H}}_{\lambda \mathbf{q}}|n\mathbf{k}⟩,$$

$$j_{\alpha }={\sigma }_{\alpha \beta }{\mathcal{E}}_{\beta }+{\sigma }_{\alpha \beta \gamma }{\mathcal{E}}_{\beta }{B}_{\gamma }+{\nu }_{\alpha \beta }{\nabla }_{\beta }T,$$

Full text

1. I Introduction

2. II Overview

3. III Atomistic Configurations

4. IV DFT Simulation Engines

5. IV.1 LCAO Representation

6. IV.2 PW Representation

7. IV.3 Pseudopotentials and LCAO Basis Sets

8. IV.4 Exchange-Correlation Methods

9. IV.4.1 Semilocal functionals

10. IV.4.2 Hybrid functionals

11. IV.4.3 Semiempirical methods

12. IV.4.4 DFT+U methods

13. IV.5 Boundary Conditions and Poisson Solvers

14. IV.5.1 Boundary conditions

15. IV.5.2 Poisson solvers

16. V Semi-Empirical Models

17. VI Empirical Force Fields

18. VII Ion Dynamics

19. VII.1 Local Structural Optimization

20. VII.2 Global Structural Optimization

21. VII.3 Reaction Pathways and Transition States

22. VII.4 Molecular Dynamics

23. VII.5 Adaptive Kinetic Monte Carlo

24. VIII Phonons

25. VIII.1 Calculating the Dynamical Matrix

26. VIII.2 Wigner–Seitz Method

27. VIII.3 Phonon Band Structure and Density of States

28. VIII.4 Electron-Phonon Coupling

29. VIII.5 Transport Coefficients

30. IX Polarization and Berry Phase

31. X Magnetic Anisotropy Energy

32. XI Quantum Transport

33. XI.1 NEGF Method

34. XI.2 Retarded Green’s Function

35. XI.3 Complex Contour Integration

36. XI.4 Bound States

37. XI.5 Spill-in Terms

38. XI.6 Device Total Energy and Forces

39. XI.7 Transmission Coefficient and Current

40. XI.8 Inelastic Transmission and Inelastic Current

41. XI.8.1 Special Thermal Displacement Method

42. XI.9 Thermoelectric Transport

43. XI.10 Photocurrent

44. XII QuantumATK Parallelization

45. XII.1 Bulk DFT and Semi-Empirical Simulations

46. XII.2 DFT-NEGF Device Simulations

47. XII.3 FF Simulations

48. XIII NanoLab Simulation Environment

49. XIII.1 Python Scripting

50. XIII.2 NanoLab Graphical User Interface

51. XIII.3 Documentation

52. XIV QuantumATK applications

53. XIV.1 Large-scale Simulations of 2D Field-effect Transistors

54. XIV.2 Phonon-limited Mobility of Metals

55. XIV.3 Multi-model Dynamics with an Applied Electric Field

56. XIV.4 Electronic Structure of Binary Alloys

57. XV Summary

58. A Computational Details

I Introduction

Atomic-scale modelling is increasingly important for industrial and academic research and development in a wide range of technology areas, including semiconductors,Shankar et al. (2008); Zographos et al. (2017) batteries,Shi et al. (2015) catalysis,Nørskov et al. (2008) renewable energy,Islam (2010) advanced materials,Saal et al. (2013) next-generation pharmaceuticals,Trau and Battersby (2001) and many others. Surveys indicate that the return on investment of atomic-scale modelling is typically around 5:1.Goldbeck (2012) 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. Nakai et al., 2014; Xiao et al., 2019 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.Goldbeck (2012) 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.Greeley et al. (2006); Saal et al. (2013); Armiento et al. (2011)

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.Bartlett and Musiał (2007) 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 pre-adjusted 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)Hohenberg and Kohn (1964); Kohn and Sham (1965); Parr and Yang (1994); Kohn et al. (1996) 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.Vogl et al. (1983) 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 electronic-structure 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 1, 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 Al2O3. 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 listed in Table 1 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 atomic-scale 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, QuantumATK 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 silicon-based electronic logics and memory elements,Thirunavukkarasu et al. (2017); Dong et al. (2018) to solar cells composed of novel materialsCrovetto et al. (2017) and next-generation electronic devices based on spintronic phenomona.Sankaran et al. (2016) 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.Zographos et al. (2017) 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,Gunst et al. (2017a) atomistic description of ferroelectricity driven by edge-absorbed polar molecules in gated graphene,Caridad et al. (2018) 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.Palsgaard et al. (2018a) 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.1), 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 QuantumATK 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 1 are invoked using ATK-Python scripting. More details are given in Section XIII.1. QuantumATK 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 shared-memory 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. QuantumATK 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 (non-periodic) 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 semi-infinite 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.Brandbyge et al. (2002)

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.Smidstrup et al. (2017)

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) formulationHohenberg and Kohn (1964); Kohn and Sham (1965); Parr and Yang (1994); Kohn et al. (1996) within the framework of the linear combination of atomic orbitals (LCAO) and plane-wave (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^{\text{eff}}[n]$,

[TABLE]

where $V^{\text{H}}$ is the Hartree potential describing the classical electrostatic interaction between the electrons, $V^{\text{xc}}$ is the exchange-correlation (XC) potential, which in practise needs to be approximated, and $V^{\text{ext}}$ is the sum of the electrostatic potential energy of the electrons in the external potential of ions and other electrostatic field sources. The total external potential is in QuantumATK given by

[TABLE]

where $V_{a}^{\text{pseudo}}$ includes the local ($V^{\text{loc}}_{a}$) and nonlocal ($V^{\text{nl}}_{a}$) contributions to the pseudopotential of the $a$-th atom. The term $V^{\text{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,

[TABLE]

and the single-electron energies ($\epsilon_{\alpha}$) and wave functions ($\psi_{\alpha}$) are solutions to eigenvalue problem

[TABLE]

The electronic ground state is found by iteratively minimizing the KS total-energy density functional, $E[n]$, with respect to the electron density,

[TABLE]

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.

IV.1 LCAO Representation

The DFT-LCAO method uses a LCAO numerical representation of the KS equations, closely resembling the SIESTA formalism.Soler et al. (2002) 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, $\psi_{\alpha}$, are expanded in a set of finite-range atomic-like basis functions $\phi_{i}$,

[TABLE]

The KS equation can then be represented as a matrix equation for determining the expansion coefficients $c_{\alpha i}$,

[TABLE]

where the Hamiltonian matrix $H_{ij}^{\text{KS}}=\langle\phi_{i}|\hat{H}^{\text{KS}}|\phi_{j}\rangle$ and overlap matrix $S_{ij}=\langle\phi_{i}|\phi_{j}\rangle$ 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 multiple-center integrals are computed on a real-space grid.Soler et al. (2002)

For molecules and bulk systems, diagonalization of the Hamiltonian matrix yields the density matrix $D_{ij}$,

[TABLE]

where $f$ is the Fermi–Dirac distribution of electrons over energy states, $\varepsilon_{\text{F}}$ the Fermi energy, $T$ the electron temperature, and $k_{\text{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,

[TABLE]

and is represented on a regular real-space grid, which is the same grid as used for the effective potential in (1).

IV.2 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,

[TABLE]

where $\alpha$ denotes both the wave vector $\mathbf{k}$ and the band index $n$, and ${\bf g}$ are reciprocal lattice vectors. The upper threshold for the reciprocal lattice-vector lengths included in the PW expansion ($g_{\text{max}}$) is determined by a kinetic-energy (wave-function) cutoff energy $E_{\text{cut}}$,

[TABLE]

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_{\text{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 state-of-the-art PW efficiency. For example, we have adopted iterative algorithms for solving the KS equations,Davidson (1975) and fast Fourier transform (FFT) techniques for applying the Hamiltonian operator and evaluating the electron density.Payne et al. (1992); Wende et al. (2016)

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.

IV.3 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.Dal Corso and Conte (2005) Projector augmented-wave (PAW) potentialsBlöchl (1994) 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 SG15Schlipf and Gygi (2015) and PseudoDojoVan Setten et al. (2018) 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 scalar-relativistic formalism.Theurich and Hill (2001); Dal Corso and Conte (2005)

For each PP, we have generated an optimized LCAO basis set, consisting of orbitals $\phi_{nlm}$,

[TABLE]

where $Y_{lm}$ are spherical harmonics, and $\chi_{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.Soler et al. (2002) For the shape of the confinement potential, we follow Ref. Blum et al., 2009.

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.Blum et al. (2009) 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:

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. 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. 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 $\Delta$-testDel ; Lejaeghere et al. (2016) to check the accuracy of the equation of state for elemental, rock-salt, and perovskite solids against all-electron reference calculations, as shown in Table 2. 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 2, the $\Delta$-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 2 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 2 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. Smidstrup et al., 2017.

IV.4 Exchange-Correlation Methods

The XC functional in (5) is the only formal approximation in KS-DFT, since the exact functional is unknown.Kohn and Sham (1965); Parr and Yang (1994); Kohn et al. (1996) 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.Marques et al. (2012) The ATK-PlaneWave engine also allows for calculations using the HSE06 screened hybrid functional.Heyd et al. (2003, 2005); Krukau et al. (2006) 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.Grimme (2006) Both DFT engines support different spin variants for each XC functional: spin-unpolarized and spin-polarized (both collinear and noncollinear). Spin-polarized noncollinear calculations may include spin-orbit interaction through the use of fully relativistic PPs.

IV.4.1 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,Perdew et al. (1996) the PBEsol (designed for solids),Perdew et al. (2008) and the revPBE/RPBE functionals (designed for chemistry applications).Hammer et al. (1999) Recently, the meta-GGA SCAN functionalSun et al. (2015) was also included in QuantumATK, often providing improved accuracy of DFT calculations as compared to PBE.

IV.4.2 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.Heyd et al. (2003); Paier et al. (2006) However, the computational cost is usually much higher than for semilocal approximations. New methodological developments based on the adaptively compressed exchange operator (ACE) methodLin (2016) allow reducing the computational burden of hybrid functionals. The ACE algorithm was recently implemented in QuantumATK for HSE06 calculations, which gives a systematically good description of the band gap of most semiconductors and insulators, see Table 3.

IV.4.3 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,Ferreira et al. (2008, 2011) the TB09 XC potential,Tran and Blaha (2009) and the pseudopotential projector-shift approach of Ref. Smidstrup et al., 2017.

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 3 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 anti-ferromagnetic transition metal oxides.Doumont et al. (2019) 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.Ferreira et al. (2008, 2011)

The Tran–Blaha meta-GGA XC functional (TB09)Tran and Blaha (2009) introduces a parameter, $c$, which can be calculated selfconsistently according to an empirical formula given in Ref. Tran and Blaha, 2009. Table 3 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 $c$-parameters 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 projector-shift (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.Wang and Zunger (1995) The PPS method is usually combined with ordinary PBE calculations.Smidstrup et al. (2017) The two main advantages of this PPS-PBE approach are that (1) for each semiconductor, the projector shifts can be fitted such that the DFT-predicted 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 4 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.

IV.4.4 DFT+U methods

QuantumATK supports the mean-field Hubbard-U correction by Dudarev et al.Dudarev et al. (1998) and Cococcioni et al.,Cococcioni and de Gironcoli (2005) 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,

[TABLE]

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.

IV.5 Boundary Conditions and Poisson Solvers

As already mentioned in Section IV.1, the electron density, $n(\mathbf{r})$ in (9), and the effective potential, $V^{\text{eff}}(\mathbf{r})$ in (3), are in QuantumATK represented on a real-space regular grid. The corresponding Hartree potential $V^{\text{H}}(\bf 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,

[TABLE]

where $e$ is the elementary charge, and $\epsilon_{0}$ is the vacuum permittivity.

In QuantumATK, one may also specify metallic or dielectric continuum regions in combination with a microscopic, atomistic structure, as demonstrated for a 2D device in Fig. 14 in Section XIV.1. This affects the solution of the Poisson equation (14) in the following way. For a metallic region denoted $\Omega$, 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

[TABLE]

For a dielectric region denoted $\Upsilon$, the right-hand side of the Poisson equation will be modified as follows:

[TABLE]

where $\epsilon_{r}$ is the relative dielectric constant, which can be specified as an external parameter in QuantumATK calculations.

IV.5.1 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}({\bf r})$ at the boundary, such that, for a facet $S$ of the simulation cell,

[TABLE]

A Neumann BC means that the normal derivative of the potential on a facet has been fixed to a given function $V_{0}^{\prime}({\bf r})$,

[TABLE]

where ${\bf 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 left- and 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.Neugebauer and Scheffler (1992)

• •

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-hand-side 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.Smidstrup et al. (2017)

• •

Dirichlet BCs for a device configuration. Two-probe 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.

IV.5.2 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 solver,Holst and Saied (1993) a parallel conjugate-gradient-based solver,Concus et al. (1976) and the MUMPS direct solver.Amestoy et al. (2006) 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).Ozaki et al. (2010)

V Semi-Empirical Models

As a computationally fast alternative to DFT, the ATK-SE engine allows for semi-empirical TB-type simulations.Stokbro et al. (2010) The TB models consist of a non-selfconsistent 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.Elstner et al. (1998) The main aspects of these TB models have been described in Ref. Stokbro et al., 2010 and below we give only a brief description of the models.

Table 5 summarizes the available models for the non-selfconsistent 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,Ammeter et al. (1978); Cerda and Soria (2000) Slater–Koster orthogonal TB models,Vogl et al. (1983); Jancu et al. (1998) and DFTB models.Elstner et al. (1998) The ATK-SE engine also supports models that take into account the position of atoms around the two centers. These currently include the environment-dependent TB models from PurdueBoykin et al. (2002) and those from the U.S. Naval Research Laboratory.Bernstein et al. (2000)

It is possible to add a selfconsistent correction to the non-selfconsistent TB models.Stokbro et al. (2010) 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, $\sigma_{l}$, can be related to an onsite repulsion, $U_{l}$, where $l$ is the angular momentum of the orbital. The relation is given byStokbro et al. (2010)

[TABLE]

This onsite repulsion can be calculated from the charge-dependent onsite energies,Elstner et al. (1998)

[TABLE]

where $\varepsilon_{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,Elstner et al. (1998)

[TABLE]

where the average is determined by the number of valence electrons of each orbital, $n_{l}$; $N=\sum_{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(\mathbf{r})$ is calculated using the same methods as used for DFT, see Section IV.5. It is added to the TB Hamiltonian through

[TABLE]

where ${\bf r}_{i}$ is the position of orbital $i$.

The ATK-SE engine also supports spin polarization through the termKöhler et al. (2007)

[TABLE]

where the sign depends on the spin. The spin splitting of shell $l$, $dE_{l_{i}}$, is calculated from the spin-dependent Mulliken populations $m_{l\uparrow},m_{l\downarrow}$ of each shell at the local site ${\mu_{l}}$:

[TABLE]

The shell-dependent spin-splitting strength $W_{ll^{\prime}}$ is calculated from a spin-polarized atomic calculation,Köhler et al. (2007) 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. Schneider et al., 2017, and we therefore only summarize some of the main features.

Table 6 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 sub-systems. 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 7 compares the computational speed of ATK-ForceField molecular dynamics simulations to that of the popular LAMMPS package.Plimpton (1995) 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.3 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. Schneider et al., 2017, so here we only summarize the main features.

VII.1 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. Sheppard et al., 2012, where the changes to the system are described as a combined vector of atomic and strain coordinates.

The default method for optimization is the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) quasi-Newton-type minimization algorithm,Liu and Nocedal (1989) but QuantumATK also implements the fast inertial relaxation engine (FIRE) method.Bitzek et al. (2006)

VII.2 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. QuantumATK 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. Glass et al., 2006. An alternative approach is to perform simulated annealing using molecular dynamics.Kirkpatrick et al. (1983)

VII.3 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.Jónsson et al. (1998) The QuantumATK platform implements state-of-the-art NEB, Henkelman and Jónsson (2000) including the climbing-image method.Henkelman et al. (2000) 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.Smidstrup et al. (2014) 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.Sheppard et al. (2008) In QuantumATK, the NEB forces for each image are combined into a single vector, $\mathbf{F}_{\mathrm{NEB}}\in\mathbb{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.Sheppard et al. (2008)

VII.4 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.Swope et al. (1982) 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,Martyna et al. (1992) an impulsive version of the Langevin thermostat,Goga et al. (2012) or the barostat proposed by Martyna et al. in Ref. Martyna et al., 1994 for isotropic and anisotropic pressure coupling.

Figure 5 also shows that one may apply so-called pre-step hooks and post-step hooks during a QuantumATK MD simulation. These hook functions are scripted in ATK-Python, and may vastly increase flexibility with respect to specialized MD simulation techniques and custom on-the-fly analysis. This makes it easy to employ predefined or user-defined custom operations during the MD simulation. The pre-step hook is called before the force calculation, and may modify atomic positions, cell 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),Müller-Plathe (1997) metadynamics, and other methods. For metadynamics, QuantumATK integrates with the PLUMED package,Tribello et al. (2014) 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.Kondati Natarajan and Behler (2017)

VII.5 Adaptive Kinetic Monte Carlo

Adaptive kinetic Monte Carlo (AKMC) is an algorithm for modelling the long-timescale kinetics of solid-state materials.Henkelman and Jónsson (2001); Xu and Henkelman (2008); Chill and Henkelman (2014) 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 high-temperature 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.Chill and Henkelman (2014); Aristoff et al. (2016)

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),Vineyard (1957)

[TABLE]

where $N$ is the number of atoms, $\nu^{\text{min}}_{i}$ and $\nu^{\ddagger}_{i}$ are the positive (stable) normal-mode frequencies at the minimum and saddle points, $E^{\text{min}}$ and $E^{\ddagger}$ the corresponding energies, $k_{\text{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. Chill and Henkelman, 2014.

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$,

[TABLE]

where $m_{a}$ ($m_{b}$) is the atomic mass of atom $a$ ($b$) and $dF_{b,\beta}/dr_{a,\alpha}$ 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.

VIII.1 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 $\alpha$, and then calculating the resulting forces on atom $b$ along directions $\beta$. 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.Alfè (2009)

VIII.2 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.Parlinski et al. (1997)

VIII.3 Phonon Band Structure and Density of States

The phonon band structure (or phonon dispersion) consists of bands with index $\lambda$ of vibrational frequencies $\omega=\omega_{\lambda\mathbf{q}}$ throughout the Brillouin zone (BZ) of phonon wave vectors $\mathbf{q}$. The phonon density of states (phonon DOS) per unit cell, $g(\omega)$, is defined as

[TABLE]

where $N$ is the number of $\mathbf{q}$-points in the sum. In practice, the phonon DOS is calculated using the tetrahedron method.Blöchl et al. (1994) 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.

VIII.4 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.5) and inelastic scattering of electrons in two-probe devices (see Section XI.8).

To obtain the EPC, we calculate the derivative of the Hamiltonian matrix with respect to the position of atom $a$, $\mathbf{r}_{a}$,

[TABLE]

where $\partial\hat{H}/\partial\mathbf{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, $n\mathbf{k}$, and Fourier transformed using the phonon polarization vectors, to finally obtain the electron-phonon couplings $g$,

[TABLE]

where $\mathbf{q}$ is the phonon momentum and $\lambda$ the phonon branch index.

Further details of how QuantumATK calculates the EPC are given in Ref. Gunst et al., 2016.

VIII.5 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 $\mathbf{j}$ to lowest order in the electric field $\mathcal{E}$, magnetic field $B$, and temperature gradient $\nabla T$,

[TABLE]

where the indices label Cartesian directions and $\sigma_{\alpha\beta}$, $\sigma_{\alpha\beta\gamma}$ and $\nu_{\alpha\beta}$ are the electronic conductivity, Hall conductivity, and thermoelectric response, respectively. Following Ref. Madsen and Singh, 2006, the band-dependent thermoelectric transport coefficients and Hall coefficients are obtained as

[TABLE]

where $\mu$ is the chemical potential and $\epsilon_{\gamma uv}$ the Levi–Civita symbol. The band group velocities $\mathbf{v}(n\mathbf{k})$ and effective mass tensors $\mathbf{M}(n\mathbf{k})$ are obtained from perturbation theory. Importantly, we may in (31) include the full scattering rate $\tau_{n\mathbf{k}}$, and thereby go beyond the constant scattering-rate approximation used in Ref. Madsen and Singh, 2006. As we will see in Section XIV.2, 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

[TABLE]

where $B$ is a temperature-dependent scattering weight,

[TABLE]

where $f$ is the Fermi function, and the scattering angle is defined by

[TABLE]

Furthermore, $P$ ($\bar{P}$) are transition rates due to phonon absorption (emission). They are obtained from Fermi’s golden rule,

[TABLE]

where $n_{\lambda\mathbf{q}}$ is the phonon occupation operator, and $g_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{q}}^{\lambda nn^{\prime}}$ the EPC constant from (29).

QuantumATK offers two different methods for performing the $\mathbf{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 $\mathbf{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 methodBlöchl et al. (1994) 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 $\mathbf{q}$-points increases, using both Gaussian and tetrahedron integration.

The tetrahedron calculation seems converged for $N_{q}=20$, that is, a $20\times 20\times 20$ $\mathbf{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. Samsonidze and Kozinsky, 2018. A two-step procedure is used for the $\mathbf{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 $\mathbf{k}$-space with a low sampling density and a well-converged $\mathbf{q}$-point sampling are used. The initial $\mathbf{k}$-point grid is automatically reduced further by including only $\mathbf{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,

[TABLE]

where we have integrated over bands $n$ and wave vectors $\mathbf{k}$, and $n(E)$ is the density of states. In the second step, we then perform a calculation on a fine $\mathbf{k}$-point grid, but using the energy-dependent isotropic scattering rate $\tau(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 $\mathbf{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.Madsen and Singh (2006) 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. Madsen and Singh, 2006. 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 $\mathbf{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.King-Smith and Vanderbilt (1993)

It is common to divide the polarization into ionic and electronic parts, $\mathbf{P}=\mathbf{P}_{\mathrm{i}}+\mathbf{P}_{\mathrm{e}}$. The ionic part can be treated as a classical electrostatic sum of point charges,

[TABLE]

where $Z_{a}^{\mathrm{ion}}$ and $\mathbf{r}_{a}$ are the valence charge and position vector of atom $a$, $\Omega$ is the unit-cell volume, and the sum runs over all ions in the unit cell.

The electronic contribution to the polarization in direction $\alpha$ is obtained asKing-Smith and Vanderbilt (1993)

[TABLE]

where $\mathbf{R}_{\alpha}$ is the lattice vector in direction $\alpha$, and the Berry phase $\Phi_{\alpha}$ is obtained as

[TABLE]

where the sum runs over $N_{\bot}$ $\mathbf{k}_{\bot}$-points in the BZ plane perpendicular to $\mathbf{R}_{\alpha}$, and

[TABLE]

with the overlap integrals

[TABLE]

and with the $J$ $\mathbf{k}$-points given by $\mathbf{k}_{j}=\mathbf{k}_{\bot}+\mathbf{k}_{\parallel,j}$ lying on a line along the $\mathbf{R}_{\alpha}$ 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\pi$. Consequently, the polarization is a periodic function and constitutes a polarization lattice itself. The polarization lattice in direction $\alpha$ is written as

[TABLE]

where $n$ is an integer labeling a polarization branch, and the polarization quantum in direction $\alpha$ is $\mathbf{P}_{Q,\alpha}=\frac{|e|}{\Omega}\mathbf{R}_{\alpha}$. 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,

[TABLE]

where Voigt notation is used for the strain component, that is, $i\in(xx,yy,zz,yz,xz,xy)$, and the Born effective charge tensor

[TABLE]

where the derivative is with respect to the position of atom $a$ in direction $\beta$.

Table 8 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 ($\parallel$) and out-of-plane ($\perp$) with respect to a crystal plane of atoms, a surface, or an interface between two materials:

[TABLE]

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,Blonski and Hafner (2009); Blanco-Rey et al. (2018) the 2PT method can lead to significantly different results.Blanco-Rey et al. (2018) 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:

Perform a selfconsistent spin-polarized calculation. 2. 2.

For each of the considered spin orientations

1. (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. 2. (b)

Calculate the band energies $\epsilon_{n}$ and projection weights $w_{n,p}$. 3. 3.

Calculate the total MAE as

[TABLE]

where $f_{n}^{\parallel}$ is the occupation factor for band $n$ (including both band and $\mathbf{k}$-point index) for the $\parallel$ spin orientation and $\epsilon_{n}^{\parallel}$ is the corresponding band energy, and likewise for the $\perp$ spin orientation.

The contribution to the total MAE for a particular projection $p$ (atom or orbital projection) is

[TABLE]

where the projection weight is

[TABLE]

with $|\psi_{n}\rangle$ being the eigenstate, $\mathbf{S}$ the overlap matrix, and $\mathbf{P}$ the projection matrix. $\mathbf{P}$ is a diagonal, singular matrix with ones in the indices corresponding to the orbitals we wish to project onto and zeros elsewhere.

Table 9 shows the calculated MAE for a number of Fe-based L$1_{0}$ alloys. Atomic structures as well as reference values calculated with SIESTA and VASP using the FT method are from Ref. Blanco-Rey et al., 2018. 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. Masuda et al., 2017. We use periodic BCs in the transverse directions. The calculated interfacial anisotropy constant $K_{1}=\textrm{MAE}/(2A)$, where $A$ is the cross-sectional area, is $K_{1}=1.41$ mJ/m2, in close agreement with a previous reported valueMasuda et al. (2017) of $K_{1}=1.40$ mJ/m2. From the atom-projected MAE (black circles) it is clear that the interface Fe atoms favor perpendicular MAE (since $\textrm{MAE}>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 $\mathbf{k}$-points, and to seamlessly connect the three regions, the same $\mathbf{k}$-point sampling is used in the transverse directions for all three regions. In the transport direction, the central-region 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 $\mathbf{k}$-point grid is therefore needed in the transport direction.

The left and right electrodes are modelled in their ground states with chemical potentials $\mu_{\text{L}}$ and $\mu_{\text{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

[TABLE]

where $V_{\mathrm{bias}}$ is the bias applied on the electrodes. It is clear that $\mu_{\text{R}}=\mu_{\text{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.

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,

[TABLE]

We now focus on the contribution from the extended states of the left ($n_{\text{L}}$) and right ($n_{\text{R}}$) electrodes, and delay the discussion of bound states ($n_{\text{B}}$) for later. The former may be obtained by calculating the scattering states incoming from the left ($\psi_{\alpha}^{\text{L}}$) and right ($\psi_{\alpha}^{\text{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 calculated by summing up the occupied scattering states,

[TABLE]

where $f(x)=(1+\mathrm{e}^{x})^{-1}$ is the Fermi–Dirac distribution.

XI.1 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.Brandbyge et al. (2002)

The electron density is given in terms of the electron density matrix. We split the density matrix into left and right contributions,

[TABLE]

The left contribution is calculated using the NEGF method asBrandbyge et al. (2002)

[TABLE]

where

[TABLE]

is the spectral density matrix, expressed in terms of the retarded Green’s function $G$ and the broadening function $\Gamma^{\text{L}}$ of the left electrode,

[TABLE]

which is given by the left electrode self-energy $\Sigma^{\text{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_{\text{L}}$.

Similar equations exist for the right density matrix contribution. The next section describes the calculation of $G$ and $\Sigma$ in more detail.

We note that the implemented NEGF method supports spintronic device simulations, using a noncollinear electronic spin representation, and possibly including spin-orbit coupling. This enables, for example, studies of spin-transfer torque driven device physics.Nikolić et al. (2018)

XI.2 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,

[TABLE]

where $\delta_{+}$ is an infinitesimal positive number.

Calculation of $G$ at a specific energy $\varepsilon$ 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,Petersen et al. (2008) which is $\mathcal{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,Sanvito et al. (1999); Sancho et al. (1985); Sørensen et al. (2008, 2009) where our preferred algorithm use the recursion method of Ref. Sancho et al., 1985, 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.

XI.3 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. Brandbyge et al., 2002 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

[TABLE]

where

[TABLE]

where $\rho^{\text{B}}$ is the density of states of any bound states in the central region. Equivalently, we could write the density matrix as

[TABLE]

where $\mathrm{L}$ and $\mathrm{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,Brandbyge et al. (2002) where (58) and (61) are weighted such that the main fraction of the integral is obtained from the equilibrium parts, $D_{\mathrm{eq}}^{\text{L}}$ and $D_{\mathrm{eq}}^{\text{R}}$, which are usually much more accurate than the non-equilibrium parts, due to the use of high-precision contour integration. We have

[TABLE]

where $W^{\text{L}}$ and $W^{\text{R}}$ are chosen according to Ref. Brandbyge et al., 2002, i.e. such that at each site, the equilibrium part of the density matrix gives the largest contribution and $W^{\text{L}}+W^{\text{R}}=1$.

XI.4 Bound States

The non-equlibrium integrals, $\Delta_{\mathrm{neq}}^{\text{L}}$ and $\Delta_{\mathrm{neq}}^{\text{R}}$, 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 $\mu_{\text{L}}<\mu_{\text{R}}$, then a bound state with energy $\varepsilon_{\text{B}}<\mu_{\text{L}}$ will be included in both $D_{\mathrm{eq}}^{\text{L}}$ and $D_{\mathrm{eq}}^{\text{R}}$, but a bound state in the bias window, $\mu_{\text{L}}<\varepsilon_{\text{B}}<\mu_{\text{R}}$, will only be included in $D_{\mathrm{eq}}^{\text{R}}$. Thus, from (62) we see that the state will be included with weight 1 if $\varepsilon_{\text{B}}<\mu_{\text{L}}$ and only with a fractional weight if $\mu_{\text{L}}<\varepsilon_{\text{B}}<\mu_{\text{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.

XI.5 Spill-in Terms

Given the density matrix $D$, the electron density is obtained from the LCAO basis functions $\phi$:

[TABLE]

The Green’s function of the central region gives the density matrix of the central region, $D^{\text{CC}}$. However, to calculate the density correctly close to the central-region boundaries towards the electrodes, the terms involving $D^{\text{LL}}$, $D^{\text{LC}}$, $D^{\text{CR}}$, and $D^{\text{RR}}$ are also needed. These are denoted spill-in terms.Stradi et al. (2016)

QuantumATK implements an accurate scheme for including all the spill-in terms, both for the electron density and for the Hamiltonian integrals.Stradi et al. (2016) This gives additional stability and well-behaved convergence in device simulations.

XI.6 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,Todorov et al. (2000)

[TABLE]

where $N_{\text{L/R}}$ is the number of electrons contributed to the central region from the left/right electrode, and $E_{\mathrm{KS}}[n]$ is the KS total energy.

Due to the screening approximation, the central region will be charge neutral, and therefore $N_{\text{L}}+N_{\text{R}}=N$, where $N$ is the ionic charge in the central region. At zero bias ($\mu_{\text{L}}=\mu_{\text{R}}$), the particle term is constant, so that $N\mu_{\text{L}}=N\mu_{\text{R}}$, and is thus independent of atom displacements in the central region. However, at finite bias ($\mu_{\text{L}}\neq\mu_{\text{R}}$), the particle terms in $\Omega$ will affect the forces.

If one neglects current-induced forces,Lü et al. (2012); Todorov et al. (2014) as done in QuantumATK simulations, the force acting on atom $a$ at position ${\bf r}_{a}$ in the device central region is given by

[TABLE]

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.Brandbyge et al. (2002); Todorov et al. (2000); Zhang (2011)

XI.7 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 $\varepsilon$ is obtained from the retarded Green’s function,Haug and Jauho (2008)

[TABLE]

and the electrical current is given by the Landauer formula,

[TABLE]

XI.8 Inelastic Transmission and Inelastic Current

QuantumATK implements the lowest-order expansion (LOE) methodLü et al. (2014) 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, $\nabla H(\mathbf{r})$. Calculation of these derivatives are described in Section VIII.

First-principles calculation of $\nabla H(\mathbf{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 $\nabla H(\mathbf{r})$ can be obtained to a good approximation by using the $\nabla H(\mathbf{r})$ of the left electrode only.Gunst et al. (2017b)

From $\nabla H(\mathbf{r})$ of the central region we get the electron-phonon matrix elements in reciprocal space,Gunst et al. (2016)

[TABLE]

where the ($mn$)-sum runs over repeated unit cells in the supercell calculation of the Hamiltonian derivatives,Gunst et al. (2016) and the subscript 0 indicates that the derivatives are only calculated for atoms in the unit cell with index ‘0’. Moreover, $|\phi_{i}\;{\bf R}_{n}\rangle$ ($|\phi_{j}\;{\bf R}_{m}\rangle$) denotes the $i$($j$)’th LCAO basis orbital in the unit cell displaced from the reference cell by the lattice vector ${\bf R}_{n}$ (${\bf R}_{m}$), while ${\bf q}$ is the phonon momentum, and ${\bf v}_{\lambda,{\bf q}}$ is the mass-scaled mode vector of phonon mode $\lambda$ with frequency $\omega_{\lambda,{\bf q}}$.

Following Ref. Lü et al., 2014, we obtain the inelastic transmission functions for a finite transfer of momentum. From these we calculate the total electrical current, including inelastic effects.Vandenberghe et al. (2011); Gunst et al. (2017b) The complete formulas for the QuantumATK implementation can be found in Ref. Gunst et al., 2017b.

XI.8.1 Special Thermal Displacement Method

In Ref. Markussen et al., 2017 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.Gunst et al. (2017c)

To obtain the STD configuration, we first calculate the phonon eigenspectrum using the dynamical matrix of the central region. We consider only $\mathbf{q}=\textbf{0}$, since only relative displacements between atoms in the cell will be important, and to account for finite $\mathbf{q}$-vectors we will have to increase the cell size. The phonon modes are labeled by $\lambda$ with frequency $\omega_{\lambda}$, eigenmode vector ${\bf e}_{\lambda}$, and characteristic length $l_{\lambda}$.

The STD vector of atomic displacements is given byGunst et al. (2017c)

[TABLE]

where $s_{\lambda}$ denotes the sign of the first non-zero element in ${\bf e}_{\lambda}$, enforcing the same choice of “gauge” for the modes. The Gaussian width $\sigma$ is related to the mean square displacement $\langle{\bf u}_{\lambda}^{2}\rangle=l^{2}_{\lambda}(2n_{B}(\frac{\hbar\omega_{\lambda}}{k_{B}T})+1)=\sigma_{\lambda}^{2}(T)$ at temperature $T$, where $n_{\text{B}}$ is the Bose–Einstein distribution.

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 ${\bf u}_{\mathrm{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.

XI.9 Thermoelectric Transport

The thermoelectric figure of merit, $\mathrm{ZT}$, quantifies how efficiently a temperature difference (heat) can be converted into a voltage difference in a thermoelectric material,

[TABLE]

where $G_{\mathrm{e}}$ is the electronic conductance, $S$ the Seebeck coefficient, $T$ the temperature, and $\kappa=\kappa_{\mathrm{e}}+\kappa_{\mathrm{ph}}$ the summed electron and phonon heat transport coefficients. Following Ref. Markussen et al., 2009, 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, $\Pi$,

[TABLE]

where $I_{Q}=dQ/dT$ is the electronic contribution to the heat current. It is calculated in a similar way as the electronic current,Sivan and Imry (1986)

[TABLE]

where $\mu=(\mu_{\mathrm{L}}+\mu_{\mathrm{R}})/2$ is the average chemical potential, and the difference to (67) is the inclusion of the factor $(\varepsilon-\mu)$ 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.

XI.10 Photocurrent

QuantumATK allows for calculating photocurrent using first-order perturbation theory within the first Born approximation.Henrickson (2002); Chen et al. (2012); Zhang et al. (2014) In brief, the electron-light interaction is added to the Hamiltonian,

[TABLE]

where $\hat{H}_{0}$ is the Hamiltonian without the electron-light interaction, $e$ the electron charge, $m_{0}$ the free-electron mass, $\hat{{\bf p}}$ the momentum operator, and ${\bf A}_{\omega}$ the electromagnetic vector potential from a single-mode monocromatic light source with frequency $\omega$.

The first-order coupling matrix is

[TABLE]

where $|j\rangle$ is an LCAO basis function.

The first Born electron-photon self-energies are

[TABLE]

where $\varepsilon^{\pm}=\varepsilon\pm\hbar\omega$, and $N$ is the number of photons. The Green’s function including electron-photon interactions to first order is then

[TABLE]

where ${\bf G}_{0}^{r,>,<}$ denote the non-interacting Green’s functions, and ${\bf\Sigma}^{>/<}_{\text{L,R}}$ are the lesser and greater self-energies due to coupling to the electrodes. The current in electrode $\alpha$ (left or right) with spin $\sigma$ is calculated as

[TABLE]

where the effective transmission coefficients are given byZhang et al. (2014)

[TABLE]

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-band-gap materials such as silicon.Palsgaard et al. (2018b)

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.

XII.1 Bulk DFT and Semi-Empirical Simulations

For bulk DFT-LCAO calculations, the basic unit of computational work to distribute in parallel is a single $\mathbf{k}$-point. QuantumATK uses the message passing interface (MPI) protocol to distribute such work units as individual computing processes on individual, or small groups 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 considered a 64-atom SiGe random-alloy supercell with $N_{\text{k}}=32$ $\mathbf{k}$-points. In this case, 2 full nodes, $N_{\text{n}}=2$, with 32 cores in total ($N_{\text{c}}=N_{\text{n}}\times 16=32$), yields full MPI parallelization over $\mathbf{k}$-points. The PW calculations were done using a blocked generalized Davidson algorithmDavidson (1975); Morgan and Scott (1986) to iteratively diagonalize the Hamiltonian matrix, which in the QuantumATK implementation parallelizes the computational work over both $\mathbf{k}$-points and plane waves. The LCAO calculations use the LAPACKAnderson et al. (1999) (when $N_{\text{c}}/N_{\text{k}}\leq 1$) or ELPAMarek et al. (2014) (when $N_{\text{c}}/N_{\text{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 $\mathbf{k}$-points in the DFT calculation (when $N_{\text{c}}/N_{\text{k}}>1$). On the contrary, MPI parallelization over both $\mathbf{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 $\mathbf{k}$-point.

XII.2 DFT-NEGF Device Simulations

As discussed in Section XI.3, the NEGF equilibrium density matrix at a single $\mathbf{k}$-point is obtained from integrating the spectral density matrix over $M_{\varepsilon}$ energy points on a complex contour. This integral must be performed at all transverse $\mathbf{k}$-points in the 2D Brillouin zone of the device central region, yielding $N_{k}\times M_{\varepsilon}$ generalized contour points. Each of these constitute a unit of computational work in equilibrium NEGF calculations, equivalent to $\mathbf{k}$-point parallelization in DFT calculations for periodic bulks.

Since we typically have $M_{\varepsilon}=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 transverse 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 MPI-only 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.

XII.3 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. Schneider et al., 2017, the MPI distribution of ATK-ForceField workload is implemented via functionality from the Tremolo-X MD package,Hamaekers et al. (2019) 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 SiO2 with 1 million atoms, using a force field from Pedone et al.Pedone et al. (2006) 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.

XIII NanoLab Simulation Environment

XIII.1 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. QuantumATK 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,The HDF Group (2019) 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 QuantumATK simulation can be stored efficently in a single file.

XIII.2 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,Kresse and Hafner (1993) Quantum ESPRESSO,Giannozzi et al. (2009) ORCA,Neese (2012) GPAW,Enkovaara et al. (2010) and CASTEPClark et al. (2005) codes.

The plugin concept also allows for many specialized functions, for example specialized Builder tools like surface builders, interface builders,Stradi et al. (2017a) NEB setups,Smidstrup et al. (2014) 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.

XIII.3 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 QuantumATK reference manual is extracted directly from the Python source code using the Sphinx documentation generator.Brandl (2019) The reference manual is available from an online platformSynopsys, Inc. (2019) 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

XIV.1 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 configurationsStradi et al. (2017b); Zhong and et al. (2016) are simulated by including dielectric regions and electrostatic gates, see Section IV.5.

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.Fiori et al. (2014); idr (2017) 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,Mueller et al. (2010) 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 MoTe2/SnS2 heterojunction.Szabó et al. (2015) We consider two different contact schemes by including atomistic metallic contacts: In the symmetrically contacted (SC) $\mathrm{M_{D}}$/MoTe2/SnS2/$\mathrm{M_{S}}$ device, Au is used for both the source (MS) and drain (MD) metallic contacts, whereas in the asymmetrically contacted (ASC) device, we set $\mathrm{M_{D}=Al}$ and $\mathrm{M_{S}=Au}$, in order to have a rather large work function difference ($\Delta\Phi$) between $\mathrm{M_{D}}$ and $\mathrm{M_{S}}$.Singh-Miller and Marzari (2009) In both devices, the metallic contacts to MoTe2 and SnS2 are represented by $\langle 110\rangle$-oriented 4-layer slabs.

The device configurations were constructed from the optimized structures of the Au(110)/MoTe2 and Au(110)/SnS2 electrodes, and the interlayer distance in the overlap region was set to 3.1 Å. Following Ref. Szabó et al., 2015, the devices were encapsulated in a high-$\kappa$ dielectric region (HfO2, $\kappa=25.0$), and a thin low-$\kappa$ dielectric region (h-BN, $\kappa=6.0$) was placed above the “exposed” MoTe2 region that is not contacted or forms part of the overlap region, hereafter denoted $\mathrm{E(MoTe_{2})}$. Electrostatic top and bottom gates were defined outside the high-$\kappa$ dielectric region, covering the overlap and half of the $\mathrm{E(MoTe_{2})}$ and $\mathrm{E(SnS_{2})}$ regions. The ASC device was constructed by replacing the Au atoms in the left electrode with Al atoms, with no further structural optimization.Che Additional computational details are given in Appendix A.

To study the impact of the contact asymmetry on the device characteristics, the reverse-bias $I_{\mathrm{DS}}$-$V_{\mathrm{GS}}$ curves (the transconductance) were simulated for both devices and for two values of the drain-source voltage, $V_{\mathrm{DS}}=-0.2$ V and $V_{\mathrm{DS}}=-0.4$ V, by grounding the top gate and by sweeping the bottom gate. The same physical picture emerges for both values of $V_{\mathrm{DS}}$, and we discuss here only the results obtained for $V_{\mathrm{DS}}=-0.4$ V. The $I_{\mathrm{DS}}$-$V_{\mathrm{GS}}$ curves in Fig. 15 show that the drain-source current is higher in the SC device than in the ASC device across the entire range of gate-source voltages. However, in the SC device, $I_{\mathrm{DS}}$ increases only by a factor of $\sim$10, from $I_{\mathrm{DS}}(V_{\mathrm{GS}}=0.05\ \mathrm{V})=4.28\times 10^{-8}\ \mathrm{A\ cm^{-1}}$ to $I_{\mathrm{DS}}(V_{\mathrm{GS}}=0.6\ \mathrm{V})=1.29\times 10^{-6}\ \mathrm{A\ cm^{-1}}$. Conversely, in the ASC device, $I_{\mathrm{DS}}$ increases by about six orders of magnitude in the same $V_{\mathrm{GS}}$ range, from $I_{\mathrm{DS}}(V_{\mathrm{GS}}=0.05\ \mathrm{V})=4.09\times 10^{-15}\ \mathrm{A\ cm^{-1}}$ to $I_{\mathrm{DS}}(V_{\mathrm{GS}}=0.6\ \mathrm{V})=1.74\times 10^{-9}\ \mathrm{A\ cm^{-1}}$.

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, $\mu_{\mathrm{D}}$ and $\mu_{\mathrm{S}}$, are aligned on a common energy scale. On the other hand, the interaction between the metallic contact and $\mathrm{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 ($\Delta V^{\mathrm{H}}$) in the two devices at $V_{\mathrm{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 $\mathrm{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 $\Phi^{M_{\mathrm{S}}}>\Phi^{M_{\mathrm{D}}}$, 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, $[\mu_{\mathrm{D}}-\mu_{\mathrm{S}}]\pm k_{\mathrm{B}}T=\Delta\mu\pm k_{\mathrm{B}}T$, is strongly inhomogeneous along the channel, as the conduction bands (CBs) of MoTe2 and SnS2 are pinned to $\mu_{\mathrm{D}}$ and $\mu_{\mathrm{S}}$, respectively (see Fig. 16(c,d)). This results in a vanishing DOS in the $\mathrm{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 $\mu_{\mathrm{D}}$, and (ii) the effect of the local electric field, which bends and depletes even more the CBs, moving 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 $\Delta V^{\text{H}}$ and of the PLDOS. The DOS within $\Delta\mu\pm k_{\mathrm{B}}T$ in the $\mathrm{E(MoTe_{2})}$ region, described in terms of an effective barrier $\phi_{\mathrm{MoTe_{2}}}$, ultimately determines the reverse-bias current in the channel. In the SC device, $\phi_{\mathrm{MoTe_{2}}}$ is lower for the case $V_{\mathrm{GS}}=0$ V, and depends only weakly on $V_{\mathrm{GS}}$, as shown in Fig. 16(e). This results in a higher absolute value of $I_{\mathrm{DS}}$, and in a lower variation of $I_{\mathrm{DS}}$ with $V_{\mathrm{GS}}$. Conversely, in the ASC device, $\phi_{\mathrm{MoTe_{2}}}$ is higher at comparable values of $V_{\mathrm{GS}}$, and varies appreciably when $V_{\mathrm{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 $\mathrm{M_{D}}$/MoTe2/SnS2/$\mathrm{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.

XIV.2 Phonon-limited Mobility of Metals

The continued downscaling of nanoelectronics makes the metal interconnects an increasingly critical part of transistor designs.Josell et al. (2009) 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.5. 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 ($\mathbf{k}$- and $\mathbf{q}$-space), and we know of only few studies of the EPC in metals that includes a full integration.Bauer et al. (1998); Gall (2016); Markussen et al. (2017) 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\times 20\times 20$ $\mathbf{q}$-points and tetrahedron integration. In addition, we apply the two-step procedure, where a $\mathbf{k}$-space isotropic but energy-dependent scattering rate is used to efficiently evaluate the resistivity.

Figure 17 shows the DFT results for the temperature-dependent phonon-limited resistivity of bulk Cu, Ag, and Au (Debye temperatures of 347, 227 and 162 K,Stewart (1983) respectively). The resistivity increases with temperature as the phonon occupation increases, and becomes linearly dependent on temperature above the Debye temperature.

Table 10 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/\rho(T)=\frac{1}{3}e^{2}v_{\text{F}}^{2}\tau(T)n(\varepsilon_{\text{F}})$. In the three FCC metals considered here, the Fermi velocity, $v_{\text{F}}$, and the DOS, $n(\varepsilon_{\text{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.

XIV.3 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 LiFePO4, a promising cathode material of the olivine family for Li-ion batteries.Islam and Fisher (2014); Zhang (2011) In this class of materials, the olivine scaffold provides natural diffusion channels for the $\mathrm{Li^{+}}$ ions, which have been shown to diffuse via a hopping mechanism, preferentially through [010]-oriented channels.Boulfelfel et al. (2011); Islam and Fisher (2015) 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 $\mathrm{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.English and Waldron (2015); Rungger et al. (2010)

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,Kima et al. (2011) or by running semi-classical dynamics on precalculated potential-energy surfaces based on DFT.Kahle et al. (2018)

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 $\mathbf{F}^{\prime}_{a}=Q_{a}\,\mathbf{D}$ on the $a$-th ion with formal charge $Q_{a}$ and $\mathbf{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.4) by defining the time-dependent formal charge $Q^{\prime}_{a}(t)$ as

[TABLE]

where $\Delta Q_{a}$ describes the time-dependent fluctuation. In principle, $\Delta Q_{a}$ can be defined arbitrarily, provided that charge neutrality is maintained in the system. In the present case, we chose a simple definition,

[TABLE]

where $\mathrm{Q}_{a}^{\mathrm{DFT}}(t)$ and $\mathrm{Q}_{a}^{\mathrm{DFT,ref}}$ are the time-dependent charge of the $i$-th atom obtained from a DFT calculation for the MD configuration at time $t$, and a time-independent 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 $\mathrm{Q}_{a}^{\mathrm{FF}}$ and $\Delta Q_{a}(t)$ does not constitute an issue, since the charge fluctuations during the dynamics are of the order $\Delta Q_{a}\sim 0.1\ e^{-}$.

We have applied this multi-model approach to investigate the interplay between Brownian and drift components of the diffusion of $\mathrm{Li^{+}}$ ions along the [010] channels in LiFePO4 in the presence of an applied displacement field. The system was described by a $1\times 2\times 1$ LiFePO4 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.,Pedone et al. (2006) which has been shown to describe qualitatively correctly the geometry and transport properties of olivine materials.Kutteh and Avdeev (2014) 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 $\mathbf{D}=[0,D_{y},0]$, with $0.0\ \text{V/\AA}\leq D_{y}\leq 0.3\ \text{V/\AA}$. 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 $\Delta Q_{a}(t)$ was recalculated every 100 MD steps, see (84), with $\mathrm{Q}_{a}^{\mathrm{DFT}}(t)$ and $\mathrm{Q}_{a}^{\mathrm{DFT,ref}}$ obtained from Mulliken population analysis. Further computational details are given in Appendix A.

Figure 18(a,c) shows the average displacement $\langle y_{\mathrm{Li^{+}}}\rangle$ of the $\mathrm{Li^{+}}$ ions along the $y$ Cartesian direction, that is, along the [010] channels of the FePO4 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 $\mathrm{Li^{+}}$-ion displacement remains constant at $\langle y_{\mathrm{Li^{+}}}\rangle=4.67\pm 0.11$ Å during the entire simulation, indicating the absence of hopping events. At 1000 K, the situation is rather similar, as $\langle y_{\mathrm{Li^{+}}}\rangle$ increases only slightly from an initial value of $7.12\pm 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\pm 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 $\langle y_{\mathrm{Li^{+}}}\rangle$ over time. This indicates that, at both temperatures, $\mathrm{Li^{+}}$ hopping due to Brownian motion is a rare event.

Applying an increasingly stronger displacement field leads to a progressive increase in the $\mathrm{Li^{+}}$ hopping probability. At 300 K, the average $\mathrm{Li^{+}}$-ion displacement increases steadily from the beginning of the MD run for $D_{y}\geq 0.20$ V/Å, indicating that, for these values of $D_{y}$, $\mathrm{Li^{+}}$ hopping is primarily due to field-induced drift. The $\mathrm{Li^{+}}$ ions accelerate until they reach a constant velocity, as shown by the tendency of the $\langle y_{\mathrm{Li^{+}}}\rangle$ 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 $\mathrm{Li^{+}}$ ion diffusion due to increased Brownian motion.Kutteh and Avdeev (2014); Islam and Fisher (2015) However, in the present case we find that the $\mathrm{Li^{+}}$ ions move less at 1000 K than at 300 K. For comparable values of $D_{y}$, the $\langle y_{\mathrm{Li^{+}}}\rangle$ vs. time curve has a smaller slope at 1000 K than those calculated at 300 K. The reason is that collision events of the $\mathrm{Li^{+}}$ ions with the LiFePO4 lattice, where phonons are considerably more excited at higher temperatures than at room temperature, limits the effective velocity of the $\mathrm{Li}^{+}$ ions.

This is evident by comparing the LiFePO4 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 LiFePO4 structure is relatively unperturbed. Consequently, the $\mathrm{Li^{+}}$ ions are able to travel through the [010] channels with relatively few scattering events with the LiFePO4 lattice. Conversely, at 1000 K, the LiFePO4 structure is significantlty perturbed, leading to a high probability of collisions between the $\mathrm{Li^{+}}$ ions and the olivine lattice.

In summary, we have studied the diffusion of $\mathrm{Li^{+}}$ in olivine LiFePO4, 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 $\mathrm{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 FePO4 lattice.

XIV.4 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.Taur et al. (1997) 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.Li and Fu (2012) Atomic-scale simulations may be used to calculate such parameters from first principles if experimental values are not available, including composition dependence.Khomyakov et al. (2015) 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) approachZunger et al. (1990) 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 Si1-xGex supercells in the full range of compositions, $0\leq x\leq 1$, by calculating composition-dependent lattice constants and band energies. The PPS-PBE method, Smidstrup et al. (2017) discussed in Section IV.4, 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.Khomyakov et al. (2015)

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. Zunger et al., 1990; van de Walle, 2009. 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 Si1-xGex 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 Si1-xGex fundamental band gap changes character at $x\sim 0.85$. The PPS-PBE and HSE06 predictions of the transition point are $x\sim 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 11, including the band-gap compositional bowing parameters. The PPS-PBE band gaps are in good agreement with room-temperature experiments (within $\sim$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 11. The bowing parameter of 0.034 Å is overestimated by $\sim$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 experimental 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 Si1-xGex random alloys, suggesting that SQS provides an accurate and efficient approach to random-alloy 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 QuantumATK 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 $\mathrm{Li^{+}}$-ion drift in a battery cathode material. We also presented applications of QuantumATK for simulating 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,Stokbro et al. (2010); Schneider et al. (2017); Stradi et al. (2016); Gunst et al. (2016, 2017c); Smidstrup et al. (2017); Palsgaard et al. (2018b); Petersen et al. (2008); Smidstrup et al. (2014); Stradi et al. (2017a) 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.

Acknowledgements.

Authors acknowledge funding from the European Union’s Horizon 2020 research and innovation programme under grant agreements No 713481 (SPICE), No 766726 (COSMICS), and No 723867 (EMMC-CSA), as well as funding from the Quantum Innovation Center (QUBIZ) and the Lundbeck Foundation (R95-A10510). CNG is sponsored by the Danish National Research Foundation (DNRF103).

Appendix A Computational Details

In the simulations presented in Fig. 1, we have considered non–crystalline a-$\mathrm{Al_{2}O_{3}}$ structures with a constant density of 2.81 g/cm3. 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,Pedone et al. (2006) 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–PackMonkhorst and Pack (1976) (MP) $\mathbf{k}$-point density of 3–4 Å. For systems with sizes between 240 and 960 atoms, 2 processes/$\mathbf{k}$-point was used, whereas for the 1920-atom system, 16 processes/$\mathbf{k}$-point was used.

For the DFT-NEGF device simulations presented in Section XIV.1, we used the PBE density functional with SG15-Medium (FHI-DZP) combinations of pseudopotentials and basis sets for MoTe2 and SnS2 (Au and Al). The real-space cutoff energy was 2721 eV (100 Ha), and MP $\mathbf{k}$-point grids of $12\times 1\times 100$ and $12\times 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.3, 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 FePO4 BZ was sampled using a $3\times 3\times 2$ MP $\mathbf{k}$-point grid.

For the electronic-structure calculations for SiGe random alloys presented in Section XIV.4, we used a $3\times 3\times 3$ MP $\mathbf{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.