Adaptively Compressed Exchange Operator for Large-Scale Hybrid Density Functional Calculations with Applications to the Adsorption of Water on Silicene

Efficient exact exchange using Wannier functions and other related developments in planewave-pseudopotential implementation of RT-TDDFT

The plane-wave pseudopotential (PW-PP) formalism is widely used for the first-principles electronic structure calculation of extended periodic systems. The PW-PP approach has also been adapted for real-time time-dependent density functional theory (RT-TDDFT) to investigate time-dependent electronic dynamical phenomena. In this work, we detail recent advances in the PW-PP formalism for RT-TDDFT, particularly how maximally localized Wannier functions (MLWFs) are used to accelerate simulations using the exact exchange. We also discuss several related developments, including an anti-Hermitian correction for the time-dependent MLWFs (TD-MLWFs) when a time-dependent electric field is applied, the refinement procedure for TD-MLWFs, comparison of the velocity and length gauge approaches for applying an electric field, and elimination of long-range electrostatic interaction, as well as usage of a complex absorbing potential for modeling isolated systems when using the PW-PP formalism.

Complex-Valued K-Means Clustering of Interpolative Separable Density Fitting Algorithm for Large-Scale Hybrid Functional Enabled Ab Initio Molecular Dynamics Simulations within Plane Waves

K-means clustering, as a classic unsupervised machine learning algorithm, is the key step to select the interpolation sampling points in interpolative separable density fitting (ISDF) decomposition for hybrid functional electronic structure calculations. Real-valued K-means clustering for accelerating the ISDF decomposition has been demonstrated for large-scale hybrid functional enabled ab initio molecular dynamics (hybrid AIMD) simulations within plane-wave basis sets where the Kohn-Sham orbitals are real-valued. However, it is unclear whether such K-means clustering works for complex-valued Kohn-Sham orbitals. Here, we propose an improved weight function defined as the sum of the square modulus of complex-valued Kohn-Sham orbitals in K-means clustering for hybrid AIMD simulations. Numerical results demonstrate that the K-means algorithm with a new weight function yields smoother and more delocalized interpolation sampling points, resulting in smoother energy potential, smaller energy drift, and longer time steps for hybrid AIMD simulations compared to the previous weight function used in the real-valued K-means algorithm. In particular, we find that this improved algorithm can obtain more accurate oxygen-oxygen radial distribution functions in liquid water molecules and a more accurate power spectrum in crystal silicon dioxide compared to the previous K-means algorithm. Finally, we describe a massively parallel implementation of this ISDF decomposition to accelerate large-scale complex-valued hybrid AIMD simulations containing thousands of atoms (2,744 atoms), which can scale up to 5,504 CPU cores on modern supercomputers.

Carbon Kagome nanotubes-quasi-one-dimensional nanostructures with flat bands

In recent years, a number of bulk materials and heterostructures have been explored due their connections with exotic materials phenomena emanating from flat band physics and strong electronic correlation. The possibility of realizing such fascinating material properties in simple realistic nanostructures is particularly exciting, especially as the investigation of exotic states of electronic matter in wire-like geometries is relatively unexplored in the literature. Motivated by these considerations, we introduce in this work carbon Kagome nanotubes (CKNTs)-a new allotrope of carbon formed by rolling up Kagome graphene, and investigate this material using specialized first principles calculations. We identify two principal varieties of CKNTs-armchair and zigzag, and find both varieties to be stable at room temperature, based on ab initio molecular dynamics simulations. CKNTs are metallic and feature dispersionless states (i.e., flat bands) near the Fermi level throughout their Brillouin zone, along with an associated singular peak in the electronic density of states. We calculate the mechanical and electronic response of CKNTs to torsional and axial strains, and show that CKNTs appear to be more mechanically compliant than conventional carbon nanotubes (CNTs). Additionally, we find that the electronic properties of CKNTs undergo significant electronic transitions-with emergent partial flat bands and tilted Dirac points-when twisted. We develop a relatively simple tight-binding model that can explain many of these electronic features. We also discuss possible routes for the synthesis of CKNTs. Overall, CKNTs appear to be unique and striking examples of realistic elemental quasi-one-dimensional materials that may display fascinating material properties due to strong electronic correlation. Distorted CKNTs may provide an interesting nanomaterial platform where flat band physics and chirality induced anomalous transport effects may be studied together.

