Performant Automatic Differentiation of Local Coupled Cluster Theories: Response Properties and Ab Initio Molecular Dynamics

Xing Zhang Chenghan Li Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA Hong-Zhou Ye Timothy C. Berkelbach Department of Chemistry, Columbia University, New York, NY 10027, USA Garnet Kin-Lic Chan Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USA gkc1000@gmail.com

Abstract

In this work, we introduce a differentiable implementation of the local natural orbital coupled cluster (LNO-CC) method within the automatic differentiation framework of the PySCFAD package. The implementation is comprehensively tuned for enhanced performance, which enables the calculation of first-order static response properties on medium-sized molecular systems using coupled cluster theory with single, double, and perturbative triple excitations [CCSD(T)]. We evaluate the accuracy of our method by benchmarking it against the canonical CCSD(T) reference for nuclear gradients, dipole moments, and geometry optimizations. In addition, we demonstrate the possibility of property calculations for chemically interesting systems through the computation of bond orders and Mössbauer spectroscopy parameters for a [NiFe]-hydrogenase active site model, along with the simulation of infrared (IR) spectra via ab initio LNO-CC molecular dynamics for a protonated water hexamer.

I Introduction

Since the pioneering work of Pulay and Sæbø, [1, 2, 3, 4, 5] local electron correlation methods have seen significant advances over the past two decades. [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45] With modern developments, calculations employing local approximations to coupled cluster theory with single, double, and perturbative triple excitations[46] [CCSD(T)] can now be routinely performed for large molecular systems [29, 36, 38, 39, 40, 42] and solids.[43, 45, 47] However, the majority of these calculations are focused on ground-state electronic energetics, while the application of local correlation methods to molecular properties remains relatively underexplored. This may be attributed to the inherent complexity of these methods, which makes the implementation of analytic derivatives or response theory more challenging, as compared to their canonical counterparts. Nevertheless, notable efforts have been made on this subject. Analytic nuclear gradients for projected atomic orbital (PAO) based local correlation methods,[6, 7, 8] including local second-order Møller–Plesset perturbation theory (LMP2) and local coupled cluster theory with single and double excitations (LCCSD), have been developed by Werner and co-workers. [48, 49, 50, 51] Meanwhile, Schütz and co-workers have implemented analytic nuclear gradients and dipole moments for the local CC2 method, applied to both ground and excited states. [52, 53] In addition, calculations of nuclear magnetic resonance (NMR) shieldings[54, 55] and magnetizabilities[56] have been reported at the LMP2 level. More recently, Neese and co-workers extended their domain-based local pair natural orbital (DLPNO) approaches[14, 15] for computing static response properties. These include first[57, 58] and second[59] derivatives of DLPNO-MP2 and orbital-unrelaxed first derivatives of DLPNO-CCSD.[60] Similarly, Yang and co-workers published analytic nuclear gradient implementations for the orbital-specific virtual (OSV) MP2 method.[61] Finally, Crawford and co-workers studied dynamic (frequency dependent) response properties, such as polarizabilities and optical rotations, by applying the local coupled cluster linear response theory. [62, 63, 64, 65, 66]

Typically, static response properties are calculated as energy (or Lagrangian) derivatives with respect to the perturbations. Computing analytic derivatives simply involves applying chain rules to the objective function. However, manually tracking the entire workflow of a complex calculation can be tedious and error-prone. Thankfully, modern automatic differentiation (AD) tools[67, 68, 69] offload the task of tracing the program’s execution to the computer,[70, 71] greatly simplifying the implementation of analytic derivatives. In recent years, there has been a growing effort to integrate AD techniques into quantum chemistry computations. [72, 73, 74, 75, 76, 77, 78, 79] Among these, the PySCFAD package is distinguished by its extensive functionality and its flexibility.[76] In an earlier publication,[76] we demonstrated a proof-of-concept showing that PySCFAD could be utilized for the rapid prototyping of new methodologies. Here, we further illustrate that PySCFAD is now also effective for production-level calculations.

In this work, we present a differentiable implementation of the local natural orbital coupled cluster (LNO-CC) method, first introduced by Rolik and Kállay,[17] within the AD framework of PySCFAD. The modular design of PySCFAD allows for the integration of new methods with virtually no changes to the existing components of the software. Furthermore, performance optimizations can be targeted to a small segment of the program, e.g., the tensor contractions for computing the triple excitation correction in the CCSD(T) method, without compromising the differentiability of the entire computational workflow. We demonstrate the resulting differentiable LNO-CC method in some non-trivial applications: calculating the bond orders and Mössbauer spectroscopy parameters for a model system of the [NiFe]-hydrogenase active site, as well as obtaining the anharmonic infrared (IR) spectra of a protonated water hexamer from ab initio molecular dynamics.

II Notations

Throughout the paper, we use Greek letters ( $\mu$ , $\nu$ , …) to denote atomic orbitals. Occupied and virtual canonical Hartree-Fock (HF) orbitals are labeled using lowercase letters ( $i$ , $j$ , … for occupied; $a$ , $b$ , … for virtual), while uppercase letters ( $I$ , $J$ , …) signify local orbitals, whether semi-canonicalized or not. The energies corresponding to canonical HF orbitals and semi-canonical local orbitals are denoted by $\varepsilon$ and $\mathcal{E}$ , respectively. Orbitals that span a local active space are distinguished by tildes ( $\tilde{i}$ , $\tilde{j}$ , …). Unless specified otherwise, equations are expressed in the spin orbital formalism, although our implementation is spin-adapted, and we adopt Dirac’s notation for the electron integrals.

III Methodology

Our differentiable LNO-CC method is built upon the recent implementation of Ye and Berkelbach,[45] with the latter limited to calculating ground-state energies. In this section, we elaborate on this realization of the LNO-CC method within the PySCFAD framework, highlighting adjustments made for efficient differentiation, and also discuss the potential issues that arise during the computation of response properties.

The procedure starts by constructing a set of occupied orthonormal local orbitals. In the original scheme by Rolik and Kállay,[17] each local orbital along with the corresponding local natural orbitals (LNOs, see the definition below) defines a local active space, within which the local electron correlation calculation is carried out. We name this the one-orbital scheme. However, the size of the computational graph grows linearly with the number of local electron correlation calculations. For systems containing hundreds of electrons, it becomes time-consuming to trace the computation and to compile the program in a just-in-time[70] (JIT) fashion. A workaround, which we call the multi-orbital scheme, is to group the local orbitals to form fragments, similar to the strategy employed by the fragment molecular orbital (FMO) approach.[80] In particular, we designate each heavy atom and its surrounding hydrogen atoms as a fragment, and the local orbitals are assigned to the corresponding fragments based on a Löwdin population analysis. The resulting fragments are sufficiently small to allow for efficient high-level local electron correlation calculations [e.g., at the CCSD(T) level], while maintaining a moderate number of local calculations, facilitating tractable computation tracing, JIT compilation, and gradient backpropagation.

With the local orbitals and fragments, we determine the local active space as follows. Suppose $|\phi_{I}^{\Omega}\rangle$ are the semi-canonical local orbitals on a fragment $\Omega$ , then the corresponding local active space is spanned by $|\phi_{I}^{\Omega}\rangle$ plus a set of LNOs determined from diagonalizing the fragment contribution to the MP2 density matrix. Following Ref. 17, the occupied-occupied (OO) and virtual-virtual (VV) blocks of this density matrix read as

