Linear-scaling implementation of molecular response theory in self-consistent field electronic-structure theory

A linear-scaling implementation of Hartree-Fock and Kohn-Sham self-consistent field theories for the calculation of frequency-dependent molecular response properties and excitation energies is presented, based on a nonredundant exponential parametrization of the one-electron density matrix in the atomic-orbital basis, avoiding the use of canonical orbitals. The response equations are solved iteratively, by an atomic-orbital subspace method equivalent to that of molecular-orbital theory. Important features of the subspace method are the use of paired trial vectors (to preserve the algebraic structure of the response equations), a nondiagonal preconditioner (for rapid convergence), and the generation of good initial guesses (for robust solution). As a result, the performance of the iterative method is the same as in canonical molecular-orbital theory, with five to ten iterations needed for convergence. As in traditional direct Hartree-Fock and Kohn-Sham theories, the calculations are dominated by the construction of the effective Fock/Kohn-Sham matrix, once in each iteration. Linear complexity is achieved by using sparse-matrix algebra, as illustrated in calculations of excitation energies and frequency-dependent polarizabilities of polyalanine peptides containing up to 1400 atoms.     

Higher order alchemical derivatives from coupled perturbed self-consistent field theory

We present an analytical approach to treat higher order derivatives of Hartree-Fock (HF) and Kohn-Sham (KS) density functional theory energy in the Born-Oppenheimer approximation with respect to the nuclear charge distribution (so-called alchemical derivatives). Modified coupled perturbed self-consistent field theory is used to calculate molecular systems response to the applied perturbation. Working equations for the second and the third derivatives of HF/KS energy are derived. Similarly, analytical forms of the first and second derivatives of orbital energies are reported. The second derivative of Kohn-Sham energy and up to the third derivative of Hartree-Fock energy with respect to the nuclear charge distribution were calculated. Some issues of practical calculations, in particular the dependence of the basis set and Becke weighting functions on the perturbation, are considered. For selected series of isoelectronic molecules values of available alchemical derivatives were computed and Taylor series expansion was used to predict energies of the "surrounding" molecules. Predicted values of energies are in unexpectedly good agreement with the ones computed using HF/KS methods. Presented method allows one to predict orbital energies with the error less than 1% or even smaller for valence orbitals.     

Communication: Monte Carlo calculation of the exchange energy

In recent generalized Kohn-Sham (GKS) schemes for density functional theory (DFT) Hartree-Fock type exchange is important. In plane waves and grid approaches the high cost of exchange energy calculations makes these GKS considerably more expensive than Kohn-Sham DFT calculations. We develop a stochastic approach for speeding up the calculation of exchange for large systems. We show that stochastic error per particle does not grow and can even decrease with system size (at a given number of iterations). We discuss several alternative approaches and explain how these ideas can be included in the GKS framework.     

A hierarchic sparse matrix data structure for large-scale Hartree-Fock/Kohn-Sham calculations

A hierarchic sparse matrix data structure for Hartree-Fock/Kohn-Sham calculations is presented. The data structure makes the implementation of matrix manipulations needed for large systems faster, easier, and more maintainable without loss of performance. Algorithms for symmetric matrix square and inverse Cholesky decomposition within the hierarchic framework are also described. The presented data structure is general; in addition to its use in Hartree-Fock/Kohn-Sham calculations, it may also be used in other research areas where matrices with similar properties are encountered. The applicability of the data structure to ab initio calculations is shown with help of benchmarks on water droplets and graphene nanoribbons.     

Electric dipole polarizabilities and C6 dipole-dipole dispersion coefficients for sodium clusters and C60

The frequency-dependent polarizabilities of closed-shell sodium clusters containing up to 20 atoms have been calculated using the linear complex polarization propagator approach in conjunction with Hartree-Fock and Kohn-Sham density functional theories. In combination with polarizabilities for C(60) from a previous work [J. Chem. Phys. 123, 124312 (2005)], the C(6) dipole-dipole dispersion coefficients for the metal-cluster-to-cluster and cluster-to-buckminster-fullerene interactions are obtained via the Casimir-Polder relation [Phys. Rev. 73, 360 (1948)]. The B3PW91 results for the polarizability of the sodium dimer and tetramer are benchmarked against coupled cluster calculations. The error bars of the reported theoretical results for the C(6) coefficients are estimated to be 5%, and the results are well within the error bars of the experiment.     

