J. Phys. Soc. Jpn. 91, 091011 (2022) [13 Pages]
SPECIAL TOPICS: Hyper-Ordered Structures: Recent Progress and Future Perspectives

Large-Scale DFT Methods for Calculations of Materials with Complex Structures

+ Affiliations
1International Centre for Materials Nanoarchitectonics (WPI-MANA), National Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0044, Japan2PRESTO, Japan Science and Technology Agency (JST), Kawaguchi, Saitama 333-0012, Japan3London Centre for Nanotechnology, University College London, 17-19 Gordon St., London WC1H 0AH, United Kingdom4Department of Physics & Astronomy, University College London, Gower St., London WC1E 6BT, United Kingdom

Large-scale density functional theory (DFT) calculations provide a powerful tool to investigate the atomic and electronic structure of materials with complex structures. This article reviews a large-scale DFT calculation method, the multi-site support function (MSSF) method, in the CONQUEST code. MSSFs are linear combinations of the basis functions which belong to a group of atoms in a local region. The method can reduce the computational time while preserving accuracy. The accuracy of MSSFs has been assessed for bulk Si, Al, Fe and NiO and hydrated DNA, which demonstrate the applicability of the MSSFs for varied materials. The applications of MSSFs on large systems with several thousand atoms, which have complex interfaces and non-periodic structures, indicate that the MSSF method is promising for precise investigations of materials with complex structures.

©2022 The Author(s)
This article is published by the Physical Society of Japan under the terms of the Creative Commons Attribution 4.0 License. Any further distribution of this work must maintain attribution to the author(s) and the title of the article, journal citation, and DOI.
1. Introduction

Atomic structure, electronic structure and the properties of materials correlate with each other strongly. The relationship of materials properties with local structures such as point defects in crystals and interatomic distances in glassy materials have been investigated widely for many years. Recently, wider structures ordered on the nanoscale, such as ring structures in amorphous glass,1) ionic and molecular positions in biomaterials,2) interfaces between metallic nanoparticle catalysts and substrates,3) and composite dopants and defects in crystals,4) have also been focused on as hyper-ordered structures, which could have significant influence on materials properties. The recent improvements in experimental measurements and computation techniques have enabled us to investigate such nano-scale complex structures.

For computation, density functional theory (DFT) calculations have been a powerful tool to investigate atomic and electronic structures of condensed-phase materials and molecules.5) However, because of the high computational cost of conventional DFT calculation methods, the system size treated by DFT has been limited, up to a thousand atoms in most cases. Therefore, efficient DFT calculation methods are desirable to treat nano-scale structures. Several methods have been proposed to overcome the size limitation of DFT calculations, and we have also proposed the large-scale DFT code, CONQUEST.610) There are two important methods which enable us to perform large-scale calculations with CONQUEST, the linear-scaling order-N, or \(O(N)\), method7,8) and the multi-site support function method.11,12) Systems as large as one million atoms can be treated by the \(O(N)\) method with a massively parallel supercomputer.9,13) The MSSF method can be used not only to decrease the computational cost but also to improve the computational accuracy, and also enables us to perform stable large DFT calculations on metallic systems.

In the next section, we first briefly review several large-scale DFT calculation methods, and then provide more details of our large-scale calculation techniques in the CONQUEST code, especially the MSSF method. In Sect. 3, by showing several examples, we demonstrate the applicability of the MSSF method to large systems with complex structures. The final section provides the conclusion.

2. Large-Scale DFT Calculation Methods
Large-scale DFT calculations

The DFT total energy E is a functional of electron density, which is the diagonal part of the density matrix. The density matrix can be written in terms of the Kohn–Sham (KS) one-electron orbitals ψ,14,15) \begin{equation} \rho (\mathbf{r},\mathbf{r}') = \sum_{n} f_{n}\psi_{n}(\mathbf{r})\psi_{n}(\mathbf{r}')^{*}, \end{equation} (1) where n runs over KS orbitals and \(f_{n}\) is the occupation number of nth KS orbital. The KS orbitals ψ and their energies ε are calculated as the eigenvectors and eigenvalues of the KS equation, \begin{align} H^{\text{KS}}\psi_{n}(\mathbf{r}) &= \left[- \frac{\hbar} {2m}\nabla^{2} + V_{\text{ext}}(\mathbf{r}) + V_{\text{H}}(\mathbf{r}) + V_{\text{XC}}(\mathbf{r}) \right]\psi_{n}(\mathbf{r}) \notag\\ &= \varepsilon_{n}\psi_{n}(\mathbf{r}). \end{align} (2) In the KS Hamiltonian \(H^{\text{KS}}\), the kinetic energy and the external potential \(V_{\text{ext}}\) are one-electron terms, and the Hartree potential \(V_{\text{H}}\) and the exchange–correlation potential \(V_{\text{XC}}\) are two-electron terms. The computational cost to solve the KS equation with conventional calculation methods such as direct diagonalization scales cubically with respect to the number of atoms N in a target system. This cubic scaling limits the system size for practical DFT calculations.

To overcome this limitation, several large-scale DFT calculation methods have been proposed. Since there are several review papers already which cover large-scale DFT calculation methods,16,17) here we briefly introduce the methods focusing on two key points, “how to solve the KS equation” and “what kind of functions are used to express the KS orbitals”. The parallelization efficiency of DFT codes, which we do not describe here, is also very important.

For linear scaling solution of the KS equation for large systems, there are several methods such as the divide-and-conquer (DC) method,1820) the orbital minimization method (OMM),21,22) the density matrix minimization (DMM) method,23,24) and the fragment molecular orbital (FMO) method.25) In the DC method, the target system is divided into small subsystems with some buffer regions whose electronic structures are calculated exactly, and then the subspace density matrices are combined to construct the density matrix of the whole target system. The DC method is used in many codes such as SIESTA26,27) and OpenMX.28) FMO is similar to the DC method, but the division is based on molecular fragments in proteins and biomolecules. The total energy of the whole system is calculated from the energy of fragments and pairs of fragments without solving for the molecular orbitals (MOs) of the whole system. Recently, a method to obtain molecular orbitals of the whole system from subspace MOs (FMO-LCMO) has been also proposed.29) FMO is implemented in GAMESS30) and ABINIT-MP.31) OMM and DMM are methods which minimize the total energy variationally as a functional of orbitals or density matrices, instead of solving KS equations directly. OMM is used in FEMTECK32) and SIESTA, and DMM is used in CONQUEST and ONETEP.33) There are also iterative calculation methods such as the Fermi-operator expansion method34) in BigDFT35) and the second-order trace-correcting (TC2) method36) in ErgoSCF.37) CP2K also uses an iterative \(O(N)\) method.38)

For the functions to represent the KS orbitals, plane-wave basis functions have been often used with the periodic boundary condition for solids and surfaces, while Gaussian- or Slater-type atomic-orbital (AO) basis functions have been popular for isolated systems such as molecules and clusters. For large-scale calculations, using local orbital functions to express the KS orbitals is crucial. Pseudo-atomic orbitals (PAOs) are atomic orbital functions constrained to go to zero at a cutoff, which avoids having long-range tails. PAOs are often used in large-scale DFT codes such as SIESTA,39,40) OpenMX,28,41) and CONQUEST.42,43) There are also basis functions defined on a regular grid similarly to plane waves, which can be systematically converged, such as B-spline functions in CONQUEST,44) periodic cardinal sine functions in ONETEP,45) and wavelet functions in BigDFT.46) CP2K uses Gaussian functions to describe KS orbitals and represents the electron density on a grid.47) There are also codes using the finite difference method, for example RSDFT48) and PARSEC.49)

Large-scale DFT calculations in CONQUEST

In this subsection, we explain more details about our large-scale DFT code CONQUEST. In CONQUEST, the density matrix in Eq. (1) is expressed by using local orbital functions ϕ, called “support functions”, so that the density matrix \(\mathbf{K}\) is represented in the support function basis, \begin{equation} \rho (\mathbf{r},\mathbf{r}') = \sum_{i\alpha,j\beta} \phi_{i\alpha} (\mathbf{r})K_{i\alpha,j\beta} \phi_{j\beta} (\mathbf{r}')^{*}, \end{equation} (3) where α and β are the indices of the support function and i and j are the indices of the atoms which the support functions belong to. The support functions are linear combinations of given basis functions, \begin{equation} \phi_{i\alpha} = \sum_{\mu} c_{i\alpha,i\mu} \chi_{i\mu}, \end{equation} (4) where \(\mathbf{c}\) is the linear-combination coefficients and χ is the μth basis functions of atom i.