D_{jk}^{\Omega}=\sum_{mn}\mathcal{P}^{\Omega\top}_{jm}\left(\frac{1}{2}\sum_{I\in\Omega}\sum_{ab}t_{Imab}t_{Inab}\right)\mathcal{P}^{\Omega}_{nk}\;,

(1)

and

D_{ab}^{\Omega}=\frac{1}{2}\sum_{I\in\Omega}\sum_{jc}t_{Ijac}t_{Ijbc}\;,

(2)

respectively, where

t_{Ijab}=\frac{\langle ab||Ij\rangle}{\mathcal{E}_{I}+\varepsilon_{j}-\varepsilon_{a}-\varepsilon_{b}}\;,

(3)

and $\bm{\mathcal{P}}^{\Omega}$ is a projection matrix that removes the contributions of $|\phi_{I}^{\Omega}\rangle$ from the density matrix. Note that if $|\phi_{I}^{\Omega}\rangle$ overlaps with the virtual space, e.g., when intrinsic atomic orbitals (IAOs) are taken as the local orbitals, a similar projection is performed in Eq. 2 as well. Diagonalizing the OO (Eq. 1) and VV (Eq. 2) blocks of the density matrix gives the occupied and virtual LNOs associated with fragment $\Omega$ , respectively. In practice, we truncate the local active space by discarding those LNOs whose occupation numbers are smaller than a predefined threshold $\zeta$ . The retained LNOs along with $|\phi_{I}^{\Omega}\rangle$ are further semi-canonicalized to simplify the subsequent local electron correlation calculations on each fragment.

Taking MP2 as an example, the local correlation energy on a fragment can be expressed as

E^{\Omega}_{\text{MP2}}=\sum_{I\in\Omega}\sum_{\tilde{m}\tilde{n}\in\Omega}U^{\Omega\top}_{I\tilde{m}}\left(\frac{1}{4}\sum_{\tilde{j}\tilde{a}\tilde{b}\in\Omega}\langle\tilde{a}\tilde{b}||\tilde{m}\tilde{j}\rangle t_{\tilde{n}\tilde{j}\tilde{a}\tilde{b}}\right)U^{\Omega}_{\tilde{n}I}\;,

(4)

where $\mathbf{U}^{\Omega}$ transforms the semi-canonical orbitals that span a local active space to the corresponding local orbitals on fragment $\Omega$ . This orbital transformation restores the fragment-based energy partitioning. Therefore, the total correlation energy of the system can be computed by summing over the fragments

E_{\text{LNO-MP2}}=\sum_{\Omega}E^{\Omega}_{\text{MP2}}\;.

(5)

The energy expression for the LNO-CC method can be derived similarly, and is detailed in the supporting information (see Sec. LABEL:sec:cc_energy). One caveat is that the orbital transformation in Eq. 4 may break certain permutation symmetries of electron integrals and CC amplitudes. In particular, this happens when computing the triple excitation correction within the LNO-CCSD(T) method,[22] which increases the computational cost by a factor of at most three for closed-shell systems, compared to the canonical CCSD(T) calculation with the same correlation domain size. Finally, it is often beneficial to perform a global electron correlation calculation at a lower level (such as MP2 in this work) to correct for the correlation effects due to the weak pair interactions.[4] This results in the correlation energy expression for the LNO-CC method being

E_{\text{LNO-CC}}=\sum_{\Omega}\left(E^{\Omega}_{\text{CC}}-E^{\Omega}_{\text{MP2}}\right)+E_{\text{MP2}}\;.

(6)

In the course of implementing the LNO-CC method, several enhancements have been applied to the core components of PySCFAD, leading to improved performance. These include: (i) incorporating permutation symmetries for electron integrals and coupled cluster amplitudes, (ii) a manually optimized implementation for the gradient of tensor contractions associated with the triple excitation correction in CCSD(T), and (iii) minimizing the memory footprint for gradient calculations by recomputing intermediate quantities during backpropagation. (More details can be found in the supporting information Sec. LABEL:sec:opt.) The optimized PySCFAD exhibits efficiency on par with its parent program PySCF,[81] especially for MP2 and CC methods (see Fig. LABEL:fig:timing). Additionally, the LNO-CC calculations have been parallelized using the Message Passing Interface (MPI),[82] whereby the computations on distinct fragments are distributed across multiple processes. (The mean-field part of the calculation is duplicated within each process for straightforward backpropagation.)

Within the LNO-CC method, local correlation domains are truncated according to a simple cutoff for the LNOs, which neither yields a continuous energy function across the potential energy surface, nor preserves the molecular point-group symmetry. As such, nuclear gradients or any other response properties computed using the LNO-CC method may inherently contain errors stemming from the energy discontinuity or the symmetry breaking. Nonetheless, it is observed in practice that these errors tend to be small, provided that the correlation domains are properly converged with respect to the total energy (see next section).

Finally, there are situations where the evaluation of the orbital localization response becomes ill-defined.[58] This occurs when a continuum of solutions that fulfill the localization criteria exists.[83, 84] (In other words, the orbital rotation Hessian is rank deficient.) A comprehensive exploration of this problem falls beyond the scope of our current work. Instead, we offer a practical remedy to prevent singular gradients, which closely follows the strategy introduced in Ref. 58 (see Sec. LABEL:sec:local_orb_autodiff).

IV Results and Discussions

In this section, we present calculations of nuclear gradients, dipole moments, geometry optimizations, bond orders, quadrupole splittings (via electric field gradients), and ab initio molecular dynamics, using the LNO-CCSD(T) method. All calculations employ the density fitting[85, 86] approximation for the two electron repulsion integral.

IV.1 Nuclear Gradient and Dipole Moment

First, we examine the correctness and convergence with threshold of nuclear gradients and dipole moments computed using the LNO-CCSD(T) method. The Pipek-Mezey (PM) procedure[87] was employed to determine the local orbitals. The Baker test set, a set of 30 molecules with main-group elements of the first three rows ranging in size from 3 to 29 atoms,[88] is considered, and the reference data were obtained at the canonical CCSD(T) level of theory (where AD was employed to compute the nuclear gradients and the dipole moments). In Fig. 1, we plot, for each molecule in the test set, the absolute energy error, the nuclear gradient root mean square deviation (RMSD), the dipole moment RMSD, and the relative active space size. (The relative active space size is defined as the ratio of the mean number of orbitals within the local active space on each fragment to the total number of orbitals.)

It can be seen that, on average, the errors in nuclear gradients and dipole moments are of the same order as those in energies. Enlarging the basis set from double zeta to triple zeta does not lead to an increase in the errors (as shown in Fig. LABEL:fig:grad_pm_TZ). We note that both nuclear gradients and dipole moments are types of first-order response properties, which can be determined using the zeroth-order wavefunction (or density matrix), according to Wigner’s $2n+1$ rule. However, for higher order response properties, such as polarizabilities, which depend on the first-order wavefunction, it may be necessary to use larger correlation domains to achieve the desired accuracy, as suggested by Crawford and co-workers.[65]