Speeding-up Hybrid Functional-Based Ab Initio Molecular Dynamics Using Multiple Time-stepping and Resonance-Free Thermostat

Ab initio molecular dynamics (AIMD) based on density functional theory (DFT) has become a workhorse for studying the structure, dynamics, and reactions in condensed matter systems. Currently, AIMD simulations are primarily carried out at the level of generalized gradient approximation (GGA), which is at the second rung of DFT functionals in terms of accuracy. Hybrid DFT functionals, which form the fourth rung in the accuracy ladder, are not commonly used in AIMD simulations as the computational cost involved is 100 times or higher. To facilitate AIMD simulations with hybrid functionals, we propose here an approach using multiple time stepping with adaptively compressed exchange operator and resonance-free thermostat, that could speed up the calculations by ∼30 times or more for systems with a few hundred of atoms. We demonstrate that by achieving this significant speed up and making the compute time of hybrid functional-based AIMD simulations at par with that of GGA functionals, we are able to study several complex condensed matter systems and model chemical reactions in solution with hybrid functionals that were earlier unthinkable to be performed.

Efficient implementation of analytical gradients for periodic hybrid functional calculations within fitted numerical atomic orbitals from NAO2GTO

The NAO2GTO scheme provides an efficient way to evaluate the electron repulsion integrals (ERIs) over numerical atomic orbitals (NAOs) with auxiliary Gaussian-type orbitals (GTOs). However, the NAO2GTO fitting will significantly impact the accuracy and convergence of hybrid functional calculations. To address this issue, here we propose to use the fitted orbitals as a new numerical basis to properly handle the mismatch between NAOs and fitted GTOs. We present an efficient and linear-scaling implementation of analytical gradients of Hartree-Fock exchange (HFX) energy for periodic HSE06 calculations with fitted NAOs in the HONPAS package. In our implementation, the ERIs and their derivatives for HFX matrix and forces are evaluated analytically with the auxiliary GTOs, while other terms are calculated using numerically discretized GTOs. Several integral screening techniques are employed to reduce the number of required ERI derivatives. We benchmark the accuracy and efficiency of our implementation and demonstrate that our results of lattice constants, bulk moduli, and band gaps of several typical semiconductors are in good agreement with the experimental values. We also show that the calculation of HFX forces based on a master-worker dynamic parallel scheme has a very high efficiency and scales linearly with respect to system size. Finally, we study the geometry optimization and polaron formation due to an excess electron in rutile TiO2 by means of HSE06 calculations to further validate the applicability of our implementation.

Exact exchange with range-separated algorithm for thermodynamic limit of periodic Hartree-Fock theory

The expensive cost of computing exact exchange in periodic systems limits the application range of density functional theory with hybrid functionals. To reduce the computational cost of exact change, we present a range-separated algorithm to compute electron repulsion integrals for Gaussian-type crystal basis. The algorithm splits the full-range Coulomb interactions into short-range and long-range parts, which are, respectively, computed in real and reciprocal space. This approach significantly reduces the overall computational cost, as integrals can be efficiently computed in both regions. The algorithm can efficiently handle large numbers of k points with limited central processing unit (CPU) and memory resources. As a demonstration, we performed an all-electron k-point Hartree-Fock calculation for LiH crystal with one million Gaussian basis functions, which was completed on a desktop computer in 1400 CPU hours.

Interpolative Separable Density Fitting for Accelerating Two-Electron Integrals: A Theoretical Perspective