Two kinds of basis functions, B-spline (blip) functions44) and PAO functions,42,43) are available in CONQUEST. Blip functions are finite-element functions akin to plane-wave basis functions, i.e., the accuracy of the support functions can be improved systematically by making the blip functions finer. However, the optimization of the linear-combination coefficients for fine blip functions can be computationally expensive. On the other hand, the computational cost of PAOs is much cheaper than blip functions. PAOs are atomic orbital functions consisting of numerical radial functions R and spherical harmonic functions Y, \begin{equation} \chi_{i\mu} (\mathbf{r}) = \chi_{i\zeta lm}(\mathbf{r}) = R_{i\zeta l}(r)Y_{lm}(\Omega). \end{equation} (5) ζ is the index of the radial function, and l and m are the azimuthal and magnetic quantum numbers. Ω is the solid angle of \(\mathbf{r}\). The accuracy of support functions with PAOs is generally improved by using multiple radial functions for each Y (i.e., multiple-ζ functions) and adding Y functions with larger angular-momentum numbers to represent polarization (i.e., polarization functions), although systematic improvement is not guaranteed.

When we use PAOs as the support functions without any contraction, which are called “primitive” PAOs, the coefficients c in Eq. (4) are 1 or 0. When we contract multiple-ζ PAOs to smaller number of support functions, the c values are optimized numerically by, for example, the conjugate gradient method. The number of support functions can be reduced to the single-ζ (SZ) size for each angular momentum functions by contraction, keeping the spatial symmetry of the primitive PAOs. For example, the primitive triple-ζ (TZ) PAOs and triple-ζ plus triple polarization (TZTP) PAOs are contracted to SZ- and SZP-sized support functions, respectively. Since the computational cost scales cubically to the number of support functions in both the exact diagonalization and the \(O(N)\) calculations, it is crucial for large systems to reduce the number of support functions by the contraction. CONQUEST uses norm-conserving pseudopotentials5052) so that the support functions are used to describe only the valence electrons.

For the KS equation solver, CONQUEST supports both exact diagonalization and the DMM to optimize electron density. With the diagonalization method, the computational cost scales cubically with the size of the electronic Hamiltonian, i.e., the number of atoms N. On the other hand, the computational cost of DMM is linear with the system size. In DMM, the electronic structure is optimized by minimizing the total energy E with respect to the auxiliary density matrix \(\mathbf{L}\), \begin{align} &\frac{\partial E}{\partial \mathbf{L}} = 6(\mathbf{LSH}- \mathbf{HLS}) - 4(\mathbf{LSLSH}+ \mathbf{LSHSL}+ \mathbf{HSLSL}). \end{align} (6) \begin{align} &\mathbf{K}= 3\mathbf{LSL}- 2\mathbf{LSLSL}, \end{align} (7) avoiding cubic-scaling exact diagonalization. The electronic KS Hamiltonian matrix \(\mathbf{H}\) and the overlap matrix \(\mathbf{S}\) are sparse because they are based on local orbital support functions. We need to set a spatial cutoff for \(\mathbf{L}\), so that \(\mathbf{L}\) also becomes a sparse matrix. It is guaranteed that the total DFT energy becomes equal to that by the exact diagonalization when we increase the \(\mathbf{L}\) cutoff to be infinite.53) Thus, Eqs. (6) and (7) are sparse-matrix multiplications whose computational cost are linear to the matrix size (i.e., system size).

Multi-site support functions in conquest
Method

Since the computational cost of both the exact diagonalization and the \(O(N)\) method scales cubically to the number of the support functions, reduction of the number of support functions is crucial for the large-scale calculations. To reduce the number of support functions without losing accuracy, we have recently proposed the multi-site support functions (MSSFs).11,12) MSSFs are defined as \begin{equation} \phi_{i\alpha} = \sum_{k}^{\textit{neighbors}}\sum_{\mu \in k}c_{i\alpha,k\mu} \chi_{k\mu}, \end{equation} (8) where μ runs over the PAOs of atom k which are the neighbor atoms of atom i within the cutoff distance \(r_{\text{MS}}\). In Eq. (8), the MSSFs of atom i consist of not only the PAOs of atom i itself but also the PAOs of atom i's neighbour atoms j, like local MOs, while the conventional support functions in Eq. (4) consist only of the PAOs of atom i. Since MSSFs are like local MOs, the number of MSSFs can be reduced to SZ size, which is enough to represent the ground-state accurately. For example, the SZ size for a Si atom is four, so that double-ζ plus polarization (DZP) PAOs of Si consisting of (2s, 2p, d) = 13 functions are contracted to four MSSFs, while TZTP PAOs consisting of (3s, 3p, 3d) = 27 functions are also contracted to four MSSFs. Although the numbers of the constructed MSSFs are the same, the accuracy of MSSFs from TZTP PAOs is higher than those from DZP because the number of degrees of freedom of \(\mathbf{c}\) in the MSSFs is larger. This means that we can improve the accuracy of MSSFs by increasing the number of primitive PAOs without increasing the number of MSSFs. The computational cost to determine \(\mathbf{c}\) increases as the number of the original PAOs increases, but the computational cost is not dominant when the whole system size is large, as shown in Sect. 3.1.1. Therefore, the MSSF method has advantages both for computational efficiency and accuracy. The \(\mathbf{H}\) and \(\mathbf{S}\) matrices for the whole system are reconstructed from PAO basis to MSSF basis by sparse matrix multiplication as \begin{equation} H_{i\alpha,j\beta} = \sum_{k,\mu} \sum_{p,\nu} c_{i\alpha,k\mu} H_{k\mu,p\nu} c_{p\nu,j\beta}^{*}. \end{equation} (9)

Ideally, the linear-combination coefficients \(\mathbf{c}\) are optimized numerically.12) Since the number of the coefficients of MSSFs in Eq. (8) is larger than that of the conventional onsite support functions in Eq. (4), starting the optimization from accurate initial values is desirable. The local filter diagonalization (LFD) method by Rayson et al.54,55) is a powerful way to construct accurate initial values from localized occupied MOs in subspaces. The LFD method is based on local orbitals ϕ which can express the occupied KS eigenstates accurately, \begin{equation} | \phi_{i\alpha} \rangle = \sum_{n}| \psi_{n} \rangle f(\varepsilon_{n})\langle \psi_{n} | t_{i\alpha} \rangle, \end{equation} (10) where \(\mathbf{t}\) is the αth trial vector localized on atom i and \(f(\varepsilon_{n})\) is the Fermi–Dirac function. In the LFD method, \(\psi_{n}\) is replaced with subspace MOs around atom i to restrict the range of \(\phi_{i\alpha}\). To obtain the subspace MOs, the subspace for each atom i is first defined with a cutoff range \(r_{\text{LFD}}\). Then the electronic Hamiltonian and the overlap matrix for the subsystem, \(\mathbf{H}_{\text{sub}}\) and \(\mathbf{S}_{\text{sub}}\), are constructed to calculate the eigenvectors (i.e., subspace MOs) \(\mathbf{C}_{\text{sub}}\) and eigenvalues \(\varepsilon_{\text{sub}}\). Using the subspace MOs, ϕ is calculated as the linear-combination of the PAOs in the subspace with the linear-combination coefficients \(\mathbf{c}\), \begin{equation} \mathbf{c}= \mathbf{C}_{\text{sub}}f(\varepsilon_{\text{sub}})\mathbf{C}_{\text{sub}}^{\text{T}}\mathbf{S}_{\text{sub}}\mathbf{t}. \end{equation} (11) In Eq. (11), \(f(\varepsilon_{\text{sub}})\) is calculated with the chemical potential of the subspace, which reduces the influence from the unphysical MOs in high-energy unoccupied region.

