Analytical Energy Gradients Research Articles

Molecular Structure Theory without the Born-Oppenheimer Approximation: Rotationless Vibrational States of LiH.

Low-lying rotationless states of the lithium hydride molecule are studied in the framework of the variational method without assuming the Born-Oppenheimer (BO) approximation. Highly accurate solutions to the six-particle (two nuclei and four electrons) Schrödinger equation are obtained by means of expanding the wave functions of the considered states in terms of many thousands of all-particle explicitly correlated Gausssians. The basis functions are optimized independently for each state using the analytic energy gradient with respect to the nonlinear parameters. The non-BO wave functions obtained in the calculations are used to evaluate the leading-order relativistic and quantum electrodynamics energy corrections in the framework of the perturbation theory. The geometric structure of the molecule in the ground and excited states is discussed based on the analysis of the nucleus-nucleus correlation functions. The non-BO energies and structural parameters obtained of this work are compared with the most accurate BO results currently available.

Analytic Gradients for Equation-of-Motion Coupled Cluster with Single, Double, and Perturbative Triple Excitations.

Understanding the process of molecular photoexcitation is crucial in various fields, including drug development, materials science, photovoltaics, and more. The electronic vertical excitation energy is a critical property, for example in determining the singlet-triplet gap of chromophores. However, a full understanding of excited-state processes requires additional explorations of the excited-state potential energy surface and electronic properties, which is greatly aided by the availability of analytic energy gradients. Owing to its robust high accuracy over a wide range of chemical problems, equation-of-motion coupled cluster with single and double excitations (EOM-CCSD) is a powerful method for predicting excited-state properties, and the implementation of analytic gradients of many EOM-CCSD variants (excitation energies, ionization potentials, electron attachment energies, etc.) along with numerous successful applications highlights the flexibility of the method. In specific cases where a higher level of accuracy is needed or in more complex electronic structures, the inclusion of triple excitations becomes essential, for example, in the EOM-CCSD* approach of Saeh and Stanton. In this work, we derive and implement for the first time the analytic gradients of EOMEE-CCSD*, which also provides a template for analytic gradients of related excited-state methods with perturbative triple excitations. The capabilities of analytic EOMEE-CCSD* gradients are illustrated by several representative examples.

Analytical Gradient Using Cluster-in-Molecule RI-MP2 Method for the Geometry Optimizations of Large Systems.

We present an efficient analytical energy gradient algorithm for the cluster-in-molecule resolution-of-identity second-order Møller-Plesset perturbation (CIM-RI-MP2) method based on the Lagrange multiplier method. Our algorithm independently constructs the Lagrangian formalism within each cluster, avoiding the solution of the coupled-perturbed Hartree-Fock (CPHF) equation for the whole system. Due to this feature, the computational cost of the CIM-RI-MP2 gradients is much lower than that of other local MP2 algorithms. Benchmark calculations of several molecules containing up to 312 atoms demonstrate the general applicability of our CIM-RI-MP2 gradient algorithm. The optimized structure of a 244-atom molecule using the CIM-RI-MP2 method with the cc-pVDZ basis set is in good agreement with the corresponding crystal structure. A single-point gradient calculation conducted for a molecular cage containing 972 atoms and 9612 basis functions takes 48 h on 25 nodes, utilizing a total of 600 CPU cores. The present CIM-RI-MP2 gradient program is applicable for obtaining the optimized geometries of large systems with hundreds of atoms.

High-Performance Multi-GPU Analytic RI-MP2 Energy Gradients.

This article presents a novel algorithm for the calculation of analytic energy gradients from second-order Møller-Plesset perturbation theory within the Resolution-of-the-Identity approximation (RI-MP2), which is designed to achieve high performance on clusters with multiple graphical processing units (GPUs). The algorithm uses GPUs for all major steps of the calculation, including integral generation, formation of all required intermediate tensors, solution of the Z-vector equation and gradient accumulation. The implementation in the EXtreme Scale Electronic Structure System (EXESS) software package includes a tailored, highly efficient, multistream scheduling system to hide CPU-GPU data transfer latencies and allows nodes with 8 A100 GPUs to operate at over 80% of theoretical peak floating-point performance. Comparative performance analysis shows a significant reduction in computational time relative to traditional multicore CPU-based methods, with our approach achieving up to a 95-fold speedup over the single-node performance of established software such as Q-Chem and ORCA. Additionally, we demonstrate that pairing our implementation with the molecular fragmentation framework in EXESS can drastically lower the computational scaling of RI-MP2 gradient calculations from quintic to subquadratic, enabling further substantial savings in runtime while retaining high numerical accuracy in the resulting gradients.