There are also some outliers in Fig. 1, such as benzene (molecule 7), for which the large errors are likely due to symmetry breaking (note the non-zero dipole moment). For these systems, the coupled perturbed localization (CPL) equations (Eq. LABEL:eq:local_zvec) may not have solutions unless both the fragmentation and the truncation of LNOs are performed in ways that preserve the point-group symmetry. These kinds of errors, however, can be significantly reduced by using projection-based local orbitals, such as intrinsic atomic orbitals[89] (IAOs), which eliminate the need to solve the CPL equations. (See results in Fig. LABEL:fig:grad_iao.)

Finally, increasing the size of the local correlation domain generally improves the accuracy of energies, nuclear gradients, and dipole moments. Many of the molecules in this test set are small (10 of them have fewer than 10 atoms), but for the majority of systems assessed, an accuracy of $10^{-3}$ a.u. may be achieved with a local active space containing no more than $50\%$ of the total orbitals, with a smaller fraction of orbitals needed in the larger molecules.

Refer to caption — Figure 1: Absolute energy error, nuclear gradient RMSD, dipole moment RMSD, and relative active space size computed at the LNO-CCSD(T)/PM/cc-pVDZ[90] level for the Baker test set. The reference data were obtained at the CCSD(T)/cc-pVDZ level, and the geometries were optimized at the MP2/cc-pVDZ level.

IV.2 Geometry Optimization

Geometry optimizations were also performed for molecules in the Baker test set, through an interface with the GeomeTRIC library.[91] The default optimization convergence criteria were employed, including a nuclear gradient root mean square less than $3\times 10^{-4}$ hartree/Bohr and an absolute energy change less than $1\times 10^{-6}$ hartree.

We compare the LNO-CCSD(T) results to the canonical CCSD(T) reference, and present the errors in Table 1. The observed errors are small for both bond lengths and various types of angles, even with the loose LNO cutoff of $\zeta=10^{-4}$ . It is worth mentioning that the geometry optimization for benzene with $\zeta=10^{-4}$ did not converge due to large gradient fluctuations caused by symmetry breaking. Nevertheless, such an issue can be avoided by using a smaller $\zeta$ or by employing the IAO local orbitals.

The present implementation of LNO-CCSD(T) has a formal scaling of $O(N^{5}+\tilde{N}^{7})$ , where $N$ is the number of basis functions, and $\tilde{N}$ represents the dimensions of the local correlation domains. The $O(N^{5})$ part of the computational complexity originates from the global MP2 correction and the LNO construction scheme, while the $O(\tilde{N}^{7})$ part is due to the local CCSD(T) calculations. For the chemical systems for which one would use CCSD(T), $\tilde{N}$ does not need to grow with the system size for a given energy accuracy per atom, making LNO-CCSD(T) orders of magnitude more efficient than canonical CCSD(T), especially for large $N$ (see Fig. LABEL:fig:time_pm_TZ and the corresponding discussion). However, employing a looser LNO cutoff (and thus a smaller $\tilde{N}$ ) might introduce more severe discontinuities into the potential energy surface, leading to more iterations to converge the geometry optimization. Overall, our tests suggest that our LNO-CCSD(T) method is capable of accurate geometry optimizations with computational scaling comparable to that of canonical MP2, albeit potentially with a large prefactor.

Table 1: Mean absolute error (MAE) and maximum error (max) in bond lengths, angles, dihedral angles and out-of-plane angles of the geometries optimized at the LNO-CCSD(T)/PM/cc-pVDZ level for the Baker test set. The reference geometries were optimized at the canonical CCSD(T) level with the same basis set.

$\zeta$		bond length (Å)	angle (°)	dihedral angle (°)	out-of-plane angle (°)
$10^{-4}$ ^a	MAE	0.0008	0.038	0.059	0.008
	max	0.0082	0.363	0.803	0.152
$10^{-5}$	MAE	0.0002	0.014	0.025	0.002
	max	0.0019	0.085	0.774	0.028

^a The geometry used for benzene was obtained from iteration 100 of the geometry optimization.

IV.3 Relaxed Density Matrix

Most first-order static response properties can be computed with the zeroth-order density matrix, which is readily available through AD as the energy derivative with respect to the unperturbed Hamiltonian. This is especially useful for the LNO-CC method, which never calculates a global wavefunction that could be used to calculate the density matrix. Specifically, the one-electron reduced density matrix may be expressed as

D_{\mu\nu}=\frac{\partial E_{\text{LNO-CC}}}{\partial h_{\mu\nu}^{\text{core}}}\;,

(7)

where $\mathbf{h}^{\text{core}}$ is the one-electron core Hamiltonian matrix. Note that orbital relaxation is incorporated here by using the Hamiltonian in the atomic orbital basis. In the following, we show an example of computing bond orders and Mössbauer spectroscopy parameters for a [NiFe]-hydrogenase active site model system using the LNO-CCSD(T) method.

The crystal structures[92] of the NAD⁺-reducing soluble [NiFe]-hydrogenase (SH) from Hydrogenophilus thermoluteolus (Ht) reveal an unusual arrangement at the oxidized active site [see Fig. 2 (a)]. There, the nickel center adopts a distorted octahedral six-coordinate configuration, featuring three bridging cysteines, one terminal cysteine, and a bidentate glutamate coordination. Meanwhile, the iron site is coordinated by two cyanide and one carbon monoxide ligand. Notably, the IR spectrum[93, 94] of the oxidized state exhibits a unique CO vibration band at 1993 cm^-1, which is distinct from all other [NiFe]-hydrogenases. Such spectral features were proposed to originate from a biologically unprecedented Ni(IV)/Fe(II) ground state, supported by density functional theory (DFT) calculations.[94] However, a subsequent study, also conducted at the DFT level, suggests that the spectral properties and the coordination geometry of the oxidized Ht-SH may be attributed to an open-shell singlet Ni(III)/Fe(III) state instead.[95]

It would be interesting to study the structural and spectral properties of the Ht-SH active site beyond mean-field theory. Here, we employ the LNO-CCSD(T) method to calculate bond orders and Mössbauer spectroscopy parameters, offering a comparative analysis with DFT results. Note that the primary aim of the current work is to demonstrate the feasibility of our LNO-CC method, while a comprehensive investigation of the electronic structure of the Ht-SH system will be deferred to future studies.

The active site of Ht-SH in the fully oxidized state is represented by the “model 20” structure [see Fig. 2 (b)] from Ref. 94. Its geometry, optimized at the DFT/TPSSh[96]/def2-TZVP[97] level, is taken from Ref. 95. Only the closed-shell singlet ground state is considered below.

