Orbital formulation for electronic-structure calculations with linear system-size scaling

Modular implementation of the linear- and cubic-scaling orbital minimization methods in electronic structure codes using atomic orbitals

We present a code modularization approach to design efficient and massively parallel cubic- and linear-scaling solvers for electronic structure calculations using atomic orbitals. The modular implementation of the orbital minimization method, in which linear algebra and parallelization issues are handled via external libraries, is demonstrated in the SIESTA code. The distributed block compressed sparse row (DBCSR) and scalable linear algebra package (ScaLAPACK) libraries are used for algebraic operations with sparse and dense matrices, respectively. The MatrixSwitch and libOMM libraries, recently developed within the Electronic Structure Library, facilitate switching between different matrix formats and implement the energy minimization. We show results comparing the performance of several cubic-scaling algorithms, and also demonstrate the parallel performance of the linear-scaling solvers, and their supremacy over the cubic-scaling solvers for insulating systems with sizes of several hundreds of atoms.

Extending multi-layer energy-based fragment method for excited-state calculations of large covalently bonded fragment systems

Recently, we developed a low-scaling Multi-Layer Energy-Based Fragment (MLEBF) method for accurate excited-state calculations and nonadiabatic dynamics simulations of nonbonded fragment systems. In this work, we extend the MLEBF method to treat covalently bonded fragment ones. The main idea is cutting a target system into many fragments according to chemical properties. Fragments with dangling bonds are first saturated by chemical groups; then, saturated fragments, together with the original fragments without dangling bonds, are grouped into different layers. The accurate total energy expression is formulated with the many-body energy expansion theory, in combination with the inclusion-exclusion principle that is used to delete the contribution of chemical groups introduced to saturate dangling bonds. Specifically, in a two-layer MLEBF model, the photochemically active and inert layers are calculated with high-level and efficient electronic structure methods, respectively. Intralayer and interlayer energies can be truncated at the two- or three-body interaction level. Subsequently, through several systems, including neutral and charged covalently bonded fragment systems, we demonstrate that MLEBF can provide accurate ground- and excited-state energies and gradients. Finally, we realize the structure, conical intersection, and path optimizations by combining our MLEBF program with commercial and free packages, e.g., ASE and SciPy. These developments make MLEBF a practical and reliable tool for studying complex photochemical and photophysical processes of large nonbonded and bonded fragment systems.

Structure optimization with stochastic density functional theory

Linear-scaling techniques for Kohn-Sham density functional theory are essential to describe the ground state properties of extended systems. Still, these techniques often rely on the localization of the density matrix or accurate embedding approaches, limiting their applicability. In contrast, stochastic density functional theory (sDFT) achieves linear- and sub-linear scaling by statistically sampling the ground state density without relying on embedding or imposing localization. In return, ground state observables, such as the forces on the nuclei, fluctuate in sDFT, making optimizing the nuclear structure a highly non-trivial problem. In this work, we combine the most recent noise-reduction schemes for sDFT with stochastic optimization algorithms to perform structure optimization within sDFT. We compare the performance of the stochastic gradient descent approach and its variations (stochastic gradient descent with momentum) with stochastic optimization techniques that rely on the Hessian, such as the stochastic Broyden-Fletcher-Goldfarb-Shanno algorithm. We further provide a detailed assessment of the computational efficiency and its dependence on the optimization parameters of each method for determining the ground state structure of bulk silicon with varying supercell dimensions.

Quantum Orbital Minimization Method for Excited States Calculation on a Quantum Computer

We propose a quantum-classical hybrid variational algorithm, the quantum orbital minimization method (qOMM), for obtaining the ground state and low-lying excited states of a Hermitian operator. Given parametrized ansatz circuits representing eigenstates, qOMM implements quantum circuits to represent the objective function in the orbital minimization method and adopts a classical optimizer to minimize the objective function with respect to the parameters in ansatz circuits. The objective function has an orthogonality constraint implicitly embedded, which allows qOMM to apply a different ansatz circuit to each input reference state. We carry out numerical simulations that seek to find excited states of H₂, LiH, and a toy model consisting of four hydrogen atoms arranged in a square lattice in the STO-3G basis with UCCSD ansatz circuits. Comparing the numerical results with existing excited states methods, qOMM is less prone to getting stuck in local minima and can achieve convergence with more shallow ansatz circuits.

