Performance of the density matrix functional theory in the quantum theory of atoms in molecules

The generalization to arbitrary molecular geometries of the energetic partitioning provided by the atomic virial theorem of the quantum theory of atoms in molecules (QTAIM) leads to an exact and chemically intuitive energy partitioning scheme, the interacting quantum atoms (IQA) approach, that depends on the availability of second-order reduced density matrices (2-RDMs). This work explores the performance of this approach in particular and of the QTAIM in general with approximate 2-RDMs obtained from the density matrix functional theory (DMFT), which rests on the natural expansion (natural orbitals and their corresponding occupation numbers) of the first-order reduced density matrix (1-RDM). A number of these functionals have been implemented in the promolden code and used to perform QTAIM and IQA analyses on several representative molecules and model chemical reactions. Total energies, covalent intra- and interbasin exchange-correlation interactions, as well as localization and delocalization indices have been determined with these functionals from 1-RDMs obtained at different levels of theory. Results are compared to the values computed from the exact 2-RDMs, whenever possible.     

Phase dilemma in density matrix functional theory

For closed-shell systems, a particular parametrization of coefficients in a configuration interaction (CI) expansion provides a convenient formulation for the search over electronic wave functions constrained by a set of natural orbitals (NOs) and the corresponding occupation numbers that are invoked in every variational construction of the density matrix functional (DMF) V(ee)(Gamma) for the electron--electron repulsion energy. It produces an explicit expression for V(ee) in terms of the Coulomb and exchange integrals over NOs, and an idempotent matrix omega, diagonal elements of which equal the occupation numbers. At the same time, it reveals a very serious bottleneck affecting any rigorous approach to the DMF theory, namely the phase dilemma that stems from the necessity to carry out minimization over a large number of possible combinations of CI coefficient signs. While underscoring its lack of variational nature, a simple approximation for the phase factor products provides a strict derivation for the recently proposed Kollmar-Hess functional.     

The adiabatic approximation in time-dependent density matrix functional theory: response properties from dynamics of phase-including natural orbitals

The adiabatic approximation is problematic in time-dependent density matrix functional theory. With pure density matrix functionals (invariant under phase change of the natural orbitals) it leads to lack of response in the occupation numbers, hence wrong frequency dependent responses, in particular α(ω→0)≠α(0) (the static polarizability). We propose to relinquish the requirement that the functional must be a pure one-body reduced density matrix (1RDM) functional, and to introduce additional variables which can be interpreted as phases of the one-particle states of the independent particle reference system formed with the natural orbitals, thus obtaining so-called phase-including natural orbital (PINO) functionals. We also stress the importance of the correct choice of the complex conjugation in the two-electron integrals in the commonly used functionals (they should not be of exchange type). We demonstrate with the Lo?wdin-Shull energy expression for two-electron systems, which is an example of a PINO functional, that for two-electron systems exact responses (polarizabilities, excitation energies) are obtained, while writing this energy expression in the usual way as a 1RDM functional yields erroneous responses.     

Analytic second-order energy derivatives in natural orbital functional theory

The analytic energy gradients in the atomic orbital representation have recently been published (Mitxelena and Piris in J Chem Phys 146:014102, 2017) within the framework of the natural orbital functional theory (NOFT). We provide here an alternative expression for them in terms of natural orbitals, and use it to derive the analytic second-order energy derivatives with respect to nuclear displacements in the NOFT. The computational burden is shifted to the calculation of perturbed natural orbitals and occupancies, since a set of linear coupled-perturbed equations obtained from the variational Euler equations must be solved to attain the analytic Hessian at the perturbed geometry. The linear response of both natural orbitals and occupation numbers to nuclear geometry displacements need only specify the reconstruction of the second-order reduced density matrix in terms of occupation numbers.     

Comparison of two genres for linear scaling in density functional theory: purification and density matrix minimization methods

Two classes of linear-scaling methods to replace diagonalization of the one-particle Hamiltonian matrix in density functional theory are compared to each other. Purification takes a density matrix with the correct eigenfunctions and corrects the occupation numbers; density matrix minimization takes a density matrix with correct occupation numbers and corrects the eigenfunctions by rotating the orbitals. Computational comparisons are performed through modification of the MondoSCF program on water clusters and the protein endothelin. A purification scheme and a density matrix minimization scheme, based on the 1,2-contracted Schrodinger equation [D. A. Mazziotti, J. Chem. Phys. 115, 8305 (2001)] are implemented in large systems.     

A novel interpretation of reduced density matrix and cumulant for electronic structure theories

We propose a novel interpretation of the reduced density matrix (RDM) and its cumulant that combines linear and exponential parametrizations of the wavefunction. Any n-particle RDM can be written as a weighted average of "configuration interaction" amplitudes. The corresponding n-particle cumulant is represented in terms of two types of contributions: "connected" (statistical averages of substitution amplitudes) and "disconnected" (cross-correlations of substitution amplitudes). A diagonal element of n-RDM represents the average occupation number of the orbital n-tuple. The diagonal elements of 2- and 3-cumulants take particularly elegant forms in the natural spin-orbital basis: they represent the covariances (correlated fluctuations) of the occupation numbers of the orbital pair and triples, respectively. Thus, the diagonal elements of the cumulants quantify the correlation between the orbital occupation numbers. Our interpretation is used to examine the weak to strong correlation transition in the "two electrons in two orbitals" problem.     