Low-rank approximations have long been considered an efficient way to accelerate electronic structure calculations associated with the evaluation of electron repulsion integrals (ERIs). As an accurate and efficient algorithm for compressing the ERI tensor, the interpolative separable density fitting (ISDF) decomposition has recently attracted great attention in this context. In this perspective, we introduce the ISDF decomposition from the theoretical aspects and technique details. The ISDF decomposition can construct a fully separable low-rank approximation (tensor hypercontraction factorization) of ERIs in real space with a cubic cost, offering great flexibility for accelerating high-scaling electronic structure calculations. We review the typical applications of ISDF in hybrid functionals, time-dependent density functional theory, and GW approximation. Finally, we discuss the promising directions for future development of ISDF.

A Chirality-Based Quantum Leap

There is increasing interest in the study of chiral degrees of freedom occurring in matter and in electromagnetic fields. Opportunities in quantum sciences will likely exploit two main areas that are the focus of this Review: (1) recent observations of the chiral-induced spin selectivity (CISS) effect in chiral molecules and engineered nanomaterials and (2) rapidly evolving nanophotonic strategies designed to amplify chiral light-matter interactions. On the one hand, the CISS effect underpins the observation that charge transport through nanoscopic chiral structures favors a particular electronic spin orientation, resulting in large room-temperature spin polarizations. Observations of the CISS effect suggest opportunities for spin control and for the design and fabrication of room-temperature quantum devices from the bottom up, with atomic-scale precision and molecular modularity. On the other hand, chiral-optical effects that depend on both spin- and orbital-angular momentum of photons could offer key advantages in all-optical and quantum information technologies. In particular, amplification of these chiral light-matter interactions using rationally designed plasmonic and dielectric nanomaterials provide approaches to manipulate light intensity, polarization, and phase in confined nanoscale geometries. Any technology that relies on optimal charge transport, or optical control and readout, including quantum devices for logic, sensing, and storage, may benefit from chiral quantum properties. These properties can be theoretically and experimentally investigated from a quantum information perspective, which has not yet been fully developed. There are uncharted implications for the quantum sciences once chiral couplings can be engineered to control the storage, transduction, and manipulation of quantum information. This forward-looking Review provides a survey of the experimental and theoretical fundamentals of chiral-influenced quantum effects and presents a vision for their possible future roles in enabling room-temperature quantum technologies.

CIDER: An Expressive, Nonlocal Feature Set for Machine Learning Density Functionals with Exact Constraints

Machine learning (ML) has recently gained attention as a means to develop more accurate exchange-correlation (XC) functionals for density functional theory, but functionals developed thus far need to be improved on several metrics, including accuracy, numerical stability, and transferability across chemical space. In this work, we introduce a set of nonlocal features of the density called the CIDER formalism, which we use to train a Gaussian process model for the exchange energy that obeys the critical uniform scaling rule for exchange. The resulting CIDER exchange functional is significantly more accurate than any semilocal functional tested here, and it has good transferability across main-group molecules. This work therefore serves as an initial step toward more accurate exchange functionals, and it also introduces useful techniques for developing robust, physics-informed XC models via ML.

Improving the scaling and performance of multiple time stepping-based molecular dynamics with hybrid density functionals

Density functionals at the level of the generalized gradient approximation (GGA) and a plane-wave basis set are widely used today to perform ab initio molecular dynamics (AIMD) simulations. Going up in the ladder of accuracy of density functionals from GGA (second rung) to hybrid density functionals (fourth rung) is much desired pertaining to the accuracy of the latter in describing structure, dynamics, and energetics of molecular and condensed matter systems. On the other hand, hybrid density functional based AIMD simulations are about two orders of magnitude slower than GGA based AIMD for systems containing ~100 atoms using ~100 compute cores. Two methods, namely MTACE and s-MTACE, based on a multiple time step integrator and adaptively compressed exchange operator formalism are able to provide a speed-up of about 7-9 in performing hybrid density functional based AIMD. In this work, we report an implementation of these methods using a task-group based parallelization within the CPMD program package, with the intention to take advantage of the large number of compute cores available on modern high-performance computing platforms. We present here the boost in performance achieved through this algorithm. This work also identifies the computational bottleneck in the s-MTACE method and proposes a way to overcome it.