Describing Chemical Reactivity with Frontier Molecular Orbitalets

Locality in physical space is critical in understanding chemical reactivity in the analysis of various phenomena and processes in chemistry, biology, and materials science, as exemplified in the concepts of reactive functional groups and active sites. Frontier molecular orbitals (FMOs) pinpoint the locality of chemical bonds that are chemically reactive because of the associated orbital energies and thus have achieved great success in describing chemical reactivity, mainly for small systems. For large systems, however, the delocalization nature of canonical molecular orbitals makes it difficult for FMOs to highlight the locality of the chemical reactivity. To obtain localized molecular orbitals that also reflect the frontier nature of the chemical processes, we develop the concept of frontier molecular orbitalets (FMOLs) for describing the reactivity of large systems. The concept of orbitalets was developed recently in the localized orbital scaling correction method, which aims for eliminating the delocalization error in common density functional approximations. Orbitalets are localized in both physical and energy spaces and thus contain both orbital locality and energy information. The FMOLs are thus the orbitalets with energies highest among occupied orbitalets and lowest among unoccupied ones. The applications of FMOLs to hexadeca-1,3,5,7,9,11,13,15-octaene in its equilibrium geometry, inter- and intra-molecular charge-transfer systems, and two transition states of a bifurcating reaction demonstrate that FMOLs can connect quantum mechanical treatments of chemical systems and chemical reactivities by locating the reactive region of large chemical systems. Therefore, FMOLs extend the role of FMOs for small systems and describe the chemical reactivity of large systems with energy and locality insight, with potentially broad applications.

Regularized Localized Molecular Orbitals in a Divide-and-Conquer Approach for Linear Scaling Calculations

Non-orthogonal localized molecular orbitals (NOLMOs) have been employed as building blocks for the divide-and-conquer (DC) linear scaling method. The NOLMOs are calculated from subsystems and used for constructing the density matrix (DM) of the entire system, instead of the subsystem DM in the original DC approach. Also, unlike the original DC method, the inverse electronic temperature parameter β is not needed anymore. Furthermore, a new regularized localization approach for NOLMOs has been developed, in which the localization cost function is a sum of the spatial spread function, as in the Boys method, and the kinetic energy, as a regularization measure to limit the oscillation of the NOLMOs. The optimal weight of the kinetic energy can be determined by optimization with analytical gradients. The resulting regularized NOLMOs have enhanced smoothness and better transferability because of reduced kinetic energies. Compared with the original DC, while NOLMO-DC has a similar computational linear scaling cost, the accuracy of NOLMO-DC is better by several orders of magnitude for large conjugated systems and by about 1 order of magnitude for other systems. The NOLMO-DC method is thus a promising development of the DC approach for linear scaling calculations.

Stochastic density functional theory: Real- and energy-space fragmentation for noise reduction

Stochastic density functional theory (sDFT) is becoming a valuable tool for studying ground-state properties of extended materials. The computational complexity of describing the Kohn-Sham orbitals is replaced by introducing a set of random (stochastic) orbitals leading to linear and often sub-linear scaling of certain ground-state observables at the account of introducing a statistical error. Schemes to reduce the noise are essential, for example, for determining the structure using the forces obtained from sDFT. Recently, we have introduced two embedding schemes to mitigate the statistical fluctuations in the electron density and resultant forces on the nuclei. Both techniques were based on fragmenting the system either in real space or slicing the occupied space into energy windows, allowing for a significant reduction in the statistical fluctuations. For chemical accuracy, further reduction of the noise is required, which could be achieved by increasing the number of stochastic orbitals. However, the convergence is relatively slow as the statistical error scales as 1/N_χ according to the central limit theorem, where N_χ is the number of random orbitals. In this paper, we combined the embedding schemes mentioned above and introduced a new approach that builds on overlapped fragments and energy windows. The new approach significantly lowers the noise for ground-state properties, such as the electron density, total energy, and forces on the nuclei, as demonstrated for a G-center in bulk silicon.

CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations

