A reformulation of the coupled perturbed self-consistent field equations entirely within a local atomic orbital density matrix-based scheme

To exploit the exponential decay found in numerical studies for the density matrix and its derivative with respect to nuclear displacements, we reformulate the coupled perturbed self-consistent field (CPSCF) equations and a quadratically convergent SCF (QCSCF) method for Hartree-Fock and density functional theory within a local density matrix-based scheme. Our D-CPSCF (density matrix-based CPSCF) and D-QCSCF schemes open the way for exploiting sparsity and to achieve asymptotically linear scaling of computational complexity with molecular size (M), in case of D-CPSCF for all (M) derivative densities. Furthermore, these methods are even for small molecules strongly competitive to conventional algorithms.     

Linear-scaling method for calculating nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals within Hartree-Fock and density-functional theory

Details of a new density matrix-based formulation for calculating nuclear magnetic resonance chemical shifts at both Hartree-Fock and density functional theory levels are presented. For systems with a nonvanishing highest occupied molecular orbital-lowest unoccupied molecular orbital gap, the method allows us to reduce the asymptotic scaling order of the computational effort from cubic to linear, so that molecular systems with 1000 and more atoms can be tackled with today's computers. The key feature is a reformulation of the coupled-perturbed self-consistent field (CPSCF) theory in terms of the one-particle density matrix (D-CPSCF), which avoids entirely the use of canonical MOs. By means of a direct solution for the required perturbed density matrices and the adaptation of linear-scaling integral contraction schemes, the overall scaling of the computational effort is reduced to linear. A particular focus of our formulation is to ensure numerical stability when sparse-algebra routines are used to obtain an overall linear-scaling behavior.     

The ab initio calculation of molecular electric, magnetic and geometric properties

We give an account of some recent advances in the development of ab initio methods for the calculation of molecular response properties, involving electric, magnetic, and geometric perturbations. Particular attention is given to properties in which the basis functions depend explicitly both on time and on the applied perturbations such as perturbations involving nuclear displacements or external magnetic fields when London atomic orbitals are used. We summarize a general framework based on the quasienergy for the calculation of arbitrary-order molecular properties using the elements of the density matrix in the atomic-orbital basis as the basic variables. We demonstrate that the necessary perturbed density matrices of arbitrary order can be determined from a set of linear equations that have the same formal structure as the set of linear equations encountered when determining the linear response equations (or time-dependent self-consistent-field equations). Additional components needed to calculate properties involving perturbation-dependent basis sets are flexible one- and two-electron integral techniques for geometric or magnetic-field differentiated integrals; in Kohn-Sham density-functional theory (KS-DFT), we also need to calculate derivatives of the exchange-correlation functional. We describe a recent proposal for evaluating these contributions based on automatic differentiation. Within this framework, it is now possible to calculate any molecular property for an arbitrary self-consistent-field reference state, including two- and four-component relativistic self-consistent-field wave functions. Examples of calculations that can be performed with this formulation are presented.     

Self-consistent, constrained linear-combination-of-atomic-potentials approach to quantum mechanics

Variational fitting gives a stationary linear-combination of atomic potentials (LCAP) approximation to the Kohn-Sham (KS) potential, V. That potential is central to density-functional theory because it generates all orbitals, occupied as well as virtual. Perturbation theory links two self-consistent field (SCF) calculations that differ by the perturbation. Using the same variational LCAP methods and basis sets in the two SCF calculations gives precise KS potentials for each order. Variational V perturbation theory, developed herein through second order, gives stationary potentials at each order and stationary even-order perturbed energies that precisely link the two SCF calculations. Iterative methods are unnecessary because the dimension of the matrix that must be inverted is the KS basis size, not the number of occupied times virtual orbitals of coupled-perturbed methods. With variational perturbation theory, the precision of derivatives and the fidelity of the LCAP KS potential are not related. Finite differences of SCF calculations allow the precision of analytic derivatives from double-precision code to be verified to roughly seven significant digits. For a simple functional, the fourth derivatives of the energy and the first and second derivative of the KS potentials with respect to orbital occupation are computed for a standard set of molecules and basis sets, with and without constraints on the fit to the KS potential. There is no significant difference between the constrained and unconstrained calculations.     