Low-Rank Approximations Accelerated Plane-Wave Hybrid Functional Calculations with k-Point Sampling

The low-rank approximations of the adaptively compressed exchange (ACE) operator and interpolative separable density fitting (ISDF) algorithms significantly reduce the computational cost and memory usage of hybrid functional calculations in real space, but the lack of k-point sampling hinders their implementation in reciprocal space for periodic systems with the plane-wave basis set. Here, we combine the ACE operator and ISDF decomposition into a new ACE-ISDF algorithm for periodic systems in reciprocal space with k-point sampling. On the basis of the ACE-ISDF algorithm with the improved reciprocal space ACE operator and k-point Fourier convolution, the time complexity of the hybrid functional calculation is reduced from O(Ne4Nk2) to O(Ne3Nklog(Nk)) (N_e and N_k are the number of electrons and k-points, respectively) with a much smaller prefactor and much lower memory consumption compared to the standard method for periodic systems with a plane-wave basis set.

Achieving an Order of Magnitude Speedup in Hybrid-Functional- and Plane-Wave-Based Ab Initio Molecular Dynamics: Applications to Proton-Transfer Reactions in Enzymes and in Solution

Ab initio molecular dynamics (MD) with hybrid density functionals and a plane wave basis is computationally expensive due to the high computational cost of exact exchange energy evaluation. Recently, we proposed a strategy to combine adaptively compressed exchange (ACE) operator formulation and a multiple time step integration scheme to reduce the computational cost significantly [J. Chem. Phys. 2019, 151, 151102 ]. However, it was found that the construction of the ACE operator, which has to be done at least once in every MD time step, is computationally expensive. In the present work, systematic improvements are introduced to further speed up by employing localized orbitals for the construction of the ACE operator. By this, we could achieve a computational speedup of an order of magnitude for a periodic system containing 32 water molecules. Benchmark calculations were carried out to show the accuracy and efficiency of the method in predicting the structural and dynamical properties of bulk water. To demonstrate the applicability, computationally intensive free-energy computations at the level of hybrid density functional theory were performed to investigate (a) methyl formate hydrolysis reaction in neutral aqueous media and (b) proton-transfer reaction within the active-site residues of the class C β-lactamase enzyme.

High performance computing of DGDFT for tens of thousands of atoms using millions of cores on Sunway TaihuLight

High performance computing (HPC) is a powerful tool to accelerate the Kohn-Sham density functional theory (KS-DFT) calculations on modern heterogeneous supercomputers. Here, we describe a massively parallel implementation of discontinuous Galerkin density functional theory (DGDFT) method on the Sunway TaihuLight supercomputer. The DGDFT method uses the adaptive local basis (ALB) functions generated on-the-fly during the self-consistent field (SCF) iteration to solve the KS equations with high precision comparable to plane-wave basis set. In particular, the DGDFT method adopts a two-level parallelization strategy that deals with various types of data distribution, task scheduling, and data communication schemes, and combines with the master-slave multi-thread heterogeneous parallelism of SW26010 processor, resulting in large-scale HPC KS-DFT calculations on the Sunway TaihuLight supercomputer. We show that the DGDFT method can scale up to 8,519,680 processing cores (131,072 core groups) on the Sunway TaihuLight supercomputer for studying the electronic structures of two-dimensional (2D) metallic graphene systems that contain tens of thousands of carbon atoms.

Efficient computation of free energy surfaces of chemical reactions using ab initio molecular dynamics with hybrid functionals and plane waves