The calculation accuracy with ϕ will depend on the choice of \(\mathbf{t}\). We have chosen \(\mathbf{t}\) from the primitive PAOs, but our previous study showed that the dependence is not very significant.11) Note that \(r_{\text{LFD}}\) should be equal to or larger than \(r_{\text{MS}}\). Once \(\mathbf{c}\) is determined, the SCF calculations of the whole systems can be performed with the matrices transformed as in Eq. (9), and ρ is obtained for the \(\mathbf{c}\). Here, we can update the subsystem \(\mathbf{H}\) matrices and construct new \(\mathbf{c}\) by using the updated ρ. Thus, the update of ρ and \(\mathbf{c}\) is iterated until self-consistency of ρ is reached. Since this update procedure is not variational, the calculated energy sometimes fluctuates, especially when \(r_{\text{MS}}\) is small.12)

By optimizing \(\mathbf{c}\) numerically subsequent to the LFD calculation, not only is the accuracy of MSSFs improved but also the energy becomes variational.12) The numerical optimization of \(\mathbf{c}\) is performed with the energy gradient with respect to \(\mathbf{c}\), \begin{equation} \frac{\partial E}{\partial c_{i\alpha,k\mu}} = \frac{\partial E}{\partial \phi_{i\alpha}}\frac{\partial \phi_{i\alpha}}{\partial c_{i\alpha,k\mu}} = \frac{\partial E}{\partial \phi_{i\alpha}}\chi_{k\mu}. \end{equation} (12) The gradient with respect to \(\phi_{i\alpha}\) is calculated as \begin{equation} \frac{\partial E}{\partial \phi_{i\alpha}} = 4\sum_{\beta}[K_{\alpha \beta} \hat{H} + G_{\alpha \beta}] \phi_{j\beta}, \end{equation} (13) where \(\mathbf{G}\) is given as the energy-weighted density matrix \begin{equation} G_{\alpha \beta} = \sum_{n}f_{n}\varepsilon_{n}u_{n\alpha} u_{n\beta}^{*} \end{equation} (14) in diagonalization calculations (\(\mathbf{u}\) is the KS coefficients in the support function basis) and as \begin{equation} G_{\alpha \beta} = 3(\mathit{LHL})_{\alpha \beta} - 2(\mathit{LSLHL}+ \mathit{LHLSL})_{\alpha \beta} \end{equation} (15) in \(O(N)\) calculations. Detailed derivations of Eqs. (13)–(15) are provided in our previous papers.7,56)

Investigation of accuracy and efficiency

In this subsection, the accuracy of MSSF is investigated by checking the \(r_{\text{MS}}\) dependence and the effect of the numerical optimization. All calculations shown in the present paper were performed with the exact diagonalization, not with the \(O(N)\) method to avoid the error coming from the cutoff of \(\mathbf{L}\), and \(r_{\text{LFD}}\) was set to be equal to \(r_{\text{MS}}\). More detailed information about the computational conditions such as the number of k-points are found in the references.

First, we investigate the accuracy of the calculated energies for bulk Si and Al systems with the local density approximation (LDA)57) exchange–correlation functional. Figure 1 shows the deviations of the DFT total energies by the MSSFs from those by the primitive PAOs.11) The TZP PAOs (3s, 3p, d) consisting of 17 functions are contracted to four MSSFs for both Si and Al atoms. The energy deviations of the MSSFs decrease exponentially as \(r_{\text{MS}}\) increase. It is worth noting that the energy of the bulk Al converges smoothly with respect to the spatial cutoff \(r_{\text{MS}}\) although the bulk Al is metallic; this is because the cutoff \(r_{\text{MS}}\) is introduced only to the support functions, not to the wave function of the whole system.


Figure 1. (Color online) Deviation of total DFT energy per atom of bulk Si and Al of MSSFs with respect to the multisite range \(r_{\text{MS}}\) from the DFT energy of primitive TZP PAOs. (Data taken with permission from Ref. 11. © 2014 American Chemical Society.)

The comparison of the energy–volume (EV) curves of bulk Si with several \(r_{\text{MS}}\) is shown in Fig. 2.12) The MSSFs are constructed from triple-ζ plus double polarization (TZDP) PAOs. The results of the MSSFs constructed only by LFD and those of the MSSFs by LFD and subsequent numerical optimization are compared in the figure. The lattice constants \(a_{0}\) calculated by fitting the EV curve with the Birch–Murnanghan equation are summarized in Table I. When only LFD is used to determine \(\mathbf{c}\), MSSFs with \(r_{\text{MS}} = 2.6\) Å, which contain only the nearest neighbor atoms in \(r_{\text{MS}}\), show an error of about 1.0%, while MSSFs with larger \(r_{\text{MS}}\), 4.2 and 8.5 Å, show much smaller errors, 0.2 and 0.0%, respectively. \(r_{\text{MS}} = 4.2\) Å contains up to the second-nearest neighbor atoms. The numerical optimization improves the accuracy of the MSSFs: when we optimize \(\mathbf{c}\) after the LFD method, even the error of MSSFs with \(r_{\text{MS}} = 2.6\) Å is small, only 0.2%.


Figure 2. (Color online) EV curve of bulk Si calculated using MSSFs with LFD, MSSFs with numerical optimization, and primitive TZDP PAOs. The multi-site ranges \(r_{\text{MS}}\) [Å] are shown in parentheses. (Reproduced from Ref. 12 with permission from the PCCP Owner Societies.)

Data table
Table I. Lattice constants \(a_{0}\) of crystalline Si calculated with MSSFs using the LFD method, MSSFs using LFD plus numerical optimization, and primitive TZDP PAOs. The percentage deviations from the result by TZDP are also shown. (Reproduced from Ref. 12 with permission from the PCCP Owner Societies.)

To check the accuracy of atomic forces, the energies and forces have been investigated for a distorted benzene molecule in which a C–H pair has been shifted away from the center of the benzene ring by 0.5 Å. Table II shows the differences of the energies and the maximum forces by the MSSFs with several \(r_{\text{MS}}\) from those by primitive DZP PAOs. It was clearly found that the increase of \(r_{\text{MS}}\) improves the accuracy of both energies and forces. When optimizing MSSF coefficients, even MSSFs with \(r_{\text{MS}} = 1.6\) Å provide comparable accuracy with DZP PAOs for forces, within 0.05 eV/Å.

Data table
Table II. Differences of the total energies and forces of a distorted benzene molecule calculated using MSSFs with the LFD method and MSSFs with LFD plus numerical optimization, from those by the primitive TZDP PAOs. (Part of this table is reproduced from Ref. 12 with permission from the PCCP Owner Societies.)

Next, the accuracy of the electronic structure was investigated by checking the density of states (DOS) of a hydrated DNA system with 3,088 atoms.12) The PBE generalized gradient approximation (GGA) functional was used.58) The difference of the DOS calculated with the primitive PAOs and the MSSFs are presented in Fig. 3. The difference is smaller in the occupied states and low-energy unoccupied states than in the high-energy unoccupied states. The difference in the high-energy unoccupied states is larger with smaller \(r_{\text{MS}}\). The accuracy in the occupied states is improved by the numerical optimization of \(\mathbf{c}\), but the description of the unoccupied states becomes less accurate. This is because the linear-combination coefficients \(\mathbf{c}\) are optimized only for the occupied states. The unoccupied states can be improved by cooperating with the Sakurai–Sugiura method (SSM),59) which will be shown in Sect. 3.4.


Figure 3. (Color online) DOS of the hydrated DNA system using primitive (black, solid) and MSSFs with local filter diagonalization (blue, dotted) and numerical optimization (red, dashed). Deviations from the primitive functions are shown for the occupied states. (Reproduced from Ref. 12 with permission from the PCCP Owner Societies.)

The computational efficiency was also investigated for the hydrated DNA system.12) Table III summarizes the times required for matrix construction, diagonalization and gradient calculation with respect to the coefficients \(\mathbf{c}\) [i.e., Eq. (12)] with primitive DZP and MSSFs. For the matrix construction, the MSSFs require additional time to construct \(\mathbf{c}\) and to reconstruct the overlap and Hamiltonian matrices in the MSSF basis. This additional time increases as \(r_{\text{MS}}\) becomes larger. On the other hand, for the diagonalization, MSSFs reduces the time dramatically with any \(r_{\text{MS}}\): about 600 s with MSSFs compared to about 12,000 s with primitive DZP. \(r_{\text{MS}}\) does not affect the time for diagonalization because the number of MSSFs does not depend on \(r_{\text{MS}}\). The times for the gradient calculations are shorter than those for the matrix constructions and the diagonalizations. Because the reduction of the diagonalization time is much larger than the increase of the times for the matrix construction and the gradient calculations, the total computational time can be reduced significantly although we need to iterate the SCF and the gradient calculations until the coefficient optimization converges. Since the number of the iterations depends on the accuracy of the initial values of the coefficients, we expect that the present optimization method will be more efficient if we can use coefficients from previous steps, such as in molecular dynamics simulations and geometry optimizations.