Efficient analytical gradients of property-based diabatic states: Geometry optimizations for localized holes.

We present efficient analytical gradients of property-based diabatic states and couplings using a Lagrangian formalism. Unlike previous formulations, the method achieves a computational scaling that is independent of the number of adiabatic states used to construct the diabats. The approach is generalizable to other property-based diabatization schemes and electronic structure methods as long as analytical energy gradients are available and integral derivatives with the property operator can be formed. We also introduce a scheme to phase and reorder diabats to ensure their continuity between molecular configurations. We demonstrate this for the specific case of Boys diabatic states obtained from state-averaged complete active space self-consistent field electronic structure calculations with GPU acceleration in the TeraChem package. The method is used to test the Condon approximation for the hole transfer in an explicitly solvated model DNA oligomer.

Analytic Gradient for Time-Dependent Density Functional Theory Combined with the Fragment Molecular Orbital Method.

The analytic energy gradient of energy with respect to nuclear coordinates is derived for the fragment molecular orbital (FMO) method combined with time-dependent density functional theory (TDDFT). The response terms arising from the use of a polarizable embedding are derived. The obtained analytic FMO-TDDFT gradient is shown to be accurate in comparison to both numerical FMO-TDDFT and unfragmented TDDFT gradients, at the level of two- and three-body expansions. The gradients are used for geometry optimizations, molecular dynamics, vibrational calculations, and simulations of IR and Raman spectra of excited states. The developed method is used to optimize the geometry of the ground and excited electronic states of the photoactive yellow protein (PDB: 2PHY).

Algorithm for analytic nuclear energy gradients of state averaged DMRG-CASSCF theory with newly derived coupled-perturbed equations.

We present an algorithm for evaluating analytic nuclear energy gradients of the state-averaged density matrix renormalization group complete-active-space self-consistent field (SA-DMRG-CASSCF) theory based on the newly derived coupled-perturbed (CP) DMRG-CASSCF equations. The Lagrangian for the conventional SA-CASSCF analytic gradient theory is extended to the SA-DMRG-CASSCF variant that can fully consider a whole set of constraints on the parameters of multi-root canonical matrix product states formed at all the DMRG block configurations. An efficient algorithm to solve the CP-DMRG-CASSCF equations for determining the multipliers was developed. The complexity of the resultant analytic gradient algorithm is overall the same as that of the unperturbed SA-DMRG-CASSCF algorithm. In addition, a reduced-scaling approach was developed to directly compute the SA reduced density matrices (SA-RDMs) and their perturbed ones without calculating separate state-specific RDMs. As part of our implementation scheme, we neglect the term associated with the constraint on the active orbitals in terms of the active-active rotation in the Lagrangian. Thus, errors from the true analytic gradients may be caused in this scheme. The proposed gradient algorithm was tested with the spin-adapted implementation by checking how accurately the computed analytic energy gradients reproduce numerical gradients of the SA-DMRG-CASSCF energies using a common number of renormalized bases. The illustrative applications show that the errors are sufficiently small when using a typical number of the renormalized bases, which is required to attain adequate accuracy in DMRG's total energies.

Automated handling of complex chemical structures in Z-matrix coordinates-The chemcoord library.

In this work, we present a fully automated method for the construction of chemically meaningful sets of hierarchical nonredundant internal coordinates (ICs; also commonly denoted as Z-matrices) from the Cartesian coordinates of a molecular system. Particular focus is placed on avoiding ill-definitions of angles and dihedrals due to linear arrangements of atoms, to consistently guarantee a well-defined transformation to Cartesian coordinates, even after structural changes. The representations thus obtained are particularly well suited for pathway construction in double-ended methods for transition state search and optimizations with nonlinear constraints. Analytical gradients for the transformation between the coordinate systems were derived for analytical geometry optimizations purely in Z-matrix coordinates. The geometry optimization was coupled with a Symbolic Algebra package to support arbitrary nonlinear constraints in Z-matrix coordinates, while retaining analytical energy gradient conversion. The difference to the commonly used nonhierarchical IC transformations is discussed. Sample applications are provided for a number of common chemical reactions and illustrative examples.

Oscillator strengths and interstate transition energies involving [formula omitted] and [formula omitted] states of the Li atom