Ab initio molecular dynamics (AIMD) simulations employing density functional theory (DFT) and plane waves are routinely carried out using density functionals at the level of generalized gradient approximation (GGA). AIMD simulations employing hybrid density functionals are of great interest as it offers a more accurate description of structural and dynamic properties than the GGA functionals. However, the computational cost for carrying out calculations using hybrid functionals and plane wave basis set is at least two orders of magnitude higher than that using GGA functionals. Recently, we proposed a strategy that combined the adaptively compressed exchange operator formulation and the multiple time step integration scheme to reduce the computational cost by an order of magnitude [J. Chem. Phys. 151, 151102 (2019)]. In this work, we demonstrate the application of this method to study chemical reactions, in particular, formamide hydrolysis in an alkaline aqueous medium. By actuating our implementation with the well-sliced metadynamics scheme, we can compute the two-dimensional free energy surface of this reaction at the level of hybrid-DFT. This work also investigates the accuracy of the PBE0 (hybrid) and the PBE (GGA) functionals in predicting the free energetics of this chemical reaction.

Large-Scale Parallel Implementation of Hartree-Fock Exchange Energy on Real-Space Grids Using 3D-Parallel Fast Fourier Transform

Two types of implementation of the Hartree-Fock (HF) exchange energy were developed with the real-space grid approach for the purpose of achieving high efficiency in the parallel execution of the hybrid exchange functional in the density functional theory. First, a parallel implementation of the three-dimensional fast Fourier transform (FFT), referred to as PFFT, was adapted to solve the Poisson equations for the electrostatic potentials of the densities of the orbital pairs. In the other approach, the Poisson equations were solved iteratively through the conjugate gradient (CG) procedures where the operation of Laplacian was parallelized by the domain decomposition scheme. Comparison of the parallel performances for the exchange energy calculation was made between these two approaches, and it was revealed that the calculation with the FFT method is faster than that with CG. The method with FFT is more advantageous than CG because a larger bandwidth can be made available in the collective message passing interface communication associated with the parallel execution of FFT. We also implemented the projection operator to circumvent the laborious calculation of the exchange energy at every self-consistent field step, which made a significant contribution to expedite the convergence. To assess the accuracy of our implementation, the association energies of a hydrated ion were computed, which showed excellent agreement with those given by the Gaussian 09 program employing sophisticated basis sets.

Porous Fe@C Composites Derived from Silkworm Excrement for Effective Separation of Anisole Compounds

Silkworm excrement is a very useful biomass waste, composed of layer-structured fats and proteins, which are great precursors for carbon composite materials. In this work, new porous composites derived from silkworm excrement were prepared for selective separation of flavor 4-methylanisole from the binary 4-methylanisole/4-anisaldehyde mixture. In particular, the silkworm excrement, possessing a unique nanosheet structure, is converted into a graphite-like carbon by a simple calcination strategy followed by a metal-ion-doping procedure. This Fe@C composite exhibits a special nano-spongy morphology, anchoring Fe₃C/Fe₅C₂ on the carbon nanosheets. Density functional theory simulations showed that 4-methylanisole presents a stronger π-π interaction and attraction forces with sp² carbon nanosheets in Fe@C composites than 4-anisaldehyde. The selective adsorption experiments further confirmed that the Fe@C composites exhibited a 4-methylanisole capacity of 7.3 mmol/g at 298 K and the highest selectivity of 17 for an equimolar 4-methylanisole/4-anisaldehyde mixture among the examined adsorbents including MOFs and commercial activated carbon materials, which demonstrates the potential of this low-cost and eco-friendly porous carbon material as a promising sustainable adsorbent.

Study of Li Adsorption on Graphdiyne Using Hybrid DFT Calculations

Promising applications of graphdiyne have often been initiated by theoretical predictions especially using DFT known as the most powerful first-principles electronic structure calculation method. However, there is no systematic study on the reliability of DFT for the prediction of the electronic properties of the graphdiyne. Here, we performed a study of Li adsorption on the graphdiyne using hybrid DFT with LC-ωPBE and compared the results with those of PBE, because accurate prediction of the Li adsorption is important for performance as a Li storage that was first theoretically suggested and then experimentally realized. Our results show that PBE overestimates the adsorption energy inside a pore and the barrier height at the transition state of in-plane diffusion compared to the those of LC-ωPBE. In particular, LC-ωPBE predicted almost barrier-less in-plane diffusion of Li on the graphdiyne because of the presence of both in-plane and out-of-plane π orbitals. Also, LC-ωPBE favors a high spin state due to the exact exchange energy when several Li atoms are adsorbed on the graphdiyne, whereas PBE favors a low spin state. Thus, the use of the hybrid DFT is critical for reliable predictions on the electronic properties of the graphdiyne.