In Table 2, the orbital resolved bond orders (defined in Eq. LABEL:eq:bond_order) are presented for the two metal centers. Minor contributions to bonding are omitted, and only representative bonds among those with similar chemical environments are shown. (See Tables LABEL:tab:NiFe_bond_order1_full and LABEL:tab:NiFe_bond_order for the complete data.) We compare the results computed by DFT with those obtained by the LNO-CCSD(T) method. For the latter, IAOs were employed as the local orbitals, and an LNO cutoff of $\zeta=2\times 10^{-5}$ was adopted, leading to fragments with correlation domains comprising approximately $80$ occupied and $200$ virtual orbitals at most. (Such choices are discussed in the supporting information Sec. LABEL:sec:HtSH.) The key findings are the following: (i) Bond orders calculated by the two methods qualitatively agree with each other. (ii) The $4s$ orbitals of both Fe and Ni do not contribute significantly to bonding. (iii) Relatively strong covalent bonding is observed between Fe and CO. (iv) The three bridging cysteines exhibit weak bonding interactions with the two metal centers, characterized by the stronger Ni–S bonds compared to the Fe–S bonds. Interestingly, the LNO-CCSD(T) method tends to predict non-bonding characters for the Fe–S bonds.

Table 2: Bond orders for the Ht-SH active site model system, computed at the DFT/TPSSh and LNO-CCSD(T)/IAO levels with the def2-TZVP basis set.

Bond	TPSSh	LNO-CCSD(T)
Fe(3d)–C¹(2p)N	0.24	0.17
Fe(3d)–C(2p)O	0.57	0.87
Fe(3d)–S¹(3p)Cys	0.21	0.09
Fe(3d)–S²(3p)Cys	0.17	0.07
Ni(3d)–S¹(3p)Cys	0.26	0.31
Ni(3d)–S²(3p)Cys	0.34	0.14
Ni(3d)–S⁴(3p)Cys	0.33	0.32
Ni(4s)–O¹(2p)	0.12	0.09

We have also computed the quadrupole splittings $\Delta_{v}$ (defined in Eq. LABEL:eq:Delta_v) for both metal centers, as shown in Table 3. It can be seen that DFT and LNO-CCSD(T) give consistent results for Fe, whereas noticeable discrepancies emerge for Ni. This suggests that the electron density distribution around the Ni atom is described differently by the two methods, warranting further investigation.

Table 3: Quadrupole splittings for the iron and nickel centers of the Ht-SH active site model system, computed at the DFT/TPSSh and LNO-CCSD(T)/IAO levels with the def2-TZVP basis set.

Method	$\Delta_{v}$ (mm/s)
	⁵⁷Fe	⁶¹Ni
TPSSh	0.46	0.17
LNO-CCSD(T)	0.43	0.35

IV.4 Ab Initio Molecular Dynamics

As the final example, we present the calculation of the IR spectrum for a protonated water hexamer, \chH^+(H2O)6, via ab initio molecular dynamics (AIMD). The IR intensity ( $A$ ) is proportional to the dipole-dipole correlation function. In the frequency domain, it can be expressed as

A(\omega)\propto\int\mathrm{d}t\langle\bm{\dot{\mu}}(0)\cdot\bm{\dot{\mu}}(t)\rangle e^{-i\omega t},

(8)

where the dipole moment, $\bm{\mu}$ , was computed through AD as the energy derivative with respect to the external electric field, and its time derivative was computed via finite difference. Additionally, the bracket in Eq. 8 indicates an ensemble average over the MD trajectories, which were computed at the LNO-CCSD(T) level using the cc-pVTZ[90] basis set. The initial geometries were prepared by Avogadro[98] to represent likely conformers of the protonated water cluster, and then loosely relaxed at the DFT/ $\omega$ B97X[99]/cc-pVDZ level. This results in four distinct structures from which we further seed the dynamics: one Zundel-like (\chH5O2^+) cation, one Eigen-like [\chH3O^+(H2O)3] cation, and two with a similar four-water-ring [\chH^+(H2O)4] structure (see Fig. LABEL:fig:water). The LNO-CCSD(T) calculations employed the PM local orbitals and an LNO cutoff of $\zeta=5\times 10^{-6}$ . Such a setup was found to reproduce the canonical CCSD(T)/cc-pVTZ energies of the DFT-relaxed geometries to well within chemical accuracy (with a mean absolute deviation of 0.04 kcal/mol and with a largest deviation of 0.06 kcal/mol). The MD equilibration in the canonical ( $NVT$ ) ensemble was then performed starting from the four DFT-relaxed geometries. The temperature was controlled by a Langevin thermostat at 50 K and the nuclear dynamics was integrated using a time step of 1 fs. The production runs using the microcanonical ( $NVE$ ) ensemble were performed for 10 ps following 1 ps of $NVT$ equilibration. For the $NVE$ simulations, we employed a longer time step and adopted the multiple time-stepping (MTS) integrator to reduce the time step error of the AIMD. For this we utilized a 2 fs outer MD time step with LNO-CCSD(T) forces, with a 0.5 fs inner time step with machine-learned forces, following the machine-learning (ML) accelerated AIMD approach[100]. We chose the Allegro framework[101] to train the ML reference potential on the $NVT$ -sampled configurations (1001 configurations for each conformer) and their energies and forces. To benchmark the error of the MTS integrator, a separate $NVE$ run (without applying the MTS method) using a time step of 0.5 fs was compared with the MTS dynamics. Agreement between the two resulting IR spectra validates our MTS approach (see Fig. LABEL:fig:mts_compare). As discussed earlier, the potential energy surface (PES) calculated by the LNO-CC method is not strictly continuous, and thus an energy drift is expected in the $NVE$ MD. However, the observed energy drift in our simulations is negligible ( $<0.05$ kcal/mol/ps), as the result of using a small enough $\zeta$ .

By computing the IR spectrum using Eq. 8 with MD simulations on the LNO-CCSD(T) PES, we have explicitly accounted for the electron correlation at the CCSD(T) level (within the local approximation) and included anharmonic vibrational effects. Assuming the electronic structure and sampling to be converged, the only remaining physical component to be included is the nuclear quantum effect (NQE). While NQEs can in principle be incorporated through path integral MD[103, 104] or by explicitly solving the vibrational Schrödinger equation, [105, 106, 107] herein we restrict ourselves to classical nuclear dynamics, and thus quantitative agreement with the experiment is not expected. As shown in Figure 3, the computed IR spectra of the Zundel- and Eigen-like conformers agree qualitatively with the experimental IR²MS² spectra, while blue-shifted signals are observed compared to the experimental peak positions. We attribute the blue shifting to the missing NQEs instead of the basis set incompleteness error, since we found that enlarging the basis set from cc-pVTZ to aug-cc-pVQZ red shifts the OH stretching bands by only $\sim$ 20 cm^-1 in a gas-phase harmonic analysis for a single water molecule at the canonical CCSD(T) level. The experimental IR²MS² signals[102] were obtained by probing the bands at 3159 cm^-1 and 3715 cm^-1 (denoted as $a_{6}$ and $a_{3}$ , respectively), which are assigned to the OH stretches in the Zundel and Eigen cations, respectively. The absence of $a_{6}$ / $a_{3}$ signals in the calculated IR spectra of the Eigen/Zundel-like conformers clearly suggests that the two frequencies are unique features of the Zundel/Eigen cations. Moreover, the IR spectra computed for the four-water-ring-like conformers show distinct signatures (see Fig. LABEL:fig:ir_all), unambiguously validating that the IR²MS² measurements correspond to the signals of interest from the Zundel- and Eigen-like conformers.

V Conclusions