We report high accuracy calculations of the ground and excited doublet S and P states of lithium atom. Overall, 24 states corresponding to dominant electronic configurations 1s2ns and 1s2np (n=2,…,13) are considered in the framework of the Ritz variational method. The nonrelativistic wave function of each of these states is generated in an independent calculation by expanding it in terms of a large number (K=11,000−17,000) of all-electron explicitly correlated Gaussian functions (ECG) whose nonlinear parameters are extensively optimized with a procedure that employs analytic energy gradient determined with respect to these parameters. The Hamiltonian used in the calculations explicitly depends on the mass of the nucleus. The leading relativistic and quantum electrodynamics (QED) corrections to the energy levels are subsequently computed using the perturbation-theory approach and the variational nonrelativistic wave functions as the zeroth-order functions. As these functions are generated in the finite-nuclear-mass (FNM) calculations, the energy corrections include the nuclear recoil effects. The obtained energy levels allow us to determine highly accurate interstate transition frequencies for both the naturally occurring stable lithium isotopes, 6Li and 7Li, and the lithium atom with an infinitely heavy nucleus, ∞Li. The nonrelativistic wave functions are used to compute the transition dipole moments and the corresponding oscillator strengths. These quantities are reported for 144 S–P transitions of each isotope. The data set generated in this work is considerably more accurate and comprehensive than the data available from the previous theoretical calculations. It can be useful in guiding future spectroscopic measurements of the lithium atom.

Solvent effects in four-component relativistic electronic structure theory based on the reference interaction-site model.

A combined method of the Dirac-Hartree-Fock (DHF) method and the reference interaction-site model (RISM) theory is reported; this is the initial implementation of the coupling of the four-component relativistic electronic structure theory and an integral equation theory of molecular liquids. In the method, the DHF and RISM equations are solved self-consistently, and therefore the electronic structure of the solute, including relativistic effects, and the solvation structure are determined simultaneously. The formulation is constructed based on the variational principle with respect to the Helmholtz energy, and analytic free energy gradients are also derived using the variational property. The method is applied to the iodine ion (I- ), methyl iodide (CH3 I), and hydrogen chalcogenide (H2 X, where X=O-Po) in aqueous solutions, and the electronic structures of the solutes, as well as the solvation free energies and their component analysis, solvent distributions, and solute-solvent interactions, are discussed.

Detection and Correction of Delocalization Errors for Electron and Hole Polarons Using Density-Corrected DFT.

Modeling polaron defects is an important aspect of computational materials science, but the description of unpaired spins in density functional theory (DFT) often suffers from delocalization error. To diagnose and correct the overdelocalization of spin defects, we report an implementation of density-corrected (DC-)DFT and its analytic energy gradient. In DC-DFT, an exchange-correlation functional is evaluated using a Hartree-Fock density, thus incorporating electron correlation while avoiding self-interaction error. Results for an electron polaron in models of titania and a hole polaron in Al-doped silica demonstrate that geometry optimization with semilocal functionals drives significant structural distortion, including the elongation of several bonds, such that subsequent single-point calculations with hybrid functionals fail to afford a localized defect even in cases where geometry optimization with the hybrid functional does localize the polaron. This has significant implications for traditional workflows in computational materials science, where semilocal functionals are often used for structure relaxation. DC-DFT calculations provide a mechanism to detect situations where delocalization error is likely to affect the results.

Cavity quantum-electrodynamical time-dependent density functional theory within Gaussian atomic basis. II. Analytic energy gradient.

Following the formulation of cavity quantum-electrodynamical time-dependent density functional theory (cQED-TDDFT) models [Flick et al., ACS Photonics 6, 2757-2778 (2019) and Yang et al., J. Chem. Phys. 155, 064107 (2021)], here, we report the derivation and implementation of the analytic energy gradient for polaritonic states of a single photochrome within the cQED-TDDFT models. Such gradient evaluation is also applicable to a complex of explicitly specified photochromes or, with proper scaling, a set of parallel-oriented, identical-geometry, and non-interacting molecules in the microcavity.

Analytical Energy Gradient for State-Averaged Orbital-Optimized Variational Quantum Eigensolvers and Its Application to a Photochemical Reaction.

Elucidating photochemical reactions is vital to understanding various biochemical phenomena and developing functional materials such as artificial photosynthesis and organic solar cells, albeit with notorious difficulty in both experiments and theories. The best theoretical way so far to analyze photochemical reactions at the level of ab initio electronic structure is the state-averaged multiconfigurational self-consistent field (SA-MCSCF) method. However, the exponential computational cost of classical computers with the increasing number of molecular orbitals hinders applications of SA-MCSCF for large systems we are interested in. Utilizing quantum computers was recently proposed as a promising approach to overcome such computational cost, dubbed as state-averaged orbital-optimized variational quantum eigensolver (SA-OO-VQE). Here, we extend a theory of SA-OO-VQE so that analytical gradients of energy can be evaluated by standard techniques that are feasible with near-term quantum computers. The analytical gradients, known only for the state-specific OO-VQE in previous studies, allow us to determine various characteristics of photochemical reactions such as the conical intersection (CI) points. We perform a proof-of-principle calculation of our methods by applying it to the photochemical cis-trans isomerization of 1,3,3,3-tetrafluoropropene. Numerical simulations of quantum circuits and measurements can correctly capture the photochemical reaction pathway of this model system, including the CI points. Our results illustrate the possibility of leveraging quantum computers for studying photochemical reactions.