Interpolative Separable Density Fitting Decomposition for Accelerating Hybrid Density Functional Calculations with Applications to Defects in Silicon

We present a new efficient way to perform hybrid density functional theory (DFT)-based electronic structure calculations. The new method uses an interpolative separable density fitting (ISDF) procedure to construct a set of numerical auxiliary basis vectors and a compact approximation of the matrix consisting of products of occupied orbitals represented in a large basis set such as the planewave basis. Such an approximation allows us to reduce the number of Poisson solves from [Formula: see text] to [Formula: see text] when we apply the exchange operator to occupied orbitals in an iterative method for solving the Kohn-Sham equations, where N_e is the number of electrons in the system to be studied. We show that the ISDF procedure can be carried out in [Formula: see text] operations, with a much smaller preconstant compared to methods used in existing approaches. When combined with the recently developed adaptively compressed exchange (ACE) operator formalism, which reduces the number of times the exchange operator needs to be updated, the resulting ACE-ISDF method significantly reduces the computational cost associated with the exchange operator by nearly 2 orders of magnitude compared to existing approaches for a large silicon system with 1000 atoms. We demonstrate that the ACE-ISDF method can produce accurate energies and forces for insulating and metallic systems and that it is possible to obtain converged hybrid functional calculation results for a 1000-atom bulk silicon within 10 min on 2000 computational cores. We also show that ACE-ISDF can scale to 8192 computational cores for a 4096-atom bulk silicon system. We use the ACE-ISDF method to geometrically optimize a 1000-atom silicon system with a vacancy defect using the HSE06 functional and computes its electronic structure. We find that that the computed energy gap from the HSE06 functional is much closer to the experimental value compared to that produced by semilocal functionals in the DFT calculations.

Projected Commutator DIIS Method for Accelerating Hybrid Functional Electronic Structure Calculations

The commutator direct inversion of the iterative subspace (commutator DIIS or C-DIIS) method developed by Pulay is an efficient and the most widely used scheme in quantum chemistry to accelerate the convergence of self-consistent field (SCF) iterations in Hartree-Fock theory and Kohn-Sham density functional theory. The C-DIIS method requires the explicit storage of the density matrix, the Fock matrix, and the commutator matrix. Hence, the method can only be used for systems with a relatively small basis set, such as the Gaussian basis set. We develop a new method that enables the C-DIIS method to be efficiently employed in electronic structure calculations with a large basis set such as planewaves for the first time. The key ingredient is the projection of both the density matrix and the commutator matrix to an auxiliary matrix called the gauge-fixing matrix. The resulting projected commutator-DIIS method (PC-DIIS) only operates on matrices of the same dimension as that consists of Kohn-Sham orbitals. The cost of the method is comparable to that of standard charge mixing schemes used in large basis set calculations. The PC-DIIS method is gauge-invariant, which guarantees that its performance is invariant with respect to any unitary transformation of the Kohn-Sham orbitals. We demonstrate that the PC-DIIS method can be viewed as an extension of an iterative eigensolver for nonlinear problems. We use the PC-DIIS method for accelerating Kohn-Sham density functional theory calculations with hybrid exchange-correlation functionals, and demonstrate its superior performance compared to the commonly used nested two-level SCF iteration procedure. Furthermore, we demonstrate that in the context of ab initio molecular dynamics (MD) simulation with hybrid functionals one can extrapolate the gauge-fixing matrix to achieve the goal of extrapolating the entire density matrix implicitly along the MD trajectory. Numerical results indicate that the new method significantly reduces the number of SCF iterations per MD step, compared to the commonly used strategy of extrapolating the electron density.