CP2K is an open source electronic structure and molecular dynamics software package to perform atomistic simulations of solid-state, liquid, molecular, and biological systems. It is especially aimed at massively parallel and linear-scaling electronic structure methods and state-of-the-art ab initio molecular dynamics simulations. Excellent performance for electronic structure calculations is achieved using novel algorithms implemented for modern high-performance computing systems. This review revisits the main capabilities of CP2K to perform efficient and accurate electronic structure simulations. The emphasis is put on density functional theory and multiple post-Hartree-Fock methods using the Gaussian and plane wave approach and its augmented all-electron extension.

Integrating Machine Learning with the Multilayer Energy-Based Fragment Method for Excited States of Large Systems

In this work we have combined machine learning techniques with our recently developed multilayer energy-based fragment method for studying excited states of large systems. The photochemically active and inert regions are separately treated with the complete active space self-consistent field method and the trained models. This method is demonstrated to provide accurate energies and gradients leading to essentially the same excited-state potential energy surfaces and nonadiabatic dynamics compared with full ab initio results. Furthermore, in conjunction with the use of machine learning models, this method is highly parallel and exhibits low-scaling computational cost. Finally, the present work could encourage the marriage of machine learning with fragment-based electronic structure methods to explore photochemistry of large systems.

A multi-layer energy-based fragment method for excited states and nonadiabatic dynamics

We developed a multi-layer energy-based fragment (MLEBF) method within the many-body energy expansion framework. It supplies accurate energies and gradients, and accurately reproduces excited-state topological structures. Moreover, MLEBF-based nonadiabatic dynamics simulations give nearly the same results compared with full ab initio ones. The present work could stimulate developing energy-based fragment methods for photochemistry of large systems.

Energy window stochastic density functional theory

Linear scaling density functional theory is important for understanding electronic structure properties of nanometer scale systems. Recently developed stochastic density functional theory can achieve linear or even sublinear scaling for various electronic properties without relying on the sparsity of the density matrix. The basic idea relies on projecting stochastic orbitals onto the occupied space by expanding the Fermi-Dirac operator and repeating this for N_χ stochastic orbitals. Often, a large number of stochastic orbitals are required to reduce the statistical fluctuations (which scale as N_χ ^-1/2) below a tolerable threshold. In this work, we introduce a new stochastic density functional theory that can efficiently reduce the statistical fluctuations for certain observable and can also be integrated with an embedded fragmentation scheme. The approach is based on dividing the occupied space into energy windows and projecting the stochastic orbitals with a single expansion onto all windows simultaneously. This allows for a significant reduction of the noise as illustrated for bulk silicon with a large supercell. We also provide theoretical analysis to rationalize why the noise can be reduced only for a certain class of ground state properties, such as the forces and electron density.

Libraries of Extremely Localized Molecular Orbitals. 3. Construction and Preliminary Assessment of the New Databanks

The fast and reliable determination of wave functions and electron densities of macromolecules has been one of the goals of theoretical chemistry for a long time, and in this context, several linear scaling techniques have been successfully devised over the years. Different approaches have been adopted to tackle this problem, and one of them exploits the fact that, according to the traditional chemical perception, molecules can be seen as constituted of recurring units (e.g., functional groups) with well-defined chemical features. This has led to the development of methods in which the global wave functions or electron densities of macromolecules are obtained by simply transferring density matrices or fuzzy electron densities associated with molecular fragments. In this context, we propose an alternative strategy that aims at quickly reconstructing wave functions and electron densities of proteins through the transfer of extremely localized molecular orbitals (ELMOs), which are orbitals strictly localized on small molecular units and, for this reason, easily transferable from molecule to molecule. To accomplish this task we have constructed original libraries of ELMOs that cover all the possible elementary fragments of the 20 natural amino acids in all their possible protonation states and forms. Our preliminary test calculations have shown that, compared to more traditional methods of quantum chemistry, the transfers from the novel ELMO databanks allow to obtain wave function and electron densities of large polypeptides and proteins at a significantly reduced computational cost. Furthermore, notwithstanding expected discrepancies, the obtained electron distributions and electrostatic potentials are in very good agreement with those obtained at Hartree-Fock and density functional theory (DFT) levels. Therefore, the results encourage to use the new libraries as alternatives to the popular pseudoatom-databases of crystallography in the refinement of crystallographic structures of macromolecules. In particular, in this context, we have already envisaged the coupling of the ELMO databanks with the promising Hirshfeld atom refinement technique to extend the applicability of the latter to very large systems.