Data table
Table III. Computational times [s] for matrix construction, diagonalization and gradient calculation with respect to the coefficients for the hydrated DNA system by MSSFs and primitive DZP PAOs. (Data taken from Ref. 12 with permission from the PCCP Owner Societies.)
Spin-dependent MSSFs

For the spin-polarized systems, there are two possibilities for how to determine MSSF linear-combination coefficients: determining the coefficients for spin-up and spin-down electrons individually; and using the same coefficients for spin-up and spin-down electrons by taking their average. The first method leads to spin-dependent MSSFs. Spin-dependent MSSFs are obtained by using spin-up and spin-down \(\mathbf{C}_{\text{sub}}\) in Eq. (11) (and subsequent numerical optimization to minimize total DFT energies with respect to both spin-up and spin-down coefficients). When we use spin-dependent MSSFs, not only the two-electron terms but also the one-electron terms in the KS Hamiltonian and the overlap matrix become spin dependent, so that additional memory for those matrices is required. The second method leads to spin-independent MSSFs so does not require the additional memory, but the accuracy might be lower than spin-dependent MSSFs.

Therefore, we compare the accuracy of spin-dependent and spin-independent MSSFs for spin-polarized systems: bulk bcc ferromagnetic Fe and cubic antiferromagnetic NiO. \(r_{\text{MS}}\) is set to be 8.5 Å for both Fe and NiO. Figure 4 shows the EV curves of the Fe and NiO, and the \(a_{0}\) values fitted from the EV curves are summarized in Table IV. In Fig. 4, the energy values and the curvature of the spin-dependent (SD) MSSFs are closer to those of the PAOs than spin-independent MSSFs are. \(a_{0}\) of the SD MSSFs is closer to that of PAOs than spin-independent MSSFs for Fe, while they are comparable for NiO. The magnetic moments \(\mu_{\text{B}}\) of bcc Fe are calculated to be 2.40, 2.26, and 2.34 by the primitive PAOs, spin-independent (SI) MSSFs and SD MSSFs, respectively. Figure 5 shows the calculated DOS of antiferromagnetic NiO. The differences of the DOS from that of PAO is also shown, in which SD MSSFs clearly reduce the differences for both spin-up and spin-down states. These comparisons indicate that the accuracy of the calculation is improved by considering spin-dependence of MSSFs.


Figure 4. (Color online) EV curve of (a) bcc ferromagnetic Fe and (b) cubic antiferromagnetic NiO calculated with primitive DZP PAOs, SI-MSSF, and SD-MSSF.


Figure 5. (Color online) DOS of cubic antiferromagnetic NiO calculated with (a) primitive DZP PAOs, (b) SI-MSSF, and (c) SD-MSSF. The absolute differences of the DOS from (a) are also shown for (b) and (c). Spin-up (red) and spin-down (blue) states are shown as positive and negative values, respectively.

Data table
Table IV. Lattice constants \(a_{0}\) [Å] of bcc ferromagnetic Fe and cubic antiferromagnetic NiO calculated using primitive DZP PAOs, SI-MSSF and SD-MSSF.
3. Applications of MSSFs on Large Systems

To demonstrate the actual applicability of MSSFs to large systems, in this section we show several examples of the applications of MSSFs to large complex structures: moiré graphene on the Rh(111) surface;60) the interfaces in YGaO3;61) and PbTiO3 films on SrTiO3.62) The applications to nonperiodic systems, hydrated DNA59) and metallic gold nanoparticles, are also shown. The LDA exchange–correlation functional is used for YGaO3 and PbTiO3 on SrTiO3, and the PBE functional is used for graphene on Rh(111), hydrated DNA and Au nanoparticles. Norm-conserving pseudopotentials have been used in the calculations with CONQUEST. The other detailed information about the computational conditions such as the number of k-points and the spatial range of PAOs are found in the corresponding reference papers. In several examples, plane-wave calculations for comparison have been performed with the PAW pseudopotential63) using the VASP software.64,65)

Graphene on Rh surface

There have been many reports showing the exotic properties of 2D materials. The structure of 2D materials and their electronic structures are often affected by the interactions with the substrates or interlayer interactions. The target system in Ref. 60 also shows the interesting property that a highly corrugated graphene layer grown on Rh(111) can be flattened by the intercalation of oxygen atoms. For this system, plane-wave DFT calculations demonstrated that strong interactions between the graphene layer and the substrate are decoupled when oxygen atoms are intercalated in the lowest moiré sites. Such DFT simulations have been performed to clarify the structural and electronic properties of 2D materials, but they were limited in the size of simulation cells and it was difficult to study the case of low concentration of oxygen atoms, large moiré structures, the effect of the edges and so on. It is expected that these obstacles can be overcome if accurate MSSF method can be applied to such 2D materials.

In Ref. 60, the accuracy and computational-time efficiency of the MSSF method with a system, graphene on Rh(111) surface (G/Rh), were investigated. First, the accuracies of the PAOs in the optimized structures of G/Rh and G/O/Rh having small supercells were confirmed by comparing to plane-wave basis functions. TZDP PAOs of carbon atoms, TZTP PAOs of oxygen atoms and DZP PAOs of rhodium atoms were contracted to form MSSFs with \(r_{\text{MS}} = 8.5\) Å. The DFT-D2 dispersion correction method66) was used to consider van der Waals interaction. At first, the structural parameters, \(a_{0}\) and bulk modulus \(B_{0}\), of graphene and bulk Rh were compared, as shown in Table V. Figures 6(a) and 6(b) show the DOS of the G/Rh(111) with and without inserted oxygen atoms. As shown in the table and the figure, the results by the PAOs and the MSSFs can be very close to those by the plane-wave basis functions.


Figure 6. (Color online) DOS for graphene on a Rh(111) substrate (460 atoms) with no oxygen atom (\(\theta_{\text{O}} = 0\)) and with 12 oxygen atoms in the interface (\(\theta_{\text{O}} = 1/2\)) calculated with plane waves (black), PAOs (green), and multi-site functions (blue dashed). The red lines in the lower panels represent the DOS difference between the PAOs and the MSSF calculations. Fermi-levels are set to be zero (black dashed). (c) Structure of corrugated graphene on Rh(111) with 3088 atoms. (Reproduced with permission from Ref. 60. © 2018 IOP Publishing.)

Data table
Table V. Lattice parameters \(a_{0}\) and bulk moduli \(B_{0}\) of graphene and bulk rhodium calculated with plane-wave basis functions and MSSFs. (Data taken from Ref. 60.)

For large-scale calculations, the computational times by the primitive PAOs and the MSSFs for the G/Rh(111) consisting of 3,088 atoms [shown in Fig. 6(c)] are compared in Table VI. These large systems were too computationally expensive to be treated with plane-waves. The calculations were performed with the supercomputer SGI ICE X [Intel Xeon E5-2680V3 (12 cores, 2.5 GHz) \(\times\) 2 and 128 GB memory per node] at NIMS. In the case of 108 processes, the computational time with the MSSFs is about 18 times shorter than that with the PAOs. Although additional time to construct the contraction coefficients \(\mathbf{c}\) is needed (the time is included in “Matrix construction” in Table VI), the time for the diagonalization of the electronic Hamiltonian of the whole system, which is conventionally dominant in the total computational time, is dramatically reduced. This is because the computational cost of the diagonalization scales cubically with the number of support functions, which is reduced dramatically by the MSSF method. When the number of cores increases eight times, the computational time with MSSFs is reduced from 2156 to 571 s, almost a quarter. Although this is not ideal scaling, it demonstrates that a large speedup can still be achieved by parallelization. It should be also emphasized that large-scale DFT calculations are available with the MSSF method even for metallic systems.

Data table
Table VI. Computational times of an SCF step with PAOs and MSSFs for graphene on Rh(111) surface with 3088 atoms. (Data taken from Ref. 60.)
Interfaces in ferroelectric YGaO3