In this work, we introduced a differentiable implementation of a local coupled cluster theory, utilizing the newly improved AD framework provided by the PySCFAD package. Calculations of first-order response properties, including nuclear gradients, dipole moments, electric field gradients and relaxed density matrices, demonstrate the feasibility of our framework for analytic derivative computations involving complex computational workflows. Moreover, the application of our method to geometry optimizations and AIMD simulations for medium-sized molecular systems further indicates that PySCFAD is not only a useful platform to rapidly prototype new methodologies but also effective for production-level calculations.

While significant progress has been made, there remain several challenges to address. Firstly, computing higher-order response properties with the current LNO-CC method can become prohibitively expensive. Our implementation has a formal scaling of $O(N^{5})$ , where $N$ is the number of basis functions. Moving to each higher order of response increases both computational and memory complexities by a factor equal to the dimension of the perturbation variable. This may be manageable for low-dimensional perturbations such as the electric field, which has a dimension of three. However, for variables like nuclear coordinates, the resultant cost elevation could be substantial. One potential solution involves employing the various domain truncation algorithms [22, 30] to mitigate the overall computational expense. These algorithms can be seamlessly integrated into our differentiable LNO-CC implementation, assuming a well defined energy expression exists.

Secondly, only static response properties are readily computable at the moment, whereas a straightforward formalism for calculating dynamic response properties via analytic derivatives is still lacking. As discussed before,[76] dynamic response properties may be computed as the derivatives of the quasienergy, which is defined as the time-averaged expectation value of $\hat{H}-i\partial/\partial t$ over the time-dependent wavefunction.[108] However, a time-dependent implementation of local correlation methods remains to be developed.

Acknowledgements.

This work was primarily supported by the United States Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, FWP LANLE3F2 awarded to Los Alamos National Laboratory under Triad National Security, LLC (‘Triad’) contract grant no. 89233218CNA000001, subaward C2448 to the California Institute of Technology. Additional support for GKC was provided by the Camille and Henry Dreyfus Foundation via a grant from the program ”Machine Learning in the Chemical Sciences and Engineering” Some of the calculations were performed at the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy, Office of Science User Facility located at the Lawrence Berkeley National Laboratory.

Author Contributions

X. Z. and C. L. contributed equally to this work. X. Z., C. L., and G. K.-L. C. designed this project. H.-Z. Y. and T. C. B. contributed to the original implementation of the LNO-CC method, based on which X. Z. developed the differentiable version used here. H.-Z. Y. and T. C. B. also provided help and advice with the methodology development in the early part of the project. C. L. developed the AIMD simulation workflow. X. Z. and C. L. performed the calculations and wrote the original draft. X. Z., C. L., and G. K.-L. C. validated the data, and all contributed to the reviewing and editing of the manuscript.

Data Availability

The data that supports the findings of this study are available within the article, and/or from the corresponding author upon reasonable request. The PySCFAD source code can be found at https://github.com/fishjojo/pyscfad.

References