Response calculations based on an independent particle system with the exact one-particle density matrix: excitation energies

Adiabatic response time-dependent density functional theory (TDDFT) suffers from the restriction to basically an occupied → virtual single excitation formulation. Adiabatic time-dependent density matrix functional theory allows to break away from this restriction. Problematic excitations for TDDFT, viz. bonding-antibonding, double, charge transfer, and higher excitations, are calculated along the bond-dissociation coordinate of the prototype molecules H(2) and HeH(+) using the recently developed adiabatic linear response phase-including (PI) natural orbital theory (PINO). The possibility to systematically increase the scope of the calculation from excitations out of (strongly) occupied into weakly occupied ("virtual") natural orbitals to larger ranges of excitations is explored. The quality of the PINO response calculations is already much improved over TDDFT even when the severest restriction is made, to virtually the size of the TDDFT diagonalization problem (only single excitation out of occupied orbitals plus all diagonal doubles). Further marked improvement is obtained with moderate extension to allow for excitation out of the lumo and lumo+1, which become fractionally occupied in particular at longer distances due to left-right correlation effects. In the second place the interpretation of density matrix response calculations is elucidated. The one-particle reduced density matrix response for an excitation is related to the transition density matrix to the corresponding excited state. The interpretation of the transition density matrix in terms of the familiar excitation character (single excitations, double excitations of various types, etc.) is detailed. The adiabatic PINO theory is shown to successfully resolve the problematic cases of adiabatic TDDFT when it uses a proper PI orbital functional such as the PILS functional.     

Iterative diagonalization for orbital optimization in natural orbital functional theory

A challenging task in natural orbital functional theory is to find an efficient procedure for doing orbital optimization. Procedures based on diagonalization techniques have confirmed its practical value since the resulting orbitals are automatically orthogonal. In this work, a new procedure is introduced, which yields the natural orbitals by iterative diagonalization of a Hermitian matrix F . The off‐diagonal elements of the latter are determined explicitly from the hermiticity of the matrix of the Lagrange multipliers. An expression for diagonal elements is absent so a generalized Fockian is undefined in the conventional sense, nevertheless, they may be determined from an aufbau principle. Thus, the diagonal elements are obtained iteratively considering as starting values those coming from a single diagonalization of the matrix of the Lagrange multipliers calculated with the Hartree‐Fock orbitals after the occupation numbers have been optimized. The method has been tested on the G2/97 set of molecules for the Piris natural orbital functional. To help the convergence, we have implemented a variable scaling factor which avoids large values of the off‐diagonal elements of F . The elapsed times of the computations required by the proposed procedure are compared with a full sequential quadratic programming optimization, so that the efficiency of the method presented here is demonstrated. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

The exchange-correlation potential in ab initio density functional theory

From coupled-cluster theory and many-body perturbation theory we derive the local exchange-correlation potential of density functional theory in an orbital dependent form. We show the relationship between the coupled-cluster approach and density functional theory, and connections and comparisons with our previous second-order correlation potential [OEP-MBPT(2) (OEP-optimized effective potential)] [I. Grabowski, S. Hirata, S. Ivanov, and R. J. Bartlett, J. Chem. Phys. 116, 4415 (2002)]. Starting from a general theoretical framework based on the density condition in Kohn-Sham theory, we define a rigorous exchange-correlation functional, potential and orbitals. Specifying initially to second-order terms, we show that our ab initio correlation potential provides the correct shape compared to those from reference quantum Monte Carlo calculations, and we demonstrate the superiority of using Fock matrix elements or more general infinite-order semicanonical transformations. This enables us to introduce a method that is guaranteed to converge to the right answer in the correlation and basis set limit, just as does ab initio wave function theory. We also demonstrate that the energies obtained from this generalized second-order method [OEP-MBPT2-f] and [OEP-MBPT2-sc] are often of coupled-cluster accuracy and substantially better than ordinary Hartree-Fock based second-order MBPT=MP2.     

Adiabatic approximation of time-dependent density matrix functional response theory