Analytic Energy Gradients for the Driven Similarity Renormalization Group Multireference Second-Order Perturbation Theory.

We derive analytic energy gradients of the driven similarity renormalization group (DSRG) multireference second-order perturbation theory (MRPT2) using the method of Lagrange multipliers. In the Lagrangian, we impose constraints for a complete-active-space self-consistent-field reference wave function and the semicanonical orthonormal molecular orbitals. Solving for the associated Lagrange multipliers is found to share the same asymptotic scaling of a single DSRG-MRPT2 energy computation. A pilot implementation of the DSRG-MRPT2 analytic gradients is used to optimize the geometry of the singlet and triplet states of p-benzyne. The equilibrium bond lengths and angles are similar to those computed via other MRPT2s and Mukherjee's multireference coupled cluster theory. An approximate DSRG-MRPT2 method that neglects the contributions of the three-body density cumulant is found to introduce negligible errors in the geometry of p-benzyne, lending itself to a promising low-cost approach for molecular geometry optimizations using large active spaces.

Third-Order Many-Body Expansion of OSV-MP2 Wave Function for Low-Order Scaling Analytical Gradient Computation.

We present a many-body expansion (MBE) formulation and implementation for efficient computation of analytical energy gradients from the orbital-specific-virtual second-order Møllet-Plesset perturbation theory (OSV-MP2) based on our earlier work (Zhou et al. J. Chem. Theory Comput. 2020, 16, 196-210). The third-order MBE(3) expansion of OSV-MP2 amplitudes and density matrices was developed to adopt the orbital-specific clustering and long-range termination schemes, which avoids term-by-term differentiations of the MBE energy bodies. We achieve better efficiency by exploiting the algorithmic sparsity that allows us to prune out insignificant fitting integrals and OSV relaxations. With these approximations, the present implementation is benchmarked on a range of molecules that show an economic scaling in the linear and quadratic regimes for computing MBE(3)-OSV-MP2 amplitude and gradient equations, respectively, and yields normal accuracy comparable to the original OSV-MP2 results. The MPI-3-based parallelism through shared memory one-sided communication is further developed for improving parallel scalability and memory accessibility by sorting the MBE(3) orbital clusters into independent tasks that are distributed on multiple processes across many nodes, supporting both global and local data locations in which selected MBE(3)-OSV-MP2 intermediates of different sizes are distinguished and accordingly placed. The accuracy and efficiency level of our MBE(3)-OSV-MP2 analytical gradient implementation is finally illustrated in two applications: we show that the subtle coordination structure differences of mechanically interlocked Cu-catenane complexes can be distinguished when tuning ligand lengths; and the porphycene molecular dynamics reveals the emergence of the vibrational signature arising from softened N-H stretching associated with hydrogen transfer, using an MP2 level of electron correlation and classical nuclei for the first time.

Energy correction and analytic energy gradients due to triples in CCSD(T) with spin–orbit coupling on graphic processing units using single-precision data

ABSTRACT Calculating the contribution of triples ((T)) to the correlation energy, the density matrices and the constant terms in the Λ equation is the most expensive steps in obtaining analytic energy gradients for the CCSD(T) approach. In this work, we report the implementation of these steps for the CCSD(T) method with spin–orbit coupling (SOC) included in post-self-consistent-field calculations (SOC-CCSD(T)) using single-precision data on a consumer GPU card to accelerate calculations. The developed program can be used for calculations on GPU with single-precision data or on CPU with either single- or double-precision data. According to our results, calculating the (T) correlation energy in SOC-CCSD(T) on GPU with single-precision data is about 7–10 times faster for the investigated molecules than that on CPU with double-precision data, and it is 5–9 times faster in calculating the (T) part of analytic energy gradients for the SOC-CCSD(T) method. Our results indicate that loss of accuracy for energy gradients and equilibrium structures using single-precision data is negligible, while a mixed-precision calculation is needed for the (T) correlation energy. In calculating harmonic frequencies based on finite difference of analytic energy gradients, a larger step size is required to achieve accurate results for certain vibrational modes with single-precision data.