Pulay [1983] P. Pulay, “Localizability of dynamic electron correlation,” Chemical physics letters 100, 151–154 (1983).
Sæbø and Pulay [1985] S. Sæbø and P. Pulay, “Local configuration interaction: An efficient approach for larger molecules,” Chemical physics letters 113, 13–18 (1985).
Pulay and Saebø [1986] P. Pulay and S. Saebø, “Orbital-invariant formulation and second-order gradient evaluation in Møller-Plesset perturbation theory,” Theoretica chimica acta 69, 357–368 (1986).
Sæbø and Pulay [1993] S. Sæbø and P. Pulay, “Local treatment of electron correlation,” Annual Review of Physical Chemistry 44, 213–236 (1993).
Boughton and Pulay [1993] J. W. Boughton and P. Pulay, “Comparison of the boys and Pipek–Mezey localizations in the local correlation approach and automatic virtual basis selection,” Journal of computational chemistry 14, 736–740 (1993).
Hampel and Werner [1996] C. Hampel and H.-J. Werner, “Local treatment of electron correlation in coupled cluster theory,” The Journal of chemical physics 104, 6286–6297 (1996).
Schütz and Werner [2000] M. Schütz and H.-J. Werner, “Local perturbative triples correction (T) with linear cost scaling,” Chemical Physics Letters 318, 370–378 (2000).
Schütz and Werner [2001] M. Schütz and H.-J. Werner, “Low-order scaling local electron correlation methods. IV. linear scaling local coupled-cluster (LCCSD),” The Journal of Chemical Physics 114, 661–681 (2001).
Schütz [2002a] M. Schütz, “Low-order scaling local electron correlation methods. V. connected triples beyond (T): Linear scaling local CCSDT-1b,” The Journal of chemical physics 116, 8772–8785 (2002a).
Schütz [2002b] M. Schütz, “A new, fast, semi-direct implementation of linear scaling local coupled cluster theory,” Physical Chemistry Chemical Physics 4, 3941–3947 (2002b).
Li, Ma, and Jiang [2002] S. Li, J. Ma, and Y. Jiang, “Linear scaling local correlation approach for solving the coupled cluster equations of large systems,” Journal of computational chemistry 23, 237–244 (2002).
Li et al. [2006] S. Li, J. Shen, W. Li, and Y. Jiang, “An efficient implementation of the “cluster-in-molecule” approach for local electron correlation calculations,” The Journal of chemical physics 125, 074109 (2006).
Li et al. [2009] W. Li, P. Piecuch, J. R. Gour, and S. Li, “Local correlation calculations using standard and renormalized coupled-cluster approaches,” The Journal of chemical physics 131, 114109 (2009).
Neese, Wennmohs, and Hansen [2009] F. Neese, F. Wennmohs, and A. Hansen, “Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method,” The Journal of chemical physics 130, 114108 (2009).
Neese, Hansen, and Liakos [2009] F. Neese, A. Hansen, and D. G. Liakos, “Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis,” The Journal of chemical physics 131, 064103 (2009).
Werner and Schütz [2011] H.-J. Werner and M. Schütz, “An efficient local coupled cluster method for accurate thermochemistry of large systems,” The Journal of Chemical Physics 135, 144116 (2011).
Rolik and Kállay [2011] Z. Rolik and M. Kállay, “A general-order local coupled-cluster method based on the cluster-in-molecule approach,” The Journal of chemical physics 135, 104111 (2011).
Yang et al. [2011] J. Yang, Y. Kurashige, F. R. Manby, and G. K. L. Chan, “Tensor factorizations of local second-order Møller–Plesset theory,” The Journal of Chemical Physics 134, 044123 (2011).
Kurashige et al. [2012] Y. Kurashige, J. Yang, G. K.-L. Chan, and F. R. Manby, “Optimization of orbital-specific virtuals in local Møller-Plesset perturbation theory,” The Journal of Chemical Physics 136, 124106 (2012).
Yang et al. [2012] J. Yang, G. K.-L. Chan, F. R. Manby, M. Schütz, and H.-J. Werner, “The orbital-specific-virtual local coupled cluster singles and doubles method,” The Journal of Chemical Physics 136, 144105 (2012).
Schütz et al. [2013] M. Schütz, J. Yang, G. K. Chan, F. R. Manby, and H.-J. Werner, “The orbital-specific virtual local triples correction: OSV-L (T),” The Journal of Chemical Physics 138, 054109 (2013).
Rolik et al. [2013] Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. Kállay, “An efficient linear-scaling CCSD (T) method based on local natural orbitals,” The Journal of chemical physics 139, 094105 (2013).
Riplinger and Neese [2013] C. Riplinger and F. Neese, “An efficient and near linear scaling pair natural orbital based local coupled cluster method,” The Journal of chemical physics 138, 034106 (2013).
Riplinger et al. [2013] C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, “Natural triple excitations in local coupled cluster calculations with pair natural orbitals,” The Journal of chemical physics 139, 134101 (2013).
Sparta and Neese [2014] M. Sparta and F. Neese, “Chemical applications carried out by local pair natural orbital based coupled-cluster methods,” Chemical Society Reviews 43, 5032–5041 (2014).
Schmitz, Hättig, and Tew [2014] G. Schmitz, C. Hättig, and D. P. Tew, “Explicitly correlated PNO-MP2 and PNO-CCSD and their application to the S66 set and large molecular systems,” Physical Chemistry Chemical Physics 16, 22167–22178 (2014).
Kállay [2015] M. Kállay, “Linear-scaling implementation of the direct random-phase approximation,” The Journal of chemical physics 142, 204105 (2015).
Liakos et al. [2015] D. G. Liakos, M. Sparta, M. K. Kesharwani, J. M. Martin, and F. Neese, “Exploring the accuracy limits of local pair natural orbital coupled-cluster theory,” Journal of chemical theory and computation 11, 1525–1539 (2015).
Liakos and Neese [2015] D. G. Liakos and F. Neese, “Is it possible to obtain coupled cluster quality energies at near density functional theory cost? domain-based local pair natural orbital coupled cluster vs modern density functional theory,” Journal of chemical theory and computation 11, 4054–4063 (2015).
Riplinger et al. [2016] C. Riplinger, P. Pinski, U. Becker, E. F. Valeev, and F. Neese, “Sparse maps—a systematic infrastructure for reduced-scaling electronic structure methods. II. linear scaling domain based pair natural orbital coupled cluster theory,” The Journal of chemical physics 144, 024109 (2016).
Schmitz and Hättig [2016] G. Schmitz and C. Hättig, “Perturbative triples correction for local pair natural orbital based explicitly correlated CCSD (F12*) using laplace transformation techniques,” The Journal of Chemical Physics 145, 234107 (2016).
Guo et al. [2016] Y. Guo, K. Sivalingam, E. F. Valeev, and F. Neese, “Sparsemaps—a systematic infrastructure for reduced-scaling electronic structure methods. III. linear-scaling multireference domain-based pair natural orbital N-electron valence perturbation theory,” The Journal of chemical physics 144, 094111 (2016).
Pavošević et al. [2017] F. Pavošević, C. Peng, P. Pinski, C. Riplinger, F. Neese, and E. F. Valeev, “Sparsemaps—a systematic infrastructure for reduced scaling electronic structure methods. V. linear scaling explicitly correlated coupled-cluster method with pair natural orbitals,” The Journal of chemical physics 146, 174108 (2017).
Ma et al. [2017] Q. Ma, M. Schwilk, C. Köppl, and H.-J. Werner, “Scalable electron correlation methods. 4. parallel explicitly correlated local coupled cluster with pair natural orbitals (PNO-LCCSD-F12),” Journal of chemical theory and computation 13, 4871–4896 (2017).
Ma and Werner [2018] Q. Ma and H.-J. Werner, “Scalable electron correlation methods. 5. parallel perturbative triples correction for explicitly correlated local coupled cluster with pair natural orbitals,” Journal of Chemical Theory and Computation 14, 198–215 (2018).
Guo, Becker, and Neese [2018] Y. Guo, U. Becker, and F. Neese, “Comparison and combination of “direct” and fragment based local correlation methods: Cluster in molecules and domain based local pair natural orbital perturbation and coupled cluster theories,” The Journal of Chemical Physics 148, 124117 (2018).
Liakos, Guo, and Neese [2019] D. G. Liakos, Y. Guo, and F. Neese, “Comprehensive benchmark results for the domain based local pair natural orbital coupled cluster method (DLPNO-CCSD (T)) for closed-and open-shell systems,” The Journal of Physical Chemistry A 124, 90–100 (2019).
Nagy, Samu, and Kállay [2018] P. R. Nagy, G. Samu, and M. Kállay, “Optimization of the linear-scaling local natural orbital CCSD (T) method: Improved algorithm and benchmark applications,” Journal of Chemical Theory and Computation 14, 4193–4215 (2018).
Nagy and Kállay [2019] P. R. Nagy and M. Kállay, “Approaching the basis set limit of CCSD (T) energies for large molecules with local natural orbital coupled-cluster methods,” Journal of Chemical Theory and Computation 15, 5275–5298 (2019).
Ni, Li, and Li [2019] Z. Ni, W. Li, and S. Li, “Fully optimized implementation of the cluster-in-molecule local correlation approach for electron correlation calculations of large systems,” Journal of Computational Chemistry 40, 1130–1140 (2019).
Wang et al. [2019] Y. Wang, Z. Ni, W. Li, and S. Li, “Cluster-in-molecule local correlation approach for periodic systems,” Journal of Chemical Theory and Computation 15, 2933–2943 (2019).
Ni et al. [2021] Z. Ni, Y. Guo, F. Neese, W. Li, and S. Li, “Cluster-in-molecule local correlation method with an accurate distant pair correction for large systems,” Journal of Chemical Theory and Computation 17, 756–766 (2021).
Wang et al. [2022] Y. Wang, Z. Ni, F. Neese, W. Li, Y. Guo, and S. Li, “Cluster-in-molecule method combined with the domain-based local pair natural orbital approach for electron correlation calculations of periodic systems,” Journal of Chemical Theory and Computation 18, 6510–6521 (2022).
Li et al. [2023] W. Li, Y. Wang, Z. Ni, and S. Li, “Cluster-in-molecule local correlation method for dispersion interactions in large systems and periodic systems,” Accounts of Chemical Research 56, 3462–3474 (2023).
Ye and Berkelbach [2024a] H.-Z. Ye and T. C. Berkelbach, “Ab initio surface chemistry with chemical accuracy,” (2024a), arXiv:2309.14640 [cond-mat.mtrl-sci] .
Raghavachari et al. [1989] K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, “A fifth-order perturbation comparison of electron correlation theories,” Chemical Physics Letters 157, 479–483 (1989).
Ye and Berkelbach [2024b] H.-Z. Ye and T. C. Berkelbach, “Adsorption and vibrational spectroscopy of CO on the surface of MgO from periodic local coupled-cluster theory,” (2024b), arXiv:2309.14651 [cond-mat.mtrl-sci] .
El Azhary et al. [1998] A. El Azhary, G. Rauhut, P. Pulay, and H.-J. Werner, “Analytical energy gradients for local second-order Møller–Plesset perturbation theory,” The Journal of chemical physics 108, 5185–5193 (1998).
Rauhut and Werner [2001] G. Rauhut and H.-J. Werner, “Analytical energy gradients for local coupled-cluster methods,” Physical Chemistry Chemical Physics 3, 4853–4862 (2001).
Schütz et al. [2004] M. Schütz, H.-J. Werner, R. Lindh, and F. R. Manby, “Analytical energy gradients for local second-order Møller–Plesset perturbation theory using density fitting approximations,” The Journal of chemical physics 121, 737–750 (2004).
Dornbach and Werner [2019] M. Dornbach and H.-J. Werner, “Analytical energy gradients for local second-order Møller-Plesset perturbation theory using intrinsic bond orbitals,” Molecular Physics 117, 1252–1263 (2019).
Ledermüller, Kats, and Schütz [2013] K. Ledermüller, D. Kats, and M. Schütz, “Local CC2 response method based on the laplace transform: Orbital-relaxed first-order properties for excited states,” The Journal of Chemical Physics 139, 084111 (2013).
Ledermüller and Schütz [2014] K. Ledermüller and M. Schütz, “Local CC2 response method based on the laplace transform: Analytic energy gradients for ground and excited states,” The Journal of Chemical Physics 140, 164113 (2014).
Gauss and Werner [2000] J. Gauss and H.-J. Werner, “NMR chemical shift calculations within local correlation methods: the GIAO-LMP2 approach,” Physical Chemistry Chemical Physics 2, 2083–2090 (2000).
Loibl and Schütz [2012] S. Loibl and M. Schütz, “NMR shielding tensors for density fitted local second-order Møller-Plesset perturbation theory using gauge including atomic orbitals,” The Journal of Chemical Physics 137, 084107 (2012).
Loibl and Schütz [2014] S. Loibl and M. Schütz, “Magnetizability and rotational g tensors for density fitted local second-order Møller-Plesset perturbation theory using gauge-including atomic orbitals,” The Journal of Chemical Physics 141, 024108 (2014).
Pinski and Neese [2018] P. Pinski and F. Neese, “Communication: Exact analytical derivatives for the domain-based local pair natural orbital mp2 method (DLPNO-MP2),” The Journal of Chemical Physics 148, 031101 (2018).
Pinski and Neese [2019] P. Pinski and F. Neese, “Analytical gradient for the domain-based local pair natural orbital second order Møller-Plesset perturbation theory method (DLPNO-MP2),” The Journal of Chemical Physics 150, 164102 (2019).
Stoychev et al. [2021] G. L. Stoychev, A. A. Auer, J. Gauss, and F. Neese, “DLPNO-MP2 second derivatives for the computation of polarizabilities and NMR shieldings,” The Journal of Chemical Physics 154, 164110 (2021).
Datta, Kossmann, and Neese [2016] D. Datta, S. Kossmann, and F. Neese, “Analytic energy derivatives for the calculation of the first-order molecular properties using the domain-based local pair-natural orbital coupled-cluster theory,” The Journal of Chemical Physics 145, 114101 (2016).
Zhou, Liang, and Yang [2019] R. Zhou, Q. Liang, and J. Yang, “Complete OSV-MP2 analytical gradient theory for molecular structure and dynamics simulations,” Journal of Chemical Theory and Computation 16, 196–210 (2019).
Russ and Crawford [2004] N. J. Russ and T. D. Crawford, “Local correlation in coupled cluster calculations of molecular response properties,” Chemical physics letters 400, 104–111 (2004).
Russ and Crawford [2008] N. J. Russ and T. D. Crawford, “Local correlation domains for coupled cluster theory: optical rotation and magnetic-field perturbations,” Physical Chemistry Chemical Physics 10, 3345–3352 (2008).
McAlexander, Mach, and Crawford [2012] H. R. McAlexander, T. J. Mach, and T. D. Crawford, “Localized optimized orbitals, coupled cluster theory, and chiroptical response properties,” Physical Chemistry Chemical Physics 14, 7830–7836 (2012).
McAlexander and Crawford [2016] H. R. McAlexander and T. D. Crawford, “A comparison of three approaches to the reduced-scaling coupled cluster treatment of non-resonant molecular response properties,” Journal of Chemical Theory and Computation 12, 209–222 (2016).
D’Cunha and Crawford [2021] R. D’Cunha and T. D. Crawford, “PNO++: Perturbed pair natural orbitals for coupled cluster linear response theory,” Journal of Chemical Theory and Computation 17, 290–301 (2021).
Maclaurin, Duvenaud, and Adams [2015] D. Maclaurin, D. Duvenaud, and R. P. Adams, “Autograd: Effortless gradients in numpy,” in ICML 2015 AutoML workshop, Vol. 238 (2015).
Bradbury et al. [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, “JAX: composable transformations of Python+NumPy programs,” (2018), available at http://github.com/google/jax.
Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al., “Pytorch: An imperative style, high-performance deep learning library,” Advances in neural information processing systems 32 (2019).
Frostig, Johnson, and Leary [2018] R. Frostig, M. J. Johnson, and C. Leary, “Compiling machine learning programs via high-level tracing,” Systems for Machine Learning 4 (2018).
Reed et al. [2022] J. Reed, Z. DeVito, H. He, A. Ussery, and J. Ansel, “torch.fx: Practical program capture and transformation for deep learning in python,” in Proceedings of Machine Learning and Systems, Vol. 4, edited by D. Marculescu, Y. Chi, and C. Wu (2022) pp. 638–651.
Tamayo-Mendoza et al. [2018] T. Tamayo-Mendoza, C. Kreisbeck, R. Lindh, and A. Aspuru-Guzik, “Automatic differentiation in quantum chemistry with applications to fully variational Hartree–Fock,” ACS central science 4, 559–566 (2018).
Song, Martínez, and Neaton [2021] C. Song, T. J. Martínez, and J. B. Neaton, “A diagrammatic approach for automatically deriving analytical gradients of tensor hyper-contracted electronic structure methods,” The Journal of Chemical Physics 155, 024108 (2021).
Abbott et al. [2021] A. S. Abbott, B. Z. Abbott, J. M. Turney, and H. F. Schaefer III, “Arbitrary-order derivatives of quantum chemical methods via automatic differentiation,” The journal of physical chemistry letters 12, 3232–3239 (2021).
Kasim, Lehtola, and Vinko [2022] M. F. Kasim, S. Lehtola, and S. M. Vinko, “DQC: A python program package for differentiable quantum chemistry,” The Journal of chemical physics 156 (2022).
Zhang and Chan [2022] X. Zhang and G. K. Chan, “Differentiable quantum chemistry with PySCF for molecules and materials at the mean-field level and beyond,” The Journal of Chemical Physics 157, 204801 (2022).
Vargas-Hernández et al. [2023] R. A. Vargas-Hernández, K. Jorner, R. Pollice, and A. Aspuru-Guzik, “Inverse molecular design and parameter optimization with Hückel theory using automatic differentiation,” The Journal of Chemical Physics 158, 104801 (2023).
Mahajan et al. [2023] A. Mahajan, J. S. Kurian, J. Lee, D. R. Reichman, and S. Sharma, “Response properties in phaseless auxiliary field quantum Monte Carlo,” The Journal of Chemical Physics 159, 184101 (2023).
M Casares et al. [2024] P. A. M Casares, J. S. Baker, M. Medvidović, R. d. Reis, and J. M. Arrazola, “GradDFT. a software library for machine learning enhanced density functional theory,” The Journal of Chemical Physics 160, 062501 (2024).
Kitaura et al. [1999] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayasi, “Fragment molecular orbital method: an approximate computational method for large molecules,” Chemical Physics Letters 313, 701–706 (1999).
Sun et al. [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, “Recent developments in the PySCF program package,” The Journal of Chemical Physics 153, 024109 (2020).
Dalcín, Paz, and Storti [2005] L. Dalcín, R. Paz, and M. Storti, “MPI for Python,” Journal of Parallel and Distributed Computing 65, 1108–1115 (2005).
Scheurer and Schwarz [2000a] P. Scheurer and W. Schwarz, “Externally localized molecular orbitals: A numerical investigation of localization degeneracy,” International Journal of Quantum Chemistry 76, 420–427 (2000a).
Scheurer and Schwarz [2000b] P. Scheurer and W. Schwarz, “Continuous degeneracy of sets of localized orbitals,” International Journal of Quantum Chemistry 76, 428–433 (2000b).
Whitten [1973] J. L. Whitten, “Coulombic potential energy integrals and approximations,” The Journal of Chemical Physics 58, 4496–4501 (1973).
Dunlap, Connolly, and Sabin [1979] B. I. Dunlap, J. Connolly, and J. Sabin, “On some approximations in applications of X $\alpha$ theory,” The Journal of Chemical Physics 71, 3396–3402 (1979).
Pipek and Mezey [1989] J. Pipek and P. G. Mezey, “A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions,” The Journal of Chemical Physics 90, 4916–4926 (1989).
Baker [1993] J. Baker, “Techniques for geometry optimization: A comparison of cartesian and natural internal coordinates,” Journal of Computational Chemistry 14, 1085–1100 (1993).
Knizia [2013] G. Knizia, “Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts,” Journal of chemical theory and computation 9, 4834–4843 (2013).
Dunning Jr [1989] T. H. Dunning Jr, “Gaussian basis sets for use in correlated molecular calculations. I. the atoms boron through neon and hydrogen,” The Journal of chemical physics 90, 1007–1023 (1989).
Wang and Song [2016] L.-P. Wang and C. Song, “Geometry optimization made simple with translation and rotation coordinates,” The Journal of chemical physics 144 (2016).
Shomura et al. [2017] Y. Shomura, M. Taketa, H. Nakashima, H. Tai, H. Nakagawa, Y. Ikeda, M. Ishii, Y. Igarashi, H. Nishihara, K.-S. Yoon, S. Ogo, S. Hirota, and Y. Higuchi, “Structural basis of the redox switches in the NAD⁺-reducing soluble [NiFe]-hydrogenase,” Science 357, 928–932 (2017).
Preissler et al. [2018] J. Preissler, S. Wahlefeld, C. Lorent, C. Teutloff, M. Horch, L. Lauterbach, S. P. Cramer, I. Zebger, and O. Lenz, “Enzymatic and spectroscopic properties of a thermostable [NiFe]‑hydrogenase performing H₂-driven NAD⁺-reduction in the presence of O₂,” Biochimica et Biophysica Acta (BBA) - Bioenergetics 1859, 8–18 (2018).
Kulka-Peschke et al. [2022] C. J. Kulka-Peschke, A.-C. Schulz, C. Lorent, Y. Rippers, S. Wahlefeld, J. Preissler, C. Schulz, C. Wiemann, C. C. M. Bernitzky, C. Karafoulidi-Retsou, S. L. D. Wrathall, B. Procacci, H. Matsuura, G. M. Greetham, C. Teutloff, L. Lauterbach, Y. Higuchi, M. Ishii, N. T. Hunt, O. Lenz, I. Zebger, and M. Horch, “Reversible glutamate coordination to high-valent nickel protects the active site of a [NiFe] hydrogenase from oxygen,” Journal of the American Chemical Society 144, 17022–17032 (2022).
Kumar and Stein [2023] R. Kumar and M. Stein, “The fully oxidized state of the glutamate coordinated O₂-tolerant [nife]-hydrogenase shows a Ni(III)/Fe(III) open-shell singlet ground state,” Journal of the American Chemical Society 145, 10954–10959 (2023).
Staroverov et al. [2003] V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. Perdew, “Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes,” The Journal of chemical physics 119, 12129–12137 (2003).
Weigend and Ahlrichs [2005] F. Weigend and R. Ahlrichs, “Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy,” Physical Chemistry Chemical Physics 7, 3297–3305 (2005).
Hanwell et al. [2012] M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek, and G. R. Hutchison, “Avogadro: an advanced semantic chemical editor, visualization, and analysis platform,” Journal of cheminformatics 4, 1–17 (2012).
Chai and Head-Gordon [2008] J.-D. Chai and M. Head-Gordon, “Systematic optimization of long-range corrected hybrid density functionals,” The Journal of chemical physics 128, 084106 (2008).
Li and Voth [2022] C. Li and G. A. Voth, “Using machine learning to greatly accelerate path integral ab initio molecular dynamics,” Journal of Chemical Theory and Computation 18, 599–604 (2022).
Musaelian et al. [2023] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, “Learning local equivariant representations for large-scale atomistic dynamics,” Nature Communications 14, 579 (2023).
Heine et al. [2013] N. Heine, M. R. Fagiani, M. Rossi, T. Wende, G. Berden, V. Blum, and K. R. Asmis, “Isomer-selective detection of hydrogen-bond vibrations in the protonated water hexamer,” Journal of the American Chemical Society 135, 8266–8273 (2013).
Rossi, Ceriotti, and Manolopoulos [2014] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, “How to remove the spurious resonances from ring polymer molecular dynamics,” The Journal of chemical physics 140, 234116 (2014).
Paesani and Voth [2010] F. Paesani and G. A. Voth, “A quantitative assessment of the accuracy of centroid molecular dynamics for the calculation of the infrared spectrum of liquid water,” The Journal of chemical physics 132, 014105 (2010).
Kjaergaard et al. [2008] H. G. Kjaergaard, A. L. Garden, G. M. Chaban, R. B. Gerber, D. A. Matthews, and J. F. Stanton, “Calculation of vibrational transition frequencies and intensities in water dimer: Comparison of different vibrational approaches,” The Journal of Physical Chemistry A 112, 4324–4335 (2008).
Yu and Bowman [2019] Q. Yu and J. M. Bowman, “Classical, thermostated ring polymer, and quantum VSCF/VCI calculations of IR spectra of $\text{H}_{7}\text{O}_{3}^{+}$ and $\text{H}_{9}\text{O}_{4}^{+}$ (Eigen) and comparison with experiment,” The Journal of Physical Chemistry A 123, 1399–1409 (2019).
Carpenter et al. [2020] W. B. Carpenter, Q. Yu, J. H. Hack, B. Dereka, J. M. Bowman, and A. Tokmakoff, “Decoding the 2D IR spectrum of the aqueous proton with high-level VSCF/VCI calculations,” The Journal of Chemical Physics 153, 124506 (2020).
Christiansen, Jørgensen, and Hättig [1998] O. Christiansen, P. Jørgensen, and C. Hättig, “Response functions from fourier component variational perturbation theory applied to a time-averaged quasienergy,” International Journal of Quantum Chemistry 68, 1–52 (1998).