Time-dependent density matrix functional theory can be formulated in terms of coupled-perturbed response equations, in which a coupling matrix K(omega) features, analogous to the well-known time-dependent density functional theory (TDDFT) case. An adiabatic approximation is needed to solve these equations, but the adiabatic approximation is much more critical since there is not a good "zero order" as in TDDFT, in which the virtual-occupied Kohn-Sham orbital energy differences serve this purpose. We discuss a simple approximation proposed earlier which uses only results from static calculations, called the static approximation (SA), and show that it is deficient, since it leads to zero response of the natural orbital occupation numbers. This leads to wrong behavior in the omega-->0 limit. An improved adiabatic approximation (AA) is formulated. The two-electron system affords a derivation of exact coupled-perturbed equations for the density matrix response, permitting analytical comparison of the adiabatic approximation with the exact equations. For the two-electron system also, the exact density matrix functional (2-matrix in terms of 1-matrix) is known, enabling testing of the static and adiabatic approximations unobscured by approximations in the functional. The two-electron HeH(+) molecule shows that at the equilibrium distance, SA consistently underestimates the frequency-dependent polarizability alpha(omega), the adiabatic TDDFT overestimates alpha(omega), while AA improves upon SA and, indeed, AA produces the correct alpha(0). For stretched HeH(+), adiabatic density matrix functional theory corrects the too low first excitation energy and overpolarization of adiabatic TDDFT methods and exhibits excellent agreement with high-quality CCSD ("exact") results over a large omega range.     

Direct minimization of the energy in density functional theory

The minimization of the energy functional of the first-order density matrix γ( r , r ') is achieved using unitary transformations applied to γ. Equivalently, such transformations can be carried out also on one-electron orbitals (natural orbitals) and their occupation (integer or non-integer) numbers. The conventional local density approximation based on the electron density p( r ) is then considered as a special case. The direct minimization of the energy functional of p with respect to the parameters of the unitary transformation leads to stationary conditions that are all equivalent to the Kohn–Sham equations. Preliminary numerical tests show that the proposed algorithms for the direct minimization of the energy work in a satisfactory manner. © John Wiley & Sons, Inc.     

Evaluation of the second-order reduced density matrix for correlated electronic wave functions

Formulas are presented for the evaluation of the second-order reduced density matrix for correlated wave functions, including wave functions constructed from pairwise nonorthogonal orbitals. Numerical results are provided for the correlation holes and conditional nuclear spin densities for the water molecule.     

Lagrange function method for energy optimization directly in the space of natural orbitals

We propose a Lagrange method to obtain the electronic energy directly in the space of natural orbitals. The Lagrange function contains constraints to force the off‐diagonal elements of the one‐particle density matrix to zero, so the molecular orbitals converge to the natural orbitals simultaneously while the energy minimizes. The recently discovered “generalized Pauli conditions” for the occupation numbers are invoked to generate an initial approximation. As a demonstration of this approach, we study the system of three electrons in six orbitals which has become a paradigm for investigating multi‐particle entanglement.     

Excitation energies from range-separated time-dependent density and density matrix functional theory

Time-dependent density functional theory (TD-DFT) in the adiabatic formulation exhibits known failures when applied to predicting excitation energies. One of them is the lack of the doubly excited configurations. On the other hand, the time-dependent theory based on a one-electron reduced density matrix functional (time-dependent density matrix functional theory, TD-DMFT) has proven accurate in determining single and double excitations of H(2) molecule if the exact functional is employed in the adiabatic approximation. We propose a new approach for computing excited state energies that relies on functionals of electron density and one-electron reduced density matrix, where the latter is applied in the long-range region of electron-electron interactions. A similar approach has been recently successfully employed in predicting ground state potential energy curves of diatomic molecules even in the dissociation limit, where static correlation effects are dominating. In the paper, a time-dependent functional theory based on the range-separation of electronic interaction operator is rigorously formulated. To turn the approach into a practical scheme the adiabatic approximation is proposed for the short- and long-range components of the coupling matrix present in the linear response equations. In the end, the problem of finding excitation energies is turned into an eigenproblem for a symmetric matrix. Assignment of obtained excitations is discussed and it is shown how to identify double excitations from the analysis of approximate transition density matrix elements. The proposed method used with the short-range local density approximation (srLDA) and the long-range Buijse-Baerends density matrix functional (lrBB) is applied to H(2) molecule (at equilibrium geometry and in the dissociation limit) and to Be atom. The method accounts for double excitations in the investigated systems but, unfortunately, the accuracy of some of them is poor. The quality of the other excitations is in general much better than that offered by TD-DFT-LDA or TD-DMFT-BB approximations if the range-separation parameter is properly chosen. The latter remains an open problem.     

Inequalities for fermion density matrices

A partial trace over the occupation numbers of all but k states in the density matrix of an ensemble with an arbitrary number of single-particle states is defined as the (reduced) k-state density matrix. This matrix is used to obtain a complete, practical solution to the problem of determining the representability of the diagonal elements of the one- and two-particle (reduced) density matrices. This solution is expressed as a series of linear inequalities involving the density-matrix elements; the inequalities are identical with those derived previously by Davidson and McCrae by a different method. In addition, our method is used to obtain nonlinear, matrix inequalities on the off-diagonal elements of the density matrices.     

Density matrix in the open shell theory

All the second-order density matrix spin components for the spin-projected Slater determinant are presented as expansions in direct products of powers of unprojected spin- and residual electron density matrices. Spin components of the second-order transition density matrices between spin-projected Slater determinants built of orthogonal orbitals are also obtained.