Spatial and Spin Symmetry Breaking in Semidefinite-Programming-Based Hartree-Fock Theory

The Hartree-Fock problem was recently recast as a semidefinite optimization over the space of rank-constrained two-body reduced-density matrices (RDMs) [ Phys. Rev. A 2014 , 89 , 010502(R) ]. This formulation of the problem transfers the nonconvexity of the Hartree-Fock energy functional to the rank constraint on the two-body RDM. We consider an equivalent optimization over the space of positive semidefinite one-electron RDMs (1-RDMs) that retains the nonconvexity of the Hartree-Fock energy expression. The optimized 1-RDM satisfies ensemble N-representability conditions, and ensemble spin-state conditions may be imposed as well. The spin-state conditions place additional linear and nonlinear constraints on the 1-RDM. We apply this RDM-based approach to several molecular systems and explore its spatial (point group) and spin ( Ŝ² and Ŝ₃) symmetry breaking properties. When imposing Ŝ² and Ŝ₃ symmetry but relaxing point group symmetry, the procedure often locates spatial-symmetry-broken solutions that are difficult to identify using standard algorithms. For example, the RDM-based approach yields a smooth, spatial-symmetry-broken potential energy curve for the well-known Be-H₂ insertion pathway. We also demonstrate numerically that, upon relaxation of Ŝ² and Ŝ₃ symmetry constraints, the RDM-based approach is equivalent to real-valued generalized Hartree-Fock theory.

Classical and quantum many-body effects on the critical properties and thermodynamic regularities of silicon

Using molecular simulation, we determine the critical properties of Si as well as the loci for several remarkable thermodynamic contours spanning the supercritical region of the phase diagram. We consider a classical three-body potential as well as a quantum (tight-binding) many-body model, and determine the loci for the ideality contours, including the Zeno line and the H line of ideal enthalpy. The two strategies (classical or quantum) lead to strongly asymmetric binodals and to critical properties in good agreement with each other. The Zeno and H lines are found to remain linear over a wide temperature interval, despite the changes in electronic structure undergone by the fluid along these contours. We also show that the classical and quantum model yield markedly different results for the parameters defining the H line, the exponents for the power-laws underlying the line of minima for the isothermal enthalpy and for the density required to achieve ideal behavior, most notably for the enthalpy.

Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. IV. Performance of many-body force fields and tight-binding schemes for the fluid phases of silicon

We extend Expanded Wang-Landau (EWL) simulations beyond classical systems and develop the EWL method for systems modeled with a tight-binding Hamiltonian. We then apply the method to determine the partition function and thus all thermodynamic properties, including the Gibbs free energy and entropy, of the fluid phases of Si. We compare the results from quantum many-body (QMB) tight binding models, which explicitly calculate the overlap between the atomic orbitals of neighboring atoms, to those obtained with classical many-body (CMB) force fields, which allow to recover the tetrahedral organization in condensed phases of Si through, e.g., a repulsive 3-body term that favors the ideal tetrahedral angle. Along the vapor-liquid coexistence, between 3000 K and 6000 K, the densities for the two coexisting phases are found to vary significantly (by 5 orders of magnitude for the vapor and by up to 25% for the liquid) and to provide a stringent test of the models. Transitions from vapor to liquid are predicted to occur for chemical potentials that are 10%-15% higher for CMB models than for QMB models, and a ranking of the force fields is provided by comparing the predictions for the vapor pressure to the experimental data. QMB models also reveal the formation of a gap in the electronic density of states of the coexisting liquid at high temperatures. Subjecting Si to a nanoscopic confinement has a dramatic effect on the phase diagram with, e.g. at 6000 K, a decrease in liquid densities by about 50% for both CMB and QMB models and an increase in vapor densities between 90% (CMB) and 170% (QMB). The results presented here provide a full picture of the impact of the strategy (CMB or QMB) chosen to model many-body effects on the thermodynamic properties of the fluid phases of Si.

Effective preconditioning for ab initio ground state energy minimization with non-orthogonal localized molecular orbitals