The next example is an investigation of topological defects in ferroelectric YGaO3.61) Ferroelectric domain walls are attracting broad attention for next-generation nanoelectronics. Although the basic properties of simple ferroelectric domain walls can be well described by small DFT calculations, complex domain patterns could not be treated since very large supercells are needed to model the structure. The target in Ref. 61 is a vortex core at which six kinds of structural domains meet. To model the complex structure with a periodic boundary condition, two pairs of vortex/antivortex cores need to be included in the calculation cell, which contains at least about 3,600 atoms. Using the MSSF method with CONQUEST, the atomic-scale structure of the vortices and their electronic structures have been investigated. TZDP PAOs were contracted to MSSFs with \(r_{\text{MS}} = 6.4\) Å in the calculations.

Before studying the complex topologically protected vortex cores, calculations of two domain walls in the \(1\times 12\times 1\) supercell containing 360 atoms (shown in Fig. 7) have been performed. The accuracy of the calculations using primitive TZDP PAOs and MSSFs was confirmed by comparing to plane-wave calculations. The phase Φ and the tilt angle Q [Fig. 7(a)] in the optimized structure are shown in Figs. 7(b) and 7(c). The profiles of Φ and Q found by TZDP PAOs and MSSFs are almost the same and are only slightly different from those found by the plane wave calculations. The average of the formation energy of the two domain walls in the two-domain model were calculated as 15.22, 14.43, and 13.72 mJ/m2 by TZDP PAOs, MSSFs and the plane waves, respectively, which are within ∼3 mJ/m2 (∼0.1 meV/Å2) difference. All these results support the robustness and accuracy of the MSSF method in this system.


Figure 7. (Color) (a) Crystal structure of ferroelectric YGaO3 (\(P6_{3}cm\) space group). The amplitude Q and phase Φ are quantified by the GaO5 tilt angle relative to the [001] direction and the GaO5 tilt direction projected onto the ab plane, respectively. (b) Φ and (c) Q in the optimized structure of 360 atoms \(1\times 12\times 1\) supercell with two-domain patterns. (Reproduced with permission from Ref. 61. © 2020 American Physical Society.)

Then the topologically protected vortex using a 3,600-atom, \(10\times 12\times 1\) supercell, have been investigated with MSSFs. The model consists of two vortex/antivortex pairs, where the domain walls and the vortices are initialized in \(P6_{3}/mmc\) symmetry, as taken by the paraelectric phase. Φ and Q in the optimized structure are visualized in Figs. 8(a) and 8(b). The optimized structure indicates that the structure around the core adopts \(P3c1\) symmetry. The DOS and the electronic density in the energy region around the conduction band minimum are illustrated in Figs. 8(c) and 8(d). In Fig. 8(d), the vortex core has larger amplitude than the other areas. According to the plane-wave calculations, the band gap of the YGaO3 unit cell is decreased from 3.19 to 2.76 eV by shifting the CBM due to the symmetry change from \(P6_{3}/mmc\) to \(P3c1\). These results suggest that the band gap of YGaO3 will be reduced by the symmetry change at the vortex core. Thus, it has been demonstrated that the large-scale DFT calculations with MSSFs help us to investigate details of complex vortex structures, leading to specific electronic structures.


Figure 8. (Color) Surface map of (a) Φ and (b) Q in the optimized structure and (c) density of states and (d) electron density around the conduction band minimum [yellow region in (c)] of 3600 atoms \(10\times 12\times 1\) supercell with two vortex/antivortex pairs. (Reproduced with permission from Ref. 61. © 2020 American Physical Society.)

PbTiO3 films on SrTiO3 substrates

We have investigated a perovskite material, PbTiO3 films on SrTiO3 substrate.62) Advanced deposition techniques allow the creation of thin film perovskite oxides and layered heterostructures, which demonstrate a wide variety of electrical polarization textures with possible applications in low dimensional functional devices. These textures arise from the interaction of different order parameters, notably anti-ferro-distortive (AFD) and ferroelectric (FE) distortions. At the surface of PbTiO3 (PTO), antiphase rotations of the TiO6 octahedra give rise to an AFD c(\(2\times 2\)) reconstruction. We used CONQUEST with MSSF to model films of PTO on SrTiO3 (STO) with a variety of polar morphologies: paraelectric; monodomain FE in-plane [both along (100) and (110) directions]; and polydomain FE films with polarization along (001). This last morphology is particularly challenging, as it requires the simulation cell to include the substrate, the in-plane width of the domain (which increases with thickness) and a doubling along the (010) direction to allow for the AFD rotations. Most simulations of thin films use a superlattice, removing the possibility to examine the important effect of the surface on the polarization textures.

We found that seven layers of STO substrate were required to avoid any influence on the PTO; we then built cells with between one and nine layers of PTO for the polar morphologies (described above/shown in Fig. 9). We used MSSFs with \(r_{\text{MS}} = 6.4\) Å (which converged all relevant parameters) and worked with the LDA, which gives excellent values of bulk polarization for PTO. In the polydomain film, domain walls lie on PbO planes.


Figure 9. (Color online) The initial supercell configurations for the \(N_{\text{z}} = 3\) films before structural relaxation. Shown here are the supercells not including AFD modes. Each configuration is, however, also treated with AFD modes following the explanation in Sect. 3.3. (a) The paraelectric supercell constrained such that spontaneous polarization cannot emerge. (b) The monodomain in-plane ferroelectric case (\(\mathbf{P}\parallel [100]\) is shown here, but we also treat \(\mathbf{P}\parallel [110]\)) constrained such that spontaneous polarization cannot develop in the out-of-plane direction. (c) The polydomain ferroelectric case with equally sized up and down domains for the ferroelectric polarization. Shown here is the \(\Lambda = 6\) case. (Used with permission from Ref. 62. © 2020 The Authors.)

We find that the polydomain arrangement is unstable for thicknesses of \(\mathrm{N}<3\), and is less stable than the monodomain (110) film with AFD distortions, though is more stable than the (100) monodomain film, up to \(\mathrm{N}=5\); after that point, polydomains are most stable. This fits well with experimental observations that polydomains start to be observed for \(\mathrm{N}\geqq 3\). The local polarization fields for two film thicknesses are shown in Fig. 10. Thick films, such as that shown in Fig. 10(a), form flux-closure domains with clear domain walls visible. For thin films, we find a polar wave morphology with small chiral cylindrical bubbles forming at the surface. As a consequence of these detailed reconstructions, we find that the surface periodicity is not c(\(2\times 2\)) as would be expected from the AFD distortions alone, but instead p(\(2\times\Lambda\)), where Λ depends on the film thickness.


Figure 10. (Color online) The local polarization vector fields in the xz plane for two film thicknesses not including AFD modes. (a) The flux-closure domains of the \(N_{\text{z}} = 9\), \(\Lambda = 12\) film. The red dashed area highlights a vortex/antivortex pair. (b) The polar wave morphology in the \(N_{\text{z}} = 3\), \(\Lambda = 6\) film. The red area indicates a cylindrical chiral bubble. (Used with permission from Ref. 62. © 2020 The Authors.)

The full exploration of these exotic polarization textures at thin film surfaces is only possible through the use of large-scale DFT, as enabled by MSSF in CONQUEST.

Hydrated DNA together with Sakurai–Sugiura method

Large-scale calculation methods are also important to investigate non-periodic materials such as glassy materials, polymers and biomaterials. For example, in this section, MSSF have been used to model a hydrated DNA system,59) shown in Fig. 11. The solvent water molecules have been treated explicitly in the calculations, therefore the system consists of 634 atoms in the DNA, 932 hydrating water molecules and 9 Mg counterions, in total 3,439 atoms. DZP PAOs with 27,883 primitive functions have been contracted to 7,447 MSSFs. Figure 11(a) compares the DOS of the hydrated DNA calculated by the primitive PAOs and the MSSFs. As discussed in Sect. 2.3.2, the MSSFs have reproduced the DOS of the primitive PAOs with high accuracy for the occupied states but not for unoccupied states. To improve the accuracy of the unoccupied states, we have introduced the SSM.59,67,68) SSM is an interior eigenproblem solver for large sparse matrices, providing the eigenvalues and eigenvectors in given energy regions with high parallel efficiency.67,68) We have performed a SCF calculation to optimize the electronic density using MSSFs and re-constructed the electronic Hamiltonian in primitive PAO basis with the optimized electronic density. Subsequently, one-shot SSM calculations have been performed for the energy region of interest. For the hydrated DNA system, we have performed the SSM calculations for the energy range [\(-1:1\)] eV with an interval of 0.027 eV (0.001 hartree). Thus, the unoccupied states of the hydrated DNA have been improved as in Fig. 11(a). Figure 11(b) shows the calculated HOMO and LUMO of the DNA system. The investigation of the KS states is important to investigate materials properties such as intramolecular electron transfers. The combination of SSM with \(O(N)\) calculations in CONQUEST is also a powerful tool to investigate the electronic structure of extremely large systems. Although the \(O(N)\) method itself does not provide information of KS eigenstates, SSM can provide the eigenstates (= KS states) for the \(O(N)\) Hamiltonian. The KS-state calculations around the Fermi level for a system with 194,573 atoms has been achieved with this combination.59)