Molecular symmetry. IV. The coupled perturbed Hartree–Fock method

Symmetry methods employed in the ab initio polyatomic program HONDO are extended to the coupled perturbed Hartree–Fock (CPHF) formalism, a key step in the analytical computation of energy first derivatives for configuration interaction (CI) wavefunctions, and energy second derivatives for Hartree–Fock (HF) wavefunctions. One possible computational strategy is to construct Fock-like matrices for each nuclear coordinate in which the one- and two-electron integrals of the usual Fock matrix are replaced by the integral first derivatives. “Skeleton” matrices are constructed from the unique blocks of electron-repulsion integral derivatives. The correct matrices are generated by applying a symmetrization operator. The analysis is valid for many wavefunctions, including closed- or open-shell spin-restricted and spin-unrestricted HF wavefunctions. To illustrate the method, we compare the computer time required for setting up the coupled perturbed HF equations for eclipsed ethane using D_3h symmetry point group and various subgroups of D_3h. Computational times are roughly inversely proportional to the order of the point group.     

Auxiliary density perturbation theory

A new approach, named auxiliary density perturbation theory, for the calculation of second energy derivatives is presented. It is based on auxiliary density functional theory in which the Coulomb and exchange-correlation potentials are expressed by auxiliary function densities. Different to conventional coupled perturbed Kohn-Sham equations the perturbed density matrix is obtained noniteratively by solving an inhomogeneous equation system with the dimension of the auxiliary function set used to expand the auxiliary function density. A prototype implementation for the analytic calculation of molecular polarizabilities is presented. It is shown that the polarizabilities obtained with the newly developed auxiliary density perturbation approach match quantitative with the ones from standard density functional theory if augmented auxiliary function sets are used. The computational advantages of auxiliary density perturbation theory are discussed, too.     

An exact reformulation of the diagonalization step in electronic structure calculations as a set of second order nonlinear equations

A new formulation of the diagonalization step in self-consistent-field (SCF) electronic structure calculations is presented. It exactly replaces the diagonalization of the effective Hamiltonian with the solution of a set of second order nonlinear equations. The density matrix and/or the new set of occupied orbitals can be directly obtained from the resulting solution. This formulation may offer interesting possibilities for new approaches to efficient SCF calculations. The working equations can be derived either from energy minimization with respect to a Cayley-type parametrization of a unitary matrix, or from a similarity transformation approach.     

Study of static and dynamic first hyperpolarizabilities using time-dependent density functional quadratic response theory with local contribution and natural bond orbital analysis

We apply time-dependent density-functional quadratic response theory to investigate the static and dynamic second-order polarizabilities (first hyperpolarizability) beta. A new implementation using Slater-type basis functions, numerical integration, and density fitting techniques is reported. The second order coupled perturbed Kohn-Sham equations are solved and the second-order perturbed charge density is obtained. It is useful to highlight atomic and bond contributions to understand the relation between molecular structure and properties. Four moderately sized molecules (para-nitroaniline and derivatives thereof) are investigated to assess the accuracy of the time-dependent density-functional theory computations and to investigate the distribution of the second-order charge density as well as the "beta density." Our results highlight the contributions from atoms and bonds on different functional groups to the total value of beta with Mulliken-type and natural bond orbital (NBO) analyses, and demonstrate in some cases how contributions from a particular bond may be identified easily by visual inspection of the beta density. In addition, the position of side group substitution on carbon-carbon bonds significantly affects the hyperpolarizability. A contribution analysis as performed here might be helpful for the design of new materials with desired properties.     

The permutational symmetry in matrix multiplications