Intracule densities in the strong-interaction limit of density functional theory

The correlation energy in density functional theory can be expressed exactly in terms of the change in the probability of finding two electrons at a given distance r(12) (intracule density) when the electron-electron interaction is multiplied by a real parameter lambda varying between 0 (Kohn-Sham system) and 1 (physical system). In this process, usually called adiabatic connection, the one-electron density is (ideally) kept fixed by a suitable local one-body potential. While an accurate intracule density of the physical system can only be obtained from expensive wavefunction-based calculations, being able to construct good models starting from Kohn-Sham ingredients would highly improve the accuracy of density functional calculations. To this purpose, we investigate the intracule density in the lambda --> infinity limit of the adiabatic connection. This strong-interaction limit of density functional theory turns out to be, like the opposite non-interacting Kohn-Sham limit, mathematically simple and can be entirely constructed from the knowledge of the one-electron density. We develop here the theoretical framework and, using accurate correlated one-electron densities, we calculate the intracule densities in the strong interaction limit for few atoms. Comparison of our results with the corresponding Kohn-Sham and physical quantities provides useful hints for building approximate intracule densities along the adiabatic connection of density functional theory.     

Relativistic calculation of indirect NMR spin-spin couplings using the Douglas-Kroll-Hess approximation

We have employed the Douglas-Kroll-Hess approximation to derive the perturbative Hamiltonians involved in the calculation of NMR spin-spin couplings in molecules containing heavy elements. We have applied this two-component quasirelativistic approach using finite perturbation theory in combination with a generalized Kohn-Sham code that includes the spin-orbit interaction self-consistently and works with Hartree-Fock and both pure and hybrid density functionals. We present numerical results for one-bond spin-spin couplings in the series of tetrahydrides CH(4), SiH(4), GeH(4), and SnH(4). Our two-component Hartree-Fock results are in good agreement with four-component Dirac-Hartree-Fock calculations, although a density-functional treatment better reproduces the available experimental data.     

Real-space post-Hartree-Fock correlation models

Real-space models of nondynamical correlation between electrons of opposite spin and of parallel spin in multicenter (molecular) systems are discussed. These models are designed to be partnered with Hartree-Fock or exact Kohn-Sham exchange. Thus the numerous and well-known problems of local density-functional exchange-correlation approximations, especially in stretched odd-electron systems, are circumvented.     

The trust-region self-consistent field method in Kohn-Sham density-functional theory

The trust-region self-consistent field (TRSCF) method is extended to the optimization of the Kohn-Sham energy. In the TRSCF method, both the Roothaan-Hall step and the density-subspace minimization step are replaced by trust-region optimizations of local approximations to the Kohn-Sham energy, leading to a controlled, monotonic convergence towards the optimized energy. Previously the TRSCF method has been developed for optimization of the Hartree-Fock energy, which is a simple quadratic function in the density matrix. However, since the Kohn-Sham energy is a nonquadratic function of the density matrix, the local energy functions must be generalized for use with the Kohn-Sham model. Such a generalization, which contains the Hartree-Fock model as a special case, is presented here. For comparison, a rederivation of the popular direct inversion in the iterative subspace (DIIS) algorithm is performed, demonstrating that the DIIS method may be viewed as a quasi-Newton method, explaining its fast local convergence. In the global region the convergence behavior of DIIS is less predictable. The related energy DIIS technique is also discussed and shown to be inappropriate for the optimization of the Kohn-Sham energy.     

Performance of the Harris functional for extended basis sets at the Hartree-Fock and density functional levels