Figure 11. (Color online) (a) DOS calculated using MSSFs (\(r_{\text{MS}} = 4.2\) Å) (blue, lower line), SSM (red, upper line), and primitive PAOs (inset) and (b) molecular orbital pictures of hydrated DNA calculated using SSM. (Reprinted with permission from Ref. 59. © 2017 American Chemical Society.)

Metallic nanoparticles

In this section, we show an example of the calculations of metallic nanoparticles with MSSFs briefly. The size-controlled metallic nanoparticles show high catalytic reactivity, and the combination of nanoparticles and substrate is one of the important factors to affect the reactivity. The interface between the nanoparticle and the substrate is a kind of hyper-ordered structure. We have investigated an Au nanoparticle with 923 atoms in octahedral (Oh) symmetry, consisting of six layers [Fig. 12(a)], using DZP PAOs. The diameter of this six-layered nanoparticle is about 3 nm, which is close to the sizes used in actual experiments.69) The nano-size calculation model enables us to investigate the site-dependence of the atomic and electronic structures of nanoparticles. Figures 13(a) and 13(b) show the intra- and inter-layer nearest neighbor atomic distances of the nanoparticle, respectively. The atomic distances of the inner layers are close to those in a bulk fcc system, while they are distributed widely in the outer layers. The wide distribution of the intra-layer distances corresponds to the site dependence, i.e., the atomic distances around the center of the faces (about 2.9 Å) are longer than those around the vertices and edges (about 2.8 Å). Figures 14(a) and 14(b) show the projected DOS (pDOS) of an Au atom in a bulk system and at the vertex of the nanoparticle. The electronic structure at the vertex of the nanoparticle is quite different from that in the bulk, where the d-band center is shifted closer to the Fermi level, which suggests high reactivity at the vertex.


Figure 12. (Color online) Optimized structures of (a) Au nanoparticle in Oh symmetry with 923 atoms and (b) Au particle on Mg(001) surface.


Figure 13. (a) Intra-layer and (b) inter-layer nearest atom distances in Au nanoparticle with 923 atoms. Black horizontal lines correspond to the Au–Au distance in a bulk fcc gold (2.95 Å). Abscissa corresponds to the indices of the layers increasing from inner to outer of the nanoparticle. The sixth layer is the surface of the nanoparticle.


Figure 14. PDOS of an Au atom (a) in bulk system and (b) at the vertex of an Au nanoparticle. Fermi levels are set to be zero (dashed line).

Not only the vertices of the nanoparticles but also the interface between the nanoparticles and the substrate have been considered as reaction active sites.3) To treat nano-scale metallic nanoparticles (i.e., not clusters) on substrates, large calculation models with several thousand of atoms are required. With MSSFs, we can treat these large models of the catalytic systems. For example, Fig. 12(b) shows the optimized structures of Au nanoparticles on the MgO(001) substrate with 2,844 atoms in total, in which we removed the bottom part of the nanoparticle as found in experimental observations.3,70) The detailed investigation of atomic and electronic structures and reactivity of metallic nanoparticles on the substrate will be provided in future studies.

4. Conclusions

We have reviewed large-scale calculation methods, especially focusing on the methods in our large-scale DFT code CONQUEST.610) To model nano-scale complex structures, we often need large calculation models with several thousand atoms or more. DFT is a powerful tool to investigate the atomic and electronic structures of materials with high accuracy, but the computational cost of conventional DFT calculation methods is quite high (cubic scaling to system sizes). Therefore, special calculation techniques to treat such large systems are required.

The MSSF method11,12) in CONQUEST makes it possible to improve both computational efficiency and accuracy. MSSFs are linear combinations of basis functions which belong not only to a target atom but also to its neighboring atoms, like local MOs. This MO-like picture of MSSFs enables us to reduce the number of the support functions to a minimal-basis size, while we can increase the number of the basis functions to improve computational accuracy, without increasing the number of the support functions. The linear-combination coefficients can be determined by using the local filter diagonalization method11,54,55) and subsequent numerical optimization.12) The investigations of accuracy for bulk Si and Al, hydrated DNA, bulk Fe and NiO demonstrate that MSSFs are applicable to varied materials such as insulating, semiconducting, metallic, and spin-polarized systems.

Examples of applications of MSSF to large systems with several thousand of atoms have been also shown. The geometry and electronic structure around complex interfaces were investigated using MSSF in these examples.5962) MSSF have been also applied to non-periodic materials such as biomolecules and metallic nanoparticle catalysts. The combination of MSSF and the Sakurai–Sugiura method,67,68) an efficient interior eigenproblem solver, enable the MSSFs to be used to investigate the excited states of large systems.59) Thus, we suggest that MSSF is now one of the most promising tools for investigation of hyper-ordered structures such as interfaces between nanoparticle catalysts and substrates.

Acknowledgment

This work is supported by the World Premier International Research Centre Initiative (WPI Initiative) on Materials Nanoarchitectonics (MANA), JSPS Grant-in-Aid for Transformative Research Areas (A) “Hyper-Ordered Structures Science” (Grant Nos. JP20H05883 and JP20H05878), JSPS Grant-in-Aid for Scientific Research (Grant No. JP18H01143), and JST PRESTO (Grant No. JPMJPR20T4). This work is also partially supported by a project, JPNP16010, commissioned by the New Energy and Industrial Technology Development Organization (NEDO). Calculations were performed on the Numerical Materials Simulator at NIMS and the supercomputer HA8000 system at Kyushu University in Japan.

The authors are grateful for computational support from the U.K. Materials and Molecular Modeling Hub, which is partially funded by EPSRC (Grant No. EP/P020194), for which access was obtained via the UKCP consortium and funded by EPSRC (Grant Ref. No. EP/P022561/1). We also acknowledge computational support from the UK national high-performance computing service, ARCHER, for which access was obtained via the UKCP consortium and funded by EPSRC (Grant Ref. No. EP/P022561/1).