Summary This article presents slight modifications to algorithms forin-core symmetric matrix multiplications in order to optimize the computational number of multiplications required. The use of thepetite liste (Pl) algorithm, a general procedure of treating spatial symmetry in molecular calculations, is extended to permutational symmetry in matrix multiplications. This implementation requires the same number of operations as a regular matrix multiplication when the dimensions of original and transformed matrices are the same. However, when the transformed space dimension is smaller, this algorithm provides savings of up to a factor of two in the overall number of multiplications involved. Such a method can be viewed as an alternative demonstration to Saunders and van Lenthe's two-index transformation technique, who developed similar expressions through the decomposition of the symmetric matrix into its upper and lower triangular parts. The final equations obtained by these authors are the same as the ones shown here. However, the present method is supported by a solid theoretical framework, permutational group theory, which makes it general and applicable over any permutational symmetry available.     

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.     

A general,recursive, and open‐ended response code

We present a new implementation of a recent open‐ended response theory formulation for time‐ and perturbation‐dependent basis sets (Thorvaldsen et al., J. Chem. Phys. 2008, 129, 214108) at the Hartree–Fock and density functional levels of theory. A novel feature of the new implementation is the use of recursive programming techniques, making it possible to write highly compact code for the analytic calculation of any response property at any valid choice of rule for the order of perturbation at which to include perturbed density matrices. The formalism is expressed in terms of the density matrix in the atomic orbital basis, allowing the recursive scheme presented here to be used in linear‐scaling formulations of response theory as well as with two‐ and four‐component relativistic wave functions. To demonstrate the new code, we present calculations of the third geometrical derivatives of the frequency‐dependent second hyperpolarizability for HSOH at the Hartree–Fock level of theory, a seventh‐order energy derivative involving basis sets that are both time and perturbation dependent. © 2014 Wiley Periodicals, Inc.     

Coupled multiconfiguration self-consistent field (MC SCF) perturbation theory

The analytical form of the perturbation theory for the MC SCF method of Veillard and Clementi is presented. The appropriate second-order energy functional which takes into account the self-consistency requirements, leads to a set of coupled first-order perturbed equations determining the perturbed configuration coefficients and orbitals. The second-order energy formula derived from this functional can be given a clear physical interpretation. The present analytical approach is compared with the finite perturbation MC SCF scheme. The possibility of the approximate solution of the coupled MC SCF perturbation equations is also discussed and the so-called uncoupled procedures are devised. In the limit of the single determinant wave function the present formulae are shown to be equivalent to the appropriate Hartree-Fock perturbation results. The differences between the one-configuration SCF and the MC SCF approach are illustrated by the calculation of the electric dipole polarizability of. HZ in the CNDO/2 approximation. It is shown that the one-configuration SCF approaches cannot account for the correct asymptotic properties of the second-order energy for large internuclear distances. This feature of the SCF perturbation theories does not depend on the specific approximations of the CNDO/2 scheme and is corrected by using the MC SCF perturbation theory.     

Efficient linear-scaling calculation of response properties: density matrix-based Laplace-transformed coupled-perturbed self-consistent field theory

A density matrix-based Laplace reformulation of coupled-perturbed self-consistent field (CPSCF) theory is presented. It allows a direct, instead of iterative, solution for the integral-independent part of the density matrix-based CPSCF (D-CPSCF) equations [J. Kussmann and C. Ochsenfeld, J. Chem. Phys. 127, 054103 (2007)]. In this way, the matrix-multiplication overhead compared to molecular orbital-based solutions is reduced to a minimum, while at the same time, the linear-scaling behavior of D-CPSCF theory is preserved. The present Laplace-based equation solver is expected to be of general applicability.     

Location of saddle points and minimum energy paths by a constrained simplex optimization procedure

Two methods are proposed, one for the location of saddle points and one for the calculation of steepest-descent paths on multidimensional surfaces. Both methods are based on a constrained simplex optimization technique that avoids the evaluation of gradients or second derivative matrices. Three chemical reactions of increasing structural complexity are studied within the PRDDO SCF approximation. Predicted properties of reaction hypersurfaces are in good overall agreement with those determined by gradient minimization and gradient following algorithms in connection with various ab initio SCF methods. Computational efforts required by the new procedures are discussed.     