The Harris functional is a noniterative variational procedure that uses an input charge density to produce an energy that is surprisingly accurate compared to the converged Kohn-Sham self-consistent result. We adapted and generalized this functional for the Hartree-Fock closed- and open-shell cases as well as examined its use for hybrid density functional methods such as B3LYP. Analysis of extended basis set calculations shows that at the B3LYP level an input density formed from a double zeta + polarization orbital basis is accurate enough to reproduce the energy of triple zeta + double polarization + diffuse orbital basis. For large molecules this translates into a computational speed that can be an order of magnitude faster. In the case of Hartree-Fock calculations a "bootstrapping technique" that employs successive applications of the Harris functional can further reduce computational times while retaining sufficient accuracy.     

Quantum chemical study of the electronic structure of NiCH2 + in its ground state and low-lying electronic excited states

The electronic structure of NiCH(2) (+), representative of transition metal carbene ions, is investigated by means of several methods of quantum chemistry. The relative stabilities of the four low-lying doublet electronic states ((2)A(1), (2)A(2), (2)B(1), and (2)B(2)) are determined at the coupled cluster singles and doubles level (CCSD) and triples level [CCSD(T) and CCSDT-3] with both a Hartree-Fock and density functional theory (Kohn-Sham) reference. The equation-of-motion coupled cluster for treatment of excited states in singles and doubles approximation (EOM-CCSD) is used to characterize the transition energies from the (2)A(1) electronic ground state to the low-lying doublet excited states. The (2)A(2) and (2)B(1) states are nearly degenerate, found to be separated by 940 cm(-1) at the EOM-CCSD level, in agreement with the CASSCF energy ordering. The (2)B(2) state is calculated to be higher in energy by more than 1.0 eV. The spin purity of the low-lying doublet and quadruplet states described by CCSD calculations based on the unrestricted open-shell Hartree-Fock reference is discussed.     

A density matrix-based method for the linear-scaling calculation of dynamic second- and third-order properties at the Hartree-Fock and Kohn-Sham density functional theory levels

A density matrix-based time-dependent self-consistent field (D-TDSCF) method for the calculation of dynamic polarizabilities and first hyperpolarizabilities using the Hartree-Fock and Kohn-Sham density functional theory approaches is presented. The D-TDSCF method allows us to reduce the asymptotic scaling behavior of the computational effort from cubic to linear for systems with a nonvanishing band gap. The linear scaling is achieved by combining a density matrix-based reformulation of the TDSCF equations with linear-scaling schemes for the formation of Fock- or Kohn-Sham-type matrices. In our reformulation only potentially linear-scaling matrices enter the formulation and efficient sparse algebra routines can be employed. Furthermore, the corresponding formulas for the first hyperpolarizabilities are given in terms of zeroth- and first-order one-particle reduced density matrices according to Wigner's (2n+1) rule. The scaling behavior of our method is illustrated for first exemplary calculations with systems of up to 1011 atoms and 8899 basis functions.     

Rydberg states with quantum Monte Carlo