References

  • 1 S. Kohara, M. Shiga, Y. Onodera, H. Masai, A. Hirata, M. Murakami, T. Morishita, K. Kimura, and K. Hayashi, Sci. Rep. 11, 22180 (2021). 10.1038/s41598-021-00965-5 CrossrefGoogle Scholar
  • 2 Y. Takano, H. X. Kondo, Y. Kanematsu, and Y. Imada, Jpn. J. Appl. Phys. 59, 010502 (2020). 10.7567/1347-4065/ab62b9 CrossrefGoogle Scholar
  • 3 H. Yoshida, Y. Kuwauchi, J. R. Jinschek, K. Sun, S. Tanaka, M. Kohyama, S. Shimada, M. Haruta, and S. Takeda, Science 335, 317 (2012). 10.1126/science.1213194 CrossrefGoogle Scholar
  • 4 M. Kawarasaki, K. Tanabe, I. Terasaki, Y. Fujii, and H. Taniguchi, Sci. Rep. 7, 5351 (2017); 10.1038/s41598-017-05651-z Crossref;, Google ScholarK. Tsutsui and Y. Morikawa, Jpn. J. Appl. Phys. 59, 010503 (2020). 10.7567/1347-4065/ab603e CrossrefGoogle Scholar
  • 5 R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, Cambridge, U.K., 2004); Crossref;, Google ScholarR. O. Jones, Rev. Mod. Phys. 87, 897 (2015). 10.1103/RevModPhys.87.897 CrossrefGoogle Scholar
  • 6 Conquest website, http://www.order-n.org (accessed February 2022). Google Scholar
  • 7 E. Hernández, M. J. Gillan, and C. M. Goringe, Phys. Rev. B 53, 7147 (1996). 10.1103/PhysRevB.53.7147 CrossrefGoogle Scholar
  • 8 D. R. Bowler, R. Choudhury, M. J. Gillan, and T. Miyazaki, Phys. Status Solidi B 243, 989 (2006). 10.1002/pssb.200541386 CrossrefGoogle Scholar
  • 9 D. R. Bowler and T. Miyazaki, J. Phys.: Condens. Matter 22, 074207 (2010). 10.1088/0953-8984/22/7/074207 CrossrefGoogle Scholar
  • 10 A. Nakata, J. S. Baker, S. Y. Mujahed, J. T. L. Poulton, S. Arapan, J. Lin, Z. Raza, S. Yadav, L. Truflandier, T. Miyazaki, and D. R. Bowler, J. Chem. Phys. 152, 164112 (2020). 10.1063/5.0005074 CrossrefGoogle Scholar
  • 11 A. Nakata, D. R. Bowler, and T. Miyazaki, J. Chem. Theory Comput. 10, 4813 (2014). 10.1021/ct5004934 CrossrefGoogle Scholar
  • 12 A. Nakata, D. R. Bowler, and T. Miyazaki, Phys. Chem. Chem. Phys. 17, 31427 (2015). 10.1039/C5CP00934K CrossrefGoogle Scholar
  • 13 M. Arita, S. Arapan, D. R. Bowler, and T. Miyazaki, J. Adv. Simulation Sci. Eng. 1, 87 (2014). 10.15748/jasse.1.87 CrossrefGoogle Scholar
  • 14 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 10.1103/PhysRev.136.B864 CrossrefGoogle Scholar
  • 15 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 10.1103/PhysRev.140.A1133 CrossrefGoogle Scholar
  • 16 D. R. Bowler and T. Miyazaki, Rep. Prog. Phys. 75, 036503 (2012). 10.1088/0034-4885/75/3/036503 CrossrefGoogle Scholar
  • 17 W. Dawson, A. Degomme, M. Stella, T. Nakajima, L. E. Ratcliff, and L. Genovese, WIREs Comput. Mol. Sci. 12, e1574 (2021); 10.1002/wcms.1574 Crossref;, Google ScholarL. E. Ratcliff, S. Mohr, G. Huhs, T. Deutsch, M. Masella, and L. Genovese, WIREs Comput. Mol. Sci. 7, e1290 (2017). 10.1002/wcms.1290 CrossrefGoogle Scholar
  • 18 W. Yang, Phys. Rev. Lett. 66, 1438 (1991). 10.1103/PhysRevLett.66.1438 CrossrefGoogle Scholar
  • 19 W. Yang, Phys. Rev. A 44, 7823 (1991). 10.1103/PhysRevA.44.7823 CrossrefGoogle Scholar
  • 20 W. Yang, J. Mol. Struct.: THEOCHEM 255, 461 (1992). 10.1016/0166-1280(92)85024-F CrossrefGoogle Scholar
  • 21 G. Galli and M. Parrinello, Phys. Rev. Lett. 69, 3547 (1992). 10.1103/PhysRevLett.69.3547 CrossrefGoogle Scholar
  • 22 E. Tsuchida, J. Phys. Soc. Jpn. 76, 034708 (2007). 10.1143/JPSJ.76.034708 LinkGoogle Scholar
  • 23 R. McWeeny, Rev. Mod. Phys. 32, 335 (1960). 10.1103/RevModPhys.32.335 CrossrefGoogle Scholar
  • 24 X. P. Li, R. W. Nunes, and D. Vanderbilt, Phys. Rev. B 47, 10891 (1993). 10.1103/PhysRevB.47.10891 CrossrefGoogle Scholar
  • 25 K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi, Chem. Phys. Lett. 313, 701 (1999). 10.1016/S0009-2614(99)00874-X CrossrefGoogle Scholar
  • 26 J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002). 10.1088/0953-8984/14/11/302 CrossrefGoogle Scholar
  • 27 A. García, N. Papior, A. Akhtar, E. Artacho, V. Blum, E. Bosoni, P. Brandimarte, M. Brandbyge, J. I. Cerdá, F. Corsetti, R. Cuadrado, V. Dikan, J. Ferrer, J. Gale, P. García-Fernández, V. M. García-Suárez, S. García, G. Huhs, S. Illera, R. Korytár, P. Koval, I. Lebedeva, L. Lin, P. López-Tarifa, S. G. Mayo, S. Mohr, P. Ordejón, A. Postnikov, Y. Pouillon, M. Pruneda, R. Robles, D. Sánchez-Portal, J. M. Soler, R. Ullah, V. W-z Yu, and J. Junquera, J. Chem. Phys. 152, 204108 (2020). 10.1063/5.0005077 CrossrefGoogle Scholar
  • 28 T. Ozaki, Phys. Rev. B 67, 155108 (2003). 10.1103/PhysRevB.67.155108 CrossrefGoogle Scholar
  • 29 S. Tsuneyuki, T. Kobori, K. Akagi, K. Sodeyama, K. Terakura, and H. Fukuyama, Chem. Phys. Lett. 476, 104 (2009). 10.1016/j.cplett.2009.05.069 CrossrefGoogle Scholar
  • 30 G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. E. Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, T. Harville, S. Irle, J. Ivanic, K. Kowalski, S. S. Leang, H. Li, W. Li, J. J. Lutz, I. Magoulas, J. Mato, V. Mironov, H. Nakata, B. Q. Pham, P. Piecuch, D. Poole, S. R. Pruitt, A. P. Rendell, L. B. Roskop, K. Ruedenberg, T. Sattasathuchana, M. W. Schmidt, J. Shen, L. Slipchenko, M. Sosonkina, V. Sundriyal, A. Tiwari, J. L. G. Vallejo, B. Westheimer, M. Wloch, P. Xu, F. Zahariev, and M. S. Gordon, J. Chem. Phys. 152, 154102 (2020). 10.1063/5.0005188 CrossrefGoogle Scholar
  • 31 S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, and K. Fukuzawa, Phys. Chem. Chem. Phys. 16, 10310 (2014). 10.1039/C4CP00316K CrossrefGoogle Scholar
  • 32 E. Tsuchida and M. Tsukada, J. Phys. Soc. Jpn. 67, 3844 (1998). 10.1143/JPSJ.67.3844 LinkGoogle Scholar
  • 33 J. C. A. Prentice, J. Aarons, J. C. Womack, A. E. A. Allen, L. Andrinopoulos, L. Anton, R. A. Bell, A. Bhandari, G. A. Bramley, R. J. Charlton, R. J. Clements, D. J. Cole, G. Constantinescu, F. Corsetti, S. M.-M. Dubois, K. K. B. Duff, J. M. Escartín, A. Greco, Q. Hill, L. P. Lee, E. Linscott, D. D. O’Regan, M. J. S. Phipps, L. E. Ratcliff, Á. Ruiz Serrano, E. W. Tait, G. Teobaldi, V. Vitale, N. Yeung, T. J. Zuehlsdorff, J. Dziedzic, P. D. Haynes, N. D. M. Hine, A. A. Mostofi, M. C. Payne, and C.-K. Skylaris, J. Chem. Phys. 152, 174111 (2020). 10.1063/5.0004445 CrossrefGoogle Scholar
  • 34 M. Ceriotti, T. D. Kühne, and M. Parrinello, J. Chem. Phys. 129, 024707 (2008). 10.1063/1.2949515 CrossrefGoogle Scholar
  • 35 L. E. Ratcliff, W. Dawson, G. Fisicaro, D. Caliste, S. Mohr, A. Degomme, B. Videau, V. Cristiglio, M. Stella, M. D’Alessandro, S. Goedecker, T. Nakajima, T. Deutsch, and L. Genovese, J. Chem. Phys. 152, 194110 (2020). 10.1063/5.0004792 CrossrefGoogle Scholar
  • 36 A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002). 10.1103/PhysRevB.66.155115 CrossrefGoogle Scholar
  • 37 E. Rudberg, E. H. Rubensson, P. Sałek, and A. Kruchinina, SoftwareX 7, 107 (2018). 10.1016/j.softx.2018.03.005 CrossrefGoogle Scholar
  • 38 T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, J. Chem. Phys. 152, 194103 (2020). 10.1063/5.0007045 CrossrefGoogle Scholar
  • 39 J. Junquera, Ó. Paz, D. Sánchez-Portal, and E. Artacho, Phys. Rev. B 64, 235111 (2001). 10.1103/PhysRevB.64.235111 CrossrefGoogle Scholar
  • 40 E. Anglada, J. M. Soler, J. Junquera, and E. Artacho, Phys. Rev. B 66, 205101 (2002). 10.1103/PhysRevB.66.205101 CrossrefGoogle Scholar
  • 41 T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004). 10.1103/PhysRevB.69.195113 CrossrefGoogle Scholar
  • 42 A. S. Torralba, M. Todorović, V. Brázdová, R. Choudhury, T. Miyazaki, M. J. Gillan, and D. R. Bowler, J. Phys.: Condens. Matter 20, 294206 (2008). 10.1088/0953-8984/20/29/294206 CrossrefGoogle Scholar
  • 43 D. R. Bowler, J. S. Baker, J. T. L. Poulton, S. Y. Mujahed, J. Lin, S. Yadav, Z. Raza, and T. Miyazaki, Jpn. J. Appl. Phys. 58, 100503 (2019). 10.7567/1347-4065/ab45af CrossrefGoogle Scholar
  • 44 E. Hernández, M. J. Gillan, and C. M. Goringe, Phys. Rev. B 55, 13485 (1997). 10.1103/PhysRevB.55.13485 CrossrefGoogle Scholar
  • 45 A. A. Mostofi, P. D. Haynes, C.-K. Skylaris, and M. C. Payne, J. Chem. Phys. 119, 8842 (2003). 10.1063/1.1613633 CrossrefGoogle Scholar
  • 46 S. Goedecker, Wavelets and Their Application: For The Solution of Partial Differential Equations in Physics (Presses Polytechniques et Universitaires Romandes, Switzerland, 1998). Google Scholar
  • 47 J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, Comput. Phys. Commun. 167, 103 (2005). 10.1016/j.cpc.2004.12.014 CrossrefGoogle Scholar
  • 48 J. Iwata, D. Takahashi, A. Oshiyama, T. Boku, K. Shiraishi, S. Okada, and K. Yabana, J. Comput. Phys. 229, 2339 (2010). 10.1016/j.jcp.2009.11.038 CrossrefGoogle Scholar
  • 49 J. R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72, 1240 (1994). 10.1103/PhysRevLett.72.1240 CrossrefGoogle Scholar
  • 50 D. R. Hamann, Phys. Rev. B 88, 085117 (2013). 10.1103/PhysRevB.88.085117 CrossrefGoogle Scholar
  • 51 D. Vanderbilt, Phys. Rev. B 32, 8412 (1985). 10.1103/PhysRevB.32.8412 CrossrefGoogle Scholar
  • 52 ATOM code for the generation of norm-conserving pseudopotentials. The version maintained by the SIESTA project, http://icmab.es/siesta/Pseudopotentials/index.html (accessed February 2022). Google Scholar
  • 53 E. Hernández and M. J. Gillan, Phys. Rev. B 51, 10157 (1995). 10.1103/PhysRevB.51.10157 CrossrefGoogle Scholar
  • 54 M. J. Rayson and P. R. Briddon, Phys. Rev. B 80, 205104 (2009). 10.1103/PhysRevB.80.205104 CrossrefGoogle Scholar
  • 55 M. J. Rayson, Comput. Phys. Commun. 181, 1051 (2010). 10.1016/j.cpc.2010.02.012 CrossrefGoogle Scholar
  • 56 T. Miyazaki, D. R. Bowler, R. Choudhury, and M. J. Gillan, J. Chem. Phys. 121, 6186 (2004). 10.1063/1.1787832 CrossrefGoogle Scholar
  • 57 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 10.1103/PhysRevB.23.5048 CrossrefGoogle Scholar
  • 58 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 10.1103/PhysRevLett.77.3865 CrossrefGoogle Scholar
  • 59 A. Nakata, Y. Futamura, T. Sakurai, D. R. Bowler, and T. Miyazaki, J. Chem. Theory Comput. 13, 4146 (2017). 10.1021/acs.jctc.7b00385 CrossrefGoogle Scholar
  • 60 C. Romero-Muñiz, A. Nakata, P. Pou, D. R. Bowler, T. Miyazaki, and R. Pérez, J. Phys.: Condens. Matter 30, 505901 (2018). 10.1088/1361-648X/aaec4c CrossrefGoogle Scholar
  • 61 D. R. Småbråten, A. Nakata, D. Meier, T. Miyazaki, and S. M. Selbach, Phys. Rev. B 102, 144103 (2020). 10.1103/PhysRevB.102.144103 CrossrefGoogle Scholar
  • 62 J. S. Baker and D. R. Bowler, Adv. Theory Simul. 3, 2000154 (2020). 10.1002/adts.202000154 CrossrefGoogle Scholar
  • 63 P. E. Blöchl, Phys. Rev. B 50, 17953 (1994). 10.1103/PhysRevB.50.17953 CrossrefGoogle Scholar
  • 64 G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996). 10.1103/PhysRevB.54.11169 CrossrefGoogle Scholar
  • 65 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 10.1103/PhysRevB.59.1758 CrossrefGoogle Scholar
  • 66 S. Grimme, J. Comput. Chem. 27, 1787 (2006). 10.1002/jcc.20495 CrossrefGoogle Scholar
  • 67 T. Sakurai and H. Sugiura, J. Comput. Appl. Math. 159, 119 (2003). 10.1016/S0377-0427(03)00565-X CrossrefGoogle Scholar
  • 68 T. Sakurai, Y. Futamura, and H. Tadano, J. Algo. & Comp. Tech. 7, 249 (2013). 10.1260/1748-3018.7.3.249 CrossrefGoogle Scholar
  • 69 T. Mitsudome, A. Noujima, T. Mizugaki, K. Jitsukawa, and K. Kaneda, Chem. Commun. 2009, 5302 (2009); 10.1039/b910208f Crossref;, Google ScholarT. Urayama, T. Mitsudome, Z. Maeno, T. Mizugaki, K. Jitsukawa, and K. Kaneda, Chem. Lett. 44, 1062 (2015). 10.1246/cl.150379 CrossrefGoogle Scholar
  • 70 S. Takeda, Y. Kuwauchi, and H. Yoshida, Ultramicroscopy 151, 178 (2015). 10.1016/j.ultramic.2014.11.017 CrossrefGoogle Scholar