The non-orthogonal localized molecular orbital (NOLMO) is the most localized representation of electronic degrees of freedom. As such, NOLMOs are thus potentially the most efficient for linear-scaling calculations of electronic structures for large systems. However, direct ab initio calculations with NOLMO have not been fully implemented and widely used, partly because of the slow convergence issue in the optimization of NOLMO. Towards realizing the potential of NOLMO for large systems, we applied an energy minimum variational principle for carrying out ab initio self-consistent-field (SCF) calculations with NOLMOs. We developed an effective preconditioning approach using the diagonal part of the second order derivatives and show that the convergence of the energy optimization is significantly improved. The speed of convergence of the energy and density are comparable with that of the conventional SCF approach, thus paving the way for the optimization of NOLMO in linear scaling calculations for large systems.

Calculation of space localized properties in correlated quantum Monte Carlo methods with reweighting: the nonlocality of statistical uncertainties

We study the efficiency of quantum Monte Carlo (QMC) methods in computing space localized ground state properties (properties which do not depend on distant degrees of freedom) as a function of the system size N. We prove that for the commonly used correlated sampling with reweighting method, the statistical fluctuations σ2(N) do not obey the locality property. σ2(N) grow at least linearly with N and with a slope that is related to the fluctuations of the reweighting factors. We provide numerical illustrations of these tendencies in the form of QMC calculations on linear chains of hydrogen atoms.

O(N) methods in electronic structure calculations

Linear-scaling methods, or O(N) methods, have computational and memory requirements which scale linearly with the number of atoms in the system, N, in contrast to standard approaches which scale with the cube of the number of atoms. These methods, which rely on the short-ranged nature of electronic structure, will allow accurate, ab initio simulations of systems of unprecedented size. The theory behind the locality of electronic structure is described and related to physical properties of systems to be modelled, along with a survey of recent developments in real-space methods which are important for efficient use of high-performance computers. The linear-scaling methods proposed to date can be divided into seven different areas, and the applicability, efficiency and advantages of the methods proposed in these areas are then discussed. The applications of linear-scaling methods, as well as the implementations available as computer programs, are considered. Finally, the prospects for and the challenges facing linear-scaling methods are discussed.

Real-space finite-difference approach for multi-body systems: path-integral renormalization group method and direct energy minimization method

The path-integral renormalization group and direct energy minimization method of practical first-principles electronic structure calculations for multi-body systems within the framework of the real-space finite-difference scheme are introduced. These two methods can handle higher dimensional systems with consideration of the correlation effect. Furthermore, they can be easily extended to the multicomponent quantum systems which contain more than two kinds of quantum particles. The key to the present methods is employing linear combinations of nonorthogonal Slater determinants (SDs) as multi-body wavefunctions. As one of the noticeable results, the same accuracy as the variational Monte Carlo method is achieved with a few SDs. This enables us to study the entire ground state consisting of electrons and nuclei without the need to use the Born-Oppenheimer approximation. Recent activities on methodological developments aiming towards practical calculations such as the implementation of auxiliary field for Coulombic interaction, the treatment of the kinetic operator in imaginary-time evolutions, the time-saving double-grid technique for bare-Coulomb atomic potentials and the optimization scheme for minimizing the total-energy functional are also introduced. As test examples, the total energy of the hydrogen molecule, the atomic configuration of the methylene and the electronic structures of two-dimensional quantum dots are calculated, and the accuracy, availability and possibility of the present methods are demonstrated.

Assessment of density matrix methods for linear scaling electronic structure calculations

Purification and minimization methods for linear scaling computation of the one-particle density matrix for a fixed Hamiltonian matrix are compared. This is done by considering the work needed by each method to achieve a given accuracy in terms of the difference from the exact solution. Numerical tests employing orthogonal as well as non-orthogonal versions of the methods are performed using both element magnitude and cutoff radius based truncation approaches. It is investigated how the convergence speed for the different methods depends on the eigenvalue distribution in the Hamiltonian matrix. An expression for the number of iterations required for the minimization methods studied is derived, taking into account the dependence on both the band gap and the chemical potential. This expression is confirmed by numerical tests. The minimization methods are found to perform at their best when the chemical potential is located near the center of the eigenspectrum. The results indicate that purification is considerably more efficient than the minimization methods studied even when a good starting guess for the minimization is available. In test calculations without a starting guess, purification is more than an order of magnitude more efficient than minimization.