Calculations on Rydberg states are performed using quantum Monte Carlo methods. Excitation energies and singlet-triplet splittings are calculated for two model systems, the carbon atom (3P and 1P) and carbon monoxide ((1Sigma and 3Sigma). Kohn-Sham wave functions constructed from open-shell localized Hartree-Fock orbitals are used as trial and guide functions. The fixed-node diffusion quantum Monte Carlo (FN-DMC) method depends strongly on the wave function's nodal hypersurface. Nodal artefacts are investigated for the ground state of the carbon atom. Their effect on the FN-DMC results can be analyzed quantitatively. FN-DMC leads to accurate excitation energies but to less accurate singlet-triplet splittings. Variational Monte Carlo calculations are able to reproduce the experimental results for both the excitation energies and the singlet-triplet splittings.     

One- and two-center physical space partitioning of the energy in the density functional theory

A conceptually new approach is introduced for the decomposition of the molecular energy calculated at the density functional theory level of theory into sum of one- and two-atomic energy components, and is realized in the "fuzzy atoms" framework. (Fuzzy atoms mean that the three-dimensional physical space is divided into atomic regions having no sharp boundaries but exhibiting a continuous transition from one to another.) The new scheme uses the new concept of "bond order density" to calculate the diatomic exchange energy components and gives them unexpectedly close to the values calculated by the exact (Hartree-Fock) exchange for the same Kohn-Sham orbitals.     

Exchange methods in Kohn-Sham theory

Differences between exchange methods in exchange-only Kohn-Sham theory are highlighted by calculations of diatomic molecule total energies, uncoupled isotropic NMR shieldings, and HOMO-LUMO eigenvalue differences. Optimised effective potential (OEP) and Wu-Yang (WY) results are very similar. Localised Hartree-Fock (LHF) and Krieger-Li-Iafrate (KLI) results are close to one another, but are different to OEP and WY. Becke 1988 exchange (B88X) is different again. Shieldings reduce from OEP/WY to LHF/KLI to B88X, which is consistent with an observed reduction in HOMO-LUMO gaps. LHF, KLI, and B88X shieldings and HOMO-LUMO gaps are closer to near-exact, correlated values, than are the OEP values. These variations arise entirely due to differences in the one-electron exchange potentials, which is clearly evident in potential difference plots, relative to OEP, for the N2 molecule. Density difference plots are also presented, which exhibit a spatial correlation with the potential differences. HOMO and LUMO probability density difference plots show a contraction of the LUMO relative to OEP, which is consistent with the NMR and HOMO-LUMO findings. Plots are also presented for near-exact, correlated Kohn-Sham calculations. The features are qualitatively similar to those observed in the LHF, KLI, and B88X plots, highlighting correlated character in these approximate exchange-only calculations.     

Density functional LCAO calculations for solids: Comparison among Hartree–Fock,DFT local density approximation,and DFT generalized gradient approximation structural properties

In this article, the results of a recently implemented DFT a posteriori and Kohn-Sham (KS ) linear combination of atomic orbital computational scheme for solids are presented. The equilibrium lattice parameters, bulk moduli, and lattice energies are calculated for eight crystallized systems. Local density approximation (LDA ) and generalized gradient approximation (GCA ) functionals and potentials are used. The maps of the Hartree-Fock (HF ) and Ks electronic densities and band structures are depicted. The KS results confirm the trend of the a posteriori scheme. Very good agreement between calculated and experimental lattice energies has been found for GGA potentials. © 1995 John Wiley & Sons, Inc.     

Comparison of the Hartree-Fock,Møller-Plesset,and Hartree-Fock-Slater method with respect to electrostatic properties of small molecules

Summary The results of various quantum chemical calculations, the Hartree-Fock (HF) method, the Møller-Plesset perturbation theory (MP2), and the Hartree-Fock-Slater (HFS) method are compared. Atomic charges, dipole moments, topological properties of the electron density distribution and polarizabilities, and first hyperpolarizabilities are calculated. Atomic charges obtained with the HFS method are found to be very close to those calculated with the MP2 method, from which we conclude that the HFS method describes to some extent electron correlation effects. Performing an MP2 calculation after an HF calculation improves the molecular dipole moments considerably, yielding values close to the experimental ones. HFS calculations are computationally less demanding than MP2 and yield comparable results for the electron density distributions, dipole moments and polarizabilities.     

Natural transition orbitals for complex two-component excited state calculations

While the natural transition orbital (NTO) method has allowed electronic excitations from time-dependent Hartree-Fock and density functional theory to be viewed in a traditional orbital picture, the extension to multicomponent molecular orbitals such as those used in relativistic two-component methods or generalized Hartree-Fock (GHF) or generalized Kohn-Sham (GKS) is less straightforward due to mixing of spin-components and the inherent inclusion of spin-flip transitions in time-dependent GHF/GKS. An extension of single-component NTOs to the two-component framework is presented, in addition to a brief discussion of the practical aspects of visualizing two-component complex orbitals. Unlike the single-component analog, the method explicitly describes the spin and frequently obtains solutions with several significant orbital pairs. The method is presented using calculations on a mercury atom and a CrO₂Cl₂ complex.