Author Biographies


Ayako Nakata (ORCiD: 0000-0002-3311-6283) is a Principal Researcher at NIMS, a researcher in JST-PREST, and an adjunctive associate professor of cooperative graduate school at University of Tsukuba. She obtained her Ph.D. degree from Waseda university in 2007. She has been a member of the CONQUEST developer team and has mainly worked on the development and applications of the multi-site method. She has also interested in theoretical investigation of catalysts (especially metallic nanoparticles) and excited-state calculations.

David R. Bowler (ORCiD: 0000-0001-7853-1520) is Professor of Physics at UCL, and a PI in both the London Centre for Nanotechnology and the World Premier Research Institute for Materials Nanoarchitectonics (MANA) in the National Institute for Materials Science (NIMS), Japan. He obtained his Ph.D. degree from Oxford University in 1997. He has driven the development of the massively-parallel linear scaling density functional theory code, CONQUEST, and collaborates extensively with experimental groups on the growth and properties of nanostructures on semiconductor surfaces.

Tsuyoshi Miyazaki (ORCiD: 0000-0003-3534-4404) obtained a Ph.D. from University of Tokyo in 1995. He is now a MANA Principal Investigator at the International Center for Nanoarchitectonics (MANA), National Research Institute for Materials Science (NIMS), and a group leader of First-principles Simulation Group in Nano-Theory Field of NIMS-MANA. He has performed theoretical studies of various materials, such as organic conductors, semiconductor surfaces, nano-structured materials, biological system and so on. He is the co-leader of the CONQUEST project, a visiting Professor of School of Engineering at Nagoya University, and a Professor of School of Integrative and Global Majors at University of Tsukuba.