NAC-TDDFT: Time-Dependent Density Functional Theory for Nonadiabatic Couplings

First-order nonadiabatic coupling (NAC) matrix elements (fo-NACMEs) are the basic quantities in theoretical descriptions of electronically nonadiabatic processes that are ubiquitous in molecular physics and chemistry. Given the large size of systems of chemical interests, time-dependent density functional theory (TDDFT) is usually the first choice of methods. However, the lack of many-electron wave functions in TDDFT renders the formulation of NAC-TDDFT for fo-NACMEs conceptually difficult. Because of this, various variants of NAC-TDDFT have been proposed in the literature from different standing points, including the Hellmann-Feynman-like expression and auxiliary/pseudo-wave function (AWF)-, equation-of-motion (EOM)-, and time-dependent perturbation theory (TDPT)-based formulations. Based on critical analyses, the following conclusions are made here: (1) The Hellmann-Feynman-like expression, which is rooted in exact wave function theory, is hardly useful due to huge demand on basis sets. (2) Although most popular, the AWF variants of NAC-TDDFT are not theoretically founded and become ambiguous particularly for the fo-NACMEs between two excited states, although they do agree with the EOM and TDPT variants under the Tamm-Dancoff approximation. (3) The TDPT variant of NAC-TDDFT is theoretically most rigorous but suffers from numerical instabilities on the one hand and does not differ to a significant extent from the EOM variant on the other hand. (4) As such, the EOM variant of NAC-TDDFT for the fo-NACMEs between the ground and excited states and between two excited states is solely the right choice in practice. These formal analyses are fully supported by numerical experimentations, taking azulene as a showcase. The proper implementation of the EOM variant of NAC-TDDFT is also highlighted, showing that the fo-NACMEs between the ground and excited states and between two excited states are computationally very much the same as the analytic energy gradients of DFT and TDDFT, respectively. Possible future developments of the EOM variant of NAC-TDDFT are also highlighted. Its extensions to spin-adapted open-shell TDDFT and proper treatment of spin-orbit couplings (which are another source of force for electronically nonadiabatic processes) are particularly warranted in the near future.

Benchmark calculations of the 2 D Rydberg spectrum of lithium

The recently implemented algorithm for calculating the leading relativistic corrections for D states of an atom with an arbitrary number of electrons is used to perform calculations of the lowest eleven Rydberg states of lithium. The corrections are obtained as expectation values of the corresponding relativistic operators with the non-relativistic wave functions expanded in terms of all-electron explicitly correlated Gaussian (ECG) functions. The wave functions are generated using the standard variational method and the finite-nuclear-mass (FNM) approach. The non-linear parameters of the ECGs are variationally optimised using a algorithm that employs the analytical energy gradient determined with respect to these parameters. Infinite-nuclear-mass calculations are also performed. The calculated FNM state energies for Li and Li are used to calculate isotopic shifts of the inter-state transition energies. The calculated interstate transitions are in good agreement with the high-quality experimental results.

High-accuracy calculations of the lowest eleven Rydberg 2 P states of the Li atom

Highly accurate calculations are reported for the eleven lowest states of the 2 P Rydberg series (1s 2 np 1, n = 2, …, 12) of the lithium atom. The nonrelativistic wave functions of the states are expanded in terms of up to 16 000 all-electron explicitly correlated Gaussian (ECG) basis functions. The ECG exponential parameters are variationally optimized using a method that employs the analytical energy gradient determined for the parameters. The finite-nuclear-mass effects of the 6Li and 7Li isotopes are explicitly included in the nonrelativistic variational calculations. The results also include the leading relativistic and quantum electrodynamics energy corrections computed using the framework of perturbation theory. The calculated interstate transition energies are compared with the available experimental data. The 6Li–7Li isotope shifts of the transition energies are determined.

Analytical energy gradient for the second-order Møller-Plesset perturbation theory coupled with the reference interaction site model self-consistent field explicitly including spatial electron density distribution.

Solvatochromic shifts of the activation free energies are important aspects to consider for reaction control. To predict the energies, the stationary points in a solution must be accurately determined along the reaction pathway. In this study, the second-order Møller-Plesset perturbation (MP2) theory combined with the reference interaction site model was applied using our fitting approach, and the MP2 analytical energy gradient was determined. The coupled-cluster energy and thermal correction were calculated using the MP2 optimized geometry with solvent effect, and the activation free energies of the Diels-Alder reaction between cyclopentadiene and methyl vinyl ketone are within an error of 2 kcal/mol compared with the experimental data.