Direct minimization of the energy functional in LCAO-MO density matrix formalism I. Closed and open shell systems

A direct method of minimization of the energy expression for closed and open shell systems in LCAO-MO density matrix formalism is presented. The method makes use of a unitary transformation acting directly on the density matrices. Expressions of the gradient and second energy derivatives are worked out. Some preliminary calculations to test the rate of minimization using a variable metric method have been made on H₂S and SO molecules and have given satisfactory results.[/p]     

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.     

Perturbative virtual SCF CI treatment for energy levels of coupled oscillator systems

A perturbative SCF CI treatment to obtain energy levels of coupled oscillator systems is proposed. The method uses the virtual SCF basis set, and the SCF equations are solved by means of a perturbative treatment that provides the diagonal matrix elements involved in the CI calculation. The off-diagonal matrix elements are calculated using a commutation relationship derived from exact quantum theorems. Numerical results for several systems are obtained and compared with those from others SCF, SCF CI , and variational treatments.     

Probabilistic evolution theory for the solution of explicit autonomous ordinary differential equations: squarified telescope matrices

Probabilistic evolution theory (PREVTH) is used for the solution of initial value problems of first order explicit autonomous ordinary differential equation sets with second degree multinomial right hand side functions. It is an approximation method based on Kronecker power series: a rewriting of multivariate Taylor series using matrices having certain flexible parameters. Kronecker power series have matrices which are called telescope matrices: \(n \times n^{j+1}\) matrices where j is the index of summation. The additive terms of each telescope matrix is formed through Kronecker product from both sides by Kronecker powers of identity matrices. Recently, squarification is proposed in order to avoid the growing of the matrices in size at each additive term of the series. This paper explains the squarification procedure: the procedure used in order to avoid Kronecker multiplications within PREVTH so that the sizes of the matrices do not grow and so that the amount of necessary computation is reduced. The recursion between squarified matrices is also given. As a numerical application, the solution of a Hénon–Heiles system is provided.     

The trust-region self-consistent field method: towards a black-box optimization in Hartree-Fock and Kohn-Sham theories

The trust-region self-consistent field (TRSCF) method is presented for optimizing the total energy E(SCF) of Hartree-Fock theory and Kohn-Sham density-functional theory. In the TRSCF method, both the Fock/Kohn-Sham matrix diagonalization step to obtain a new density matrix and the step to determine the optimal density matrix in the subspace of the density matrices of the preceding diagonalization steps have been improved. The improvements follow from the recognition that local models to E(SCF) may be introduced by carrying out a Taylor expansion of the energy about the current density matrix. At the point of expansion, the local models have the same gradient as E(SCF) but only an approximate Hessian. The local models are therefore valid only in a restricted region-the trust region-and steps can only be taken with confidence within this region. By restricting the steps of the TRSCF model to be inside the trust region, a monotonic and significant reduction of the total energy is ensured in each iteration of the TRSCF method. Examples are given where the TRSCF method converges monotonically and smoothly, but where the standard DIIS method diverges.     

Bond fukui indices: Comparison of frozen molecular orbital and finite differences through mulliken populations

Bond Fukui functions and matrices are introduced for ab initio levels of theory using a Mulliken atoms in molecules model. It is shown how these indices may be obtained from first‐order density matrix derivatives without need for going to second‐order density matrices as in a previous work. The importance of taking into account the nonorthogonality of the basis in ab initio calculations is shown, contrasting the present results with previous work based on Hückel theory. It is shown how the extension of Fukui functions to Fukui matrices allows getting more insight into the nature of bond Fukui functions. All presently introduced indices respect the necessary normalization conditions and include the classical single atom condensed Fukui functions. © 2013 Wiley Periodicals, Inc.