Equivalent potential of water for the electronic structure of glycine

First-principles, all-electron, ab initio calculations have been performed to construct an equivalent potential of water for the electronic structure of glycine (Gly) in solution. The calculation involved three steps. The first step was to search for the minimum-energy geometric structure of the Gly + nH₂O system. The second step was to calculate the electronic structure of Gly with the potential of water molecules via the self-consistent cluster-embedding method (SCCE), based on the result obtained in the first step. The last step was to calculate the electronic structure of Gly with the potential of dipoles after replacing the water molecules with dipoles. The results show that the occupied molecular orbitals of Gly are raised by about 0.0524 Ry on average due to the effect of water. The effect of water can be simulated well using the dipole potential. The equivalent potential obtained can be applied directly to electronic structure calculations of proteins in solution using the SCCE method.

Novel Computational Methods for Nanostructure Electronic Structure Calculations

Experimentally relevant nanocrystals often contain a few thousands to hundreds of thousands of atoms. Yet, to understand their electronic structures, surface and impurity effects, atomic relaxations, interior electric fields, carrier dynamics, and transports, it is often necessary to carry out atomistic simulations. Owing to the advance of recent algorithm developments and improved supercomputer powers, it is now possible to calculate such nanocrystals based on ab initio methods. In this review, we discuss the numerical algorithms (the plane-wave pseudopotential method and the real-space finite-difference method) used in conventional density-functional-theory calculations, which enable the simulations of systems up to one or two thousand atoms. We also introduce methods designed specifically for nanostructure calculations. These methods [the charge-patching method (CPM) and the linear scaling three-dimensional fragment method (LS3DF)] can be used to calculate systems with hundreds of thousands of atoms. Whereas CPM is an approximation with ab initio quality, the LS3DF method is an O(N) method with essentially the same results as the direct methods. The computational aspects of the algorithms, especially for their parallelization scalability, are also emphasized in the review.

Calculations for millions of atoms with density functional theory: linear scaling shows its potential

An overview of the CONQUEST linear scaling density functional theory (DFT) code is given, focusing particularly on the scaling behaviour on modern high-performance computing platforms. We demonstrate that essentially perfect linear scaling and weak parallel scaling (with fixed number of atoms per processor core) can be achieved, and that DFT calculations on millions of atoms are now possible.

A first-principles divide-and-conquer approach for electronic structure of large systems and its application to graphene nanoribbons

We demonstrate an efficient and accurate first-principles method to calculate the electronic structure of a large system using a divide-and-conquer strategy based on localized quasi-atomic minimal basis set orbitals recently developed. Tight-binding Hamiltonian and overlap matrices of a large system can be constructed by extracting the matrix elements for a given pair of atoms from first-principles calculations of smaller systems that represent the local bonding environment of the particular atom pair. The approach is successfully applied to the studies of electronic structure in graphene nanoribbons. This provides a promising way to do the electronic simulation for large systems directly from first principles.

Total-energy minimization of few-body electron systems in the real-space finite-difference scheme

A practical and high-accuracy computation method to search for ground states of few-electron systems is presented on the basis of the real-space finite-difference scheme. A linear combination of Slater determinants is employed as a many-electron wavefunction, and the total-energy functional is described in terms of overlap integrals of one-electron orbitals without the constraints of orthogonality and normalization. In order to execute a direct energy minimization process of the energy functional, the steepest-descent method is used. For accurate descriptions of integrals which include bare-Coulomb potentials of ions, the time-saving double-grid technique is introduced. As an example of the present method, calculations for the ground state of the hydrogen molecule are demonstrated. An adiabatic potential curve is illustrated, and the accessibility and accuracy of the present method are discussed.

Theory of the low frequency mechanical modes and Raman spectra of the M13 bacteriophage capsid with atomic detail

We present a theoretical study of the low frequency vibrational modes of the M13 bacteriophage using a fully atomistic model. Using ideas from electronic structure theory, the few lowest vibrational modes of the M13 bacteriophage are determined using classical harmonic analysis. The relative Raman intensity is estimated for each of the mechanical modes using a bond polarizability model. Comparison of the atomic mechanical modes calculated here with modes derived from elastic continuum theory shows that a much richer spectrum emerges from an atomistic picture.