Frozen natural orbital coupled-cluster theory: Forces and application to decomposition of nitroethane

Taube, Andrew G.; Bartlett, Rodney J.

doi:10.1063/1.2902285

The frozen natural orbital (FNO) coupled-cluster method increases the speed of coupled-cluster (CC) calculations by an order of magnitude with no consequential error along a potential energy surface. This method allows the virtual space of a correlated calculation to be reduced by about half, significantly reducing the time spent performing the coupled-cluster (CC) calculation. This paper reports the derivation and implementation of analytical gradients for FNO-CC, including all orbital relaxation for both noncanonical and semicanonical perturbed orbitals. These derivatives introduce several new orbital relaxation contributions to the CC density matrices. $FNO - CCSD (T)$ and $FNO - Λ CCSD (T)$ are applied to a test set of equilibrium structures, verifying that these methods are capable of reproducing geometries and vibrational frequencies accurately, as well as energies. Several decomposition pathways of nitroethane are investigated using $CCSD (T)$ and $Λ CCSD (T)$ with 60% of the FNO virtual orbitals in a cc-pVTZ basis, and find differences on the order of $5 kcal ∕ mol$ with reordering of the transition state energies when compared to B3LYP $6 - 311 + G (3 d f, 2 p)$ ⁠.

I. INTRODUCTION

Prediction of structures of equilibria and transition states is among the most important tasks for computational chemistry. Coupled-cluster theory has shown itself to be well-suited for determining equilibrium structures, especially in the coupled-cluster singles, doubles, and perturbative triples $[CCSD (T)]$ form.^1–3 A significant drawback of the coupled-cluster approach is the high computational scaling with respect to the size of the system. For $CCSD (T)$ ⁠, for example, the most expensive step is $O (n^{3} N^{4})$ where $n$ is the number of occupied orbitals and $N$ is the number of unoccupied orbitals and $M = N + n$ is the total number of basis functions. It is well-known, however, that standard basis sets are not optimal; one can reduce the size of the virtual (unoccupied) space without adversely affecting the numerical results. In particular, in large basis sets there are combinations of virtual orbitals that do not contribute significantly to the coupled-cluster (CC) energy. One can, therefore, reduce the computational cost by identifying and removing these irrelevant functions from the basis set. There is a long history of trying to generate such spaces for configuration interaction^4–11 and many-body perturbation theory^12–14 (MBPT) and, more recently for coupled-cluster theory.^15–22 Perhaps the most powerful method of doing so is to use frozen natural orbitals (FNOs).^{18,20,23–28} These orbitals use information from an approximate one-particle reduced density matrix to choose the best subset of one-particle orbitals within which to perform a correlated calculation. When using FNOs based on the MBPT(2) density matrix,^29,30 this truncated orbital set has been shown to be surprisingly effective at truncating larger basis sets, allowing $\sim 50 %$ of a modified unoccupied orbital set to be removed without significant changes to ground state CC energies and density matrices.^18–20 The FNO procedure could be combined with further reductions in the underlying contracted Gaussian basis^31–33 for further savings.

The frozen natural orbitals are able to achieve this speed-up by tailoring the basis set to a molecule at a particular geometry. Therefore, when one wants to calculate forces there is a nontrivial component due to the changes in the underlying structure of the frozen natural orbitals. To account for these changes, one must introduce a set of coupled-perturbed frozen natural orbital equations, similar to the coupled-perturbed Hartree-Fock (CPHF) equations, and rearrange the terms in a computationally efficient manner, to avoid calculating many perturbation-dependent quantities.

High-energy gas phase chemistry is an area that both tests and takes advantage of computational chemistry techniques. These reactions are especially susceptible to small errors in the treatment of correlation both in energetics and structures.³⁴ To be predictive, one must use high-levels of correlation and large basis sets, which tax computational resources. Nitroethane is a prototype for the decomposition of the nitroalkane class of high-energy materials.^35,36 When compared to the better studied nitromethane,^37–39 the decomposition of nitroethane introduces one important additional pathway: The elimination of HONO.^40,41 This pathway is apparently the kinetically favored one for thermal decomposition of nitroethane, and is important in the decomposition of more complicated materials, such as 1,3,5-trinitrohexahydro-1,3,5-triazine (RDX). The energetics of decomposition for nitroethane have been studied using the density functional theory (DFT) (using the B3LYP functional) (Ref. 42) but to be able to understand the relative importance of the various pathways with confidence high-level correlated calculations should be performed at the appropriately optimized stationary points. Nitroethane has five heavy atoms and, to use a large enough basis set to be definitive, one needs several hundred basis functions. To perform multiple optimizations and transition state searches with $CCSD (T)$ in that size of basis is a computational challenge; instead, we apply the FNO procedure, with analytical gradients to calculate the potential energy surface for nitroethane decomposition. We use both $CCSD (T)$ and the more recent $Λ CCSD (T)$ ⁠,^39,43–45 which, by using information from the left-hand CCSD eigenvector, improves the description of bond breaking.

II. THEORY

A. Frozen natural orbital CC

The FNO-CC method has been summarized before in Ref. 20. Below, $i$ ⁠, $j$ ⁠, and $k$ indicate occupied orbitals, $a$ ⁠, $b$ ⁠, and $c$ indicate virtual orbitals, and $p$ ⁠, $q$ ⁠, and $r$ indicate arbitrary orbitals. A set of improved virtual orbitals is generated through a series of relatively simple operations. First, the conventional HF equations are solved and then the MBPT(2) density matrix is computed in the resulting virtual orbital space. This density matrix is defined as

D_{a b}^{(2)} = \frac{1}{2} \sum_{c i j} \frac{⟨ c b ∥ i j ⟩ ⟨ i j ∥ c a ⟩}{ϵ_{i j}^{c b} ϵ_{i j}^{c a}},

(1)

where the denominator

ϵ_{i j}^{a b} = f_{i i} + f_{j j} - f_{a a} - f_{b b}

(2)

is composed of diagonal Fock matrix elements $f_{p p}$ ⁠. The density matrix is diagonalized yielding a set of natural orbitals whose occupied space has been frozen to the original Hartree-Fock orbitals (hence the name “frozen natural orbitals”). Associated with each natural orbital is an approximate occupation number. Then, based on keeping those orbitals with highest occupation, the virtual space is partitioned into two subspaces: A set of kept orbitals and a set of dropped orbitals. The Fock matrix is formed in these new orbitals, and is separately diagonalized in each of the two subspaces. Therefore, at the end of the process, one has three sets of orbitals: Canonical Hartree-Fock occupied orbitals, a kept set of Hartree-Fock virtual orbitals that are canonical among themselves, and a dropped set of Hartree-Fock virtual orbitals that are also canonical among themselves. Orbitals in the original Hartree-Fock basis will be uncapitalized, while orbitals after the FNO transformation will be denoted by capitals. Kept virtuals are indicated by $A^{'}$ ⁠, $B^{'}$ ⁠, and $C^{'}$ ⁠, dropped virtuals are $A^{″}$ ⁠, $B^{″}$ ⁠, and $C^{″}$ ⁠, and arbitrary virtuals (the union of the kept and dropped orbitals) are $A$ ⁠, $B$ ⁠, and $C$ ⁠.

Hartree-Fock orbitals in a given basis set are completely defined by the following Brillouin condition:

f_{I A^{'}} = f_{I A^{″}} = 0 .

(3)

Similarly, at the end of the FNO procedure, the additional “Brillouin” condition is

D_{A^{'} A^{″}}^{(2)} = 0 .

(4)

This equation is similar to the Hartree-Fock condition for noncanonical orbitals as the off-diagonal elements $D_{A^{'} B^{'}}^{(2)}$ and $D_{A^{″} B^{″}}^{(2)}$ can be nonzero. The combination of the Brillouin condition and the equivalent FNO condition will be used to calculate gradients.

The overall set of FNO molecular orbital coefficients can be expressed in the following series of equalities. The elements $V_{μ P}$ are the overall transformations, while $C_{μ p}$ is the original Hartree-Fock transformation and $U_{a B}$ is the additional FNO transformation in the virtual space, where $μ, ν, \dots$ are atomic orbitals,

V_{μ I} = C_{μ i},

(5a)

V_{μ A} = \sum_{b} C_{μ b} U_{b A} .

(5b)

The final molecular orbitals are defined by

∣ P ⟩ = \sum_{μ} V_{μ P} ∣ μ ⟩ .

(6)

B. Gradients

A general expression for CC gradients is^46–48

\frac{\partial E_{CC}}{\partial χ} = \sum_{P^{'} Q^{'}} γ_{P^{'} Q^{'}} \frac{\partial f_{P^{'} Q^{'}}}{\partial χ} + \sum_{P^{'} Q^{'} R^{'} S^{'}} Γ_{P^{'} Q^{'}, R^{'} S^{'}} \frac{\partial ⟨ P^{'} Q^{'} ∥ R^{'} S^{'} ⟩}{\partial χ},

(7)

where $P^{'}$ ⁠, $Q^{'}$ ⁠, $R^{'}$ ⁠, and $S^{'}$ run over all correlated orbitals. $γ_{P^{'} Q^{'}}$ and $Γ_{P^{'} Q^{'}, R^{'} S^{'}}$ are the one- and two-particle coupled-cluster response density matrices in the active space. The derivatives $\partial f_{P^{'} Q^{'}} ∕ \partial χ$ and $\partial ⟨ P^{'} Q^{'} ∥ R^{'} S^{'} ⟩ ∕ \partial χ$ are total derivatives of the molecular orbital Fock operator and two-electron integrals. These total derivatives can be separated into a piece due to the atomic orbitals and a piece due to the molecular orbital coefficients; to calculate the gradient efficiently, it is necessary to distinguish between these two.

The derivative of an active FNO orbital with respect to an external perturbation is

\frac{\partial ∣ P^{'} ⟩}{\partial χ} = \sum_{μ} \frac{\partial V_{μ P^{'}}}{\partial χ} ∣ μ ⟩ + \sum_{μ} V_{μ P^{'}} \frac{\partial ∣ μ ⟩}{\partial χ} .

(8)

Focusing on the first (molecular orbital) term one can parametrize the response as

\frac{\partial V_{μ P^{'}}}{\partial χ} = \sum_{Q} V_{μ Q} V_{Q P^{'}}^{χ} .

(9)

The coefficient $V_{Q P^{'}}^{χ}$ is the equivalent of a coupled-perturbed Hartree-Fock coefficient for the FNOs, which we will refer to as a CPFNO coefficient. It is important to note that the CPFNO coefficients have contributions from all orbitals, including those dropped during the FNO procedure.

The atomic orbital piece from Eq. (8) contributes to several perturbed integrals that are transformed into the FNO basis (here and below we are assuming real orbitals),

S_{P Q}^{χ} = \sum_{μ ν} V_{μ P} V_{ν Q} \frac{\partial S_{μ ν}}{\partial χ} = \sum_{μ ν} V_{μ P} V_{ν Q} \frac{\partial ⟨ μ ∣ ν ⟩}{\partial χ},

(10)

h_{P Q}^{χ} = \sum_{μ ν} V_{μ P} V_{ν Q} \frac{\partial h_{μ ν}}{\partial χ} = \sum_{μ ν} V_{μ P} V_{ν Q} \frac{\partial ⟨ μ ∣ h ∣ ν ⟩}{\partial χ},

(11)

{⟨ P Q ∥ R S ⟩}^{χ} = \sum_{μ ν λ σ} V_{μ P} V_{ν Q} V_{λ R} V_{σ S} \frac{\partial ⟨ μ ν ∥ λ σ ⟩}{\partial χ},

(12)

where the atomic orbital part of the Fock matrix derivative is

f_{P Q}^{(χ)} = h_{P Q}^{χ} + \sum_{i} {⟨ P I ∥ Q I ⟩}^{χ} .

(13)

Then the full partial derivatives of the Fock matrix and the two-electron integrals can be written as

\frac{\partial f_{P Q}}{\partial χ} = f_{P Q}^{(χ)} + \sum_{R} [f_{R Q} V_{R P}^{χ} + f_{R P} V_{R Q}^{χ}] + \sum_{R I} [⟨ P R ∥ Q I ⟩ + ⟨ P I ∥ Q R ⟩] V_{R I}^{χ},

(14)

\frac{\partial ⟨ P Q ∥ R S ⟩}{\partial χ} = {⟨ P Q ∥ R S ⟩}^{χ} + \sum_{U} [V_{U P}^{χ} ⟨ U Q ∥ R S ⟩ + V_{U Q}^{χ} ⟨ P U ∥ R S ⟩ + V_{U R}^{χ} ⟨ P Q ∥ U S ⟩ + V_{U S}^{χ} ⟨ P Q ∥ R U ⟩] .

(15)

Substituting these definitions into Eq. (7),

\frac{\partial E_{CC}}{\partial χ} = \sum_{P^{'} Q^{'}} γ_{P^{'} Q^{'}} f_{P^{'} Q^{'}}^{(χ)} + \sum_{P^{'} Q^{'} R^{'} S^{'}} Γ_{P^{'} Q^{'}, R^{'} S^{'}} {⟨ P^{'} Q^{'} ∥ R^{'} S^{'} ⟩}^{χ} - 2 \sum_{P P^{'}} I_{P P^{'}}^{″} V_{P P^{'}}^{χ},

(16)

where the intermediate matrix $I_{P P^{'}}^{″}$ is

I_{P P^{'}}^{″} = - \frac{1}{2} [\sum_{Q^{'}} (γ_{P^{'} Q^{'}} + γ_{Q^{'} P^{'}}) f_{P Q^{'}} + \sum_{Q^{'} R^{'} S^{'}} (⟨ P Q^{'} ∥ R^{'} S^{'} ⟩ Γ_{P^{'} Q^{'}, R^{'} S^{'}} + ⟨ Q^{'} P ∥ R^{'} S^{'} ⟩ Γ_{Q^{'} P^{'}, R^{'} S^{'}} + ⟨ Q^{'} R^{'} ∥ P S^{'} ⟩ Γ_{Q^{'} R^{'}, P^{'} S^{'}} + ⟨ Q^{'} R^{'} ∥ S^{'} P ⟩ Γ_{Q^{'} R^{'}, S^{'} P^{'}}) + \sum_{Q^{'} R^{'}} (⟨ Q^{'} P ∥ R^{'} P^{'} ⟩ + ⟨ Q^{'} P^{'} ∥ R^{'} P ⟩) γ_{Q^{'} R^{'}} δ_{P^{'}, occ}],

(17)

with $δ_{P^{'}, occ} = 0$ if $P^{'}$ is not an occupied orbital and $= 1$ if $P^{'}$ is occupied. Unlike the energy, the derivative depends on orbitals that are dropped by the FNO procedure due to the presence of terms such as $f_{P Q^{'}}$ and $⟨ P Q^{'} ∥ R^{'} S^{'} ⟩$ ⁠.

By requiring orthonormality of the perturbed orbitals, the CPFNO coefficients satisfy

V_{Q P}^{χ} + S_{P Q}^{χ} + V_{P Q}^{χ} = 0 .

(18)

The perturbed integral $S_{P Q}^{X}$ is known, so one can solve for $V_{Q P}^{χ}$ ⁠,

V_{Q P}^{χ} = - S_{P Q}^{χ} - V_{P Q}^{χ}, V_{P P}^{χ} = - \frac{1}{2} S_{P P}^{χ} .

(19)

Therefore, there are only independent equations for $P > Q$ ⁠. Expanding the last term of Eq. (16),

- 2 \sum_{P P^{'}} I_{P P^{'}}^{″} V_{P P^{'}}^{χ} = 2 \sum_{P > P^{'}} X_{P P^{'}} V_{P P^{'}}^{χ} + \sum_{P P^{'}} I_{P P^{'}}^{'} S_{P P^{'}}^{χ},

(20)

where

X_{P P^{'}} = I_{P^{'} P}^{'} - I_{P P^{'}}^{'},

(21)

I_{P P^{'}}^{'} = {\begin{matrix} I_{P P^{'}}^{″} & for P ⩽ P^{'}, \\ I_{P^{'} P}^{″} & for P > P^{'} . \end{matrix}

(22)

Now, we must address how to calculate the CPFNO coefficients. As is the case for CPHF,⁴⁹ in CC theory their direct calculation is avoided by using the Dalgarno-Stewart interchange theorem,^50,51 and sometimes called the $z$ -vector method for CPHF.^52,53 The governing equations of the FNOs are those expressed in Eqs. (3) and (4). Differentiating these equations, one has the requirements that

\frac{\partial f_{A^{'} I}}{\partial χ} = \frac{\partial f_{A^{″} I}}{\partial χ} = 0,

(23a)

\frac{\partial D_{A^{″} A^{'}}^{(2)}}{\partial χ} = 0 .

(23b)

Taking into account that the overall molecular orbital matrix

V

is a composite of the Hartree-Fock molecular orbital matrix

C

and the FNO transformation matrix

U

⁠, differentiating Eq. ,

\frac{\partial V_{μ I}}{\partial χ} = \frac{\partial C_{μ i}}{\partial χ},

(24a)

\frac{\partial V_{μ A}}{\partial χ} = \sum_{b} [\frac{\partial C_{μ b}}{\partial χ} U_{b A} + C_{μ b} \frac{\partial U_{b A}}{\partial χ}] .

(24b)

One can use a similar parametrization of the responses of the

C

and

U

matrices as was used for

V

⁠,

\frac{\partial C_{μ p}}{\partial χ} = \sum_{q} C_{μ q} C_{q p}^{χ},

(25a)

\frac{\partial U_{a B}}{\partial χ} = \sum_{C} U_{a C} U_{C B}^{χ} .

(25b)

Note that, by construction, the matrix

U

can be expanded purely within the virtual space. Relating the CPFNO coefficients

V_{P Q}^{χ}

to

C_{p q}^{χ}

and

U_{A B}^{χ}

⁠,

V_{I J}^{χ} = C_{i j}^{χ},

(26a)

V_{A I}^{χ} = \sum_{b} U_{b A} C_{b i}^{χ},

(26b)

V_{I A}^{χ} = \sum_{b} C_{i b}^{χ} U_{b A},

(26c)

V_{A B}^{χ} = U_{A B}^{χ} + \sum_{c d} U_{c B} C_{c d}^{χ} U_{d A} .

(26d)

The coefficients

C_{p q}^{χ}

are CPHF coefficients, whose solutions will be implicit; they obey an orthonormality condition. The perturbed FNO transformation obeys a slightly different condition because

U

acts in the orthonormal Hartree-Fock basis,

C_{q p}^{χ} + S_{p q}^{χ} + C_{p q}^{χ} = 0,

(27a)

U_{B A}^{χ} + U_{A B}^{χ} = 0 .

(27b)

The relations expressed in Eq. can be rearranged to reduce the number of independent variables,

C_{q p}^{χ} = - S_{p q}^{χ} - C_{p q}^{χ}, C_{p p}^{χ} = - \frac{1}{2} S_{p p}^{χ},

(28a)

U_{B A}^{χ} = - U_{A B}^{χ}, U_{A A}^{χ} = 0 .

(28b)

For all choices of perturbed FNO orbitals, choose the underlying CPHF coefficients in the virtual space to be noncanonical (as will be discussed further below), to define

C_{a b}^{χ} = - \frac{1}{2} S_{a b}^{χ} .

(29)

This choice then can be inserted into Eq. (26d), yielding

V_{A B}^{χ} = U_{A B}^{χ} - \frac{1}{2} S_{A B}^{χ} .

(30)

The definition of $C_{i j}^{χ}$ will be dependent on choosing canonical or noncanonical perturbed orbitals. Satisfying the Hartree-Fock perturbed Brillouin conditions is unchanged due to the FNO procedure, because the FNOs are still a (noncanonical) set of Hartree-Fock orbitals. Because the Brillouin condition is still satisfied, we do not need to include single excitation contributions to the MBPT(2) density matrix. On the other hand, focusing on the perturbed density matrix reveals some additional complexities.

The form of the density matrix illustrated in Eq. (1) only holds for canonical Hartree-Fock orbitals. Therefore, to directly use that equation to derive a perturbed density matrix, as is needed to impose the CPFNO condition [Eq. (23b)] one needs to require canonical perturbed underlying. Hartree-Fock orbitals. This restriction would necessitate solving the CPHF equations before constructing the perturbed density matrix for each perturbation. We instead construct the response independently of solving the CPHF coefficients by first working in the original Hartree-Fock basis (that defined by the coefficients $C$ ⁠). The general form of the MBPT(2) density matrix for noncanonical Hartree-Fock orbitals is

D_{a b}^{(2)} = \frac{1}{2} \sum_{i j c} t_{i j}^{a b (1)} t_{i j}^{c a (1)},

(31)

where the first-order $T_{2}$ amplitudes satisfy the following equation:

ϵ_{i j}^{a b} t_{i j}^{a b (1)} = ⟨ a b ∥ i j ⟩ + P_{-} (a b) \sum_{c \neq a} f_{a c} t_{i j}^{c b (1)} - P_{-} (i j) \sum_{k \neq i} f_{i k} t_{k j}^{a b (1)},

(32)

with $P_{\pm} (p q) = 1 \pm P_{p q}$ where $P_{p q}$ interchanges orbitals $p$ and $q$ ⁠. In the case that the Hartree-Fock orbitals are canonical, the last two terms vanish, allowing the following solution:

t_{i j}^{a b (1)} = \frac{⟨ a b ∥ i j ⟩}{ϵ_{i j}^{a b}},

(33)

which, when inserted into the general expression for the density matrix, returns the original result from Eq. (1). Introducing an external perturbation yields

\frac{\partial D_{a b}^{(2)}}{\partial χ} = \frac{1}{2} \sum_{i j c} [\frac{\partial t_{i j}^{c b (1)}}{\partial χ} t_{i j}^{c a (1)} + t_{i j}^{c b (1)} \frac{\partial t_{i j}^{c a (1)}}{\partial χ}] = \frac{1}{2} P_{+} (a b) \sum_{i j c} \frac{\partial t_{i j}^{c b (1)}}{\partial χ} t_{i j}^{c a (1)},

(34)

where the perturbed $T$ amplitudes are defined by the following perturbed amplitude equation:⁵⁴

\frac{\partial ϵ_{i j}^{a b}}{\partial χ} t_{i j}^{a b (1)} + ϵ_{i j}^{a b} \frac{\partial t_{i j}^{a b (1)}}{\partial χ} = \frac{\partial ⟨ a b ∥ i j ⟩}{\partial χ} + P_{-} (a b) \sum_{c \neq a} (\frac{\partial f_{a c}}{\partial χ} t_{i j}^{c b (1)} + f_{a c} \frac{\partial t_{i j}^{c b (1)}}{\partial χ}) - P_{-} (i j) \sum_{k \neq i} (\frac{\partial f_{i k}}{\partial χ} t_{k j}^{a b (1)} + f_{i k} \frac{\partial t_{k j}^{a b (1)}}{\partial χ}) .

(35)

Using the fact that the underlying unperturbed orbitals are canonical (even if the perturbed orbitals are not), one can simplify the amplitude equation to

\frac{\partial ϵ_{i j}^{a b}}{\partial χ} t_{i j}^{a b (1)} + ϵ_{i j}^{a b} \frac{\partial t_{i j}^{a b (1)}}{\partial χ} = \frac{\partial ⟨ a b ∥ i j ⟩}{\partial χ} + P_{-} (a b) \sum_{c \neq a} \frac{\partial f_{a c}}{\partial χ} t_{i j}^{c b (1)} - P_{-} (i j) \sum_{k \neq i} \frac{\partial f_{i k}}{\partial χ} t_{k j}^{a b (1)} .

(36)

The unperturbed $T$ amplitudes are known (since they correspond to canonical orbitals), leading to the following final expression:

ϵ_{i j}^{a b} \frac{\partial t_{i j}^{a b (1)}}{\partial χ} = \frac{\partial ⟨ a b ∥ i j ⟩}{\partial χ} + P_{-} (a b) \sum_{c} \frac{\partial f_{a c}}{\partial χ} \frac{⟨ c b ∥ i j ⟩}{ϵ_{i j}^{c b}} - P_{-} (i j) \sum_{k} \frac{\partial f_{i k}}{\partial χ} \frac{⟨ a b ∥ k j ⟩}{ϵ_{k j}^{a b}} .

(37)

This set of equations depends on the occupied-virtual block of CPHF coefficients through the perturbed two-electron integrals and the perturbed Fock operator, as can be seen by expanding the full perturbed density matrix,

\frac{\partial D_{a b}^{(2)}}{\partial χ} = P_{+} (a b) [\frac{1}{2} \sum_{i j c} \frac{\partial ⟨ c b ∥ i j ⟩}{\partial χ} \frac{⟨ i j ∥ c a ⟩}{ϵ_{i j}^{c b} ϵ_{i j}^{c a}} - \sum_{i j k c} \frac{\partial f_{i k}}{\partial χ} \frac{⟨ c b ∥ k j ⟩ ⟨ i j ∥ c a ⟩}{ϵ_{k j}^{c b} ϵ_{i j}^{c b} ϵ_{i j}^{c a}} + \frac{1}{2} \sum_{i j c d} (\frac{\partial f_{c d}}{\partial χ} \frac{⟨ d b ∥ i j ⟩ ⟨ i j ∥ c a ⟩}{ϵ_{i j}^{c b} ϵ_{i j}^{d b} ϵ_{i j}^{c a}} + \frac{\partial f_{b d}}{\partial χ} \frac{⟨ c b ∥ i j ⟩ ⟨ i j ∥ c a ⟩}{ϵ_{i j}^{c b} ϵ_{i j}^{c d} ϵ_{i j}^{c a}})] .

(38)

Expanding the integral derivatives,

\frac{\partial D_{a b}^{(2)}}{\partial χ} = P_{+} (a b) [\sum_{c d} (Y_{a b, c d} S_{c d}^{χ} + M_{a b, c d} f_{c d}^{(χ)}) + \sum_{i > j} (Y_{a b, i j} - Y_{a b, j i}) C_{i j}^{χ} + \sum_{i j} M_{a b, i j} f_{i j}^{(χ)} - \sum_{i > j} Y_{a b, j i} S_{i j}^{χ} + \sum_{i c} (Y_{a b, c i} C_{c i}^{χ} + M_{a b, c i} S_{c i}^{χ})] + D_{a b}^{(2), χ},

(39)

where the intermediate quantities $M$ and $Y$ are defined in Table I and the perturbed density matrix $D_{a b}^{(2), χ}$ is defined as

D_{a b}^{(2), χ} = \frac{1}{2} P_{+} (a b) \sum_{i j c} \frac{{⟨ c b ∥ i j ⟩}^{χ} ⟨ i j ∥ c a ⟩}{ϵ_{i j}^{c b} ϵ_{i j}^{c a}} .

(40)

This term will be further discussed below.

Table I.

Intermediate quantities for FNO gradients.


$M_{a b, i j} = - \sum_{k c} \frac{⟨ c b ∥ j k ⟩ ⟨ i k ∥ c a ⟩}{ϵ_{j k}^{c b} ϵ_{i k}^{c b} ϵ_{i k}^{c a}}$
$M_{a b, c i} = - \frac{1}{2} \sum_{j k} \frac{⟨ i b ∥ j k ⟩ ⟨ j k ∥ c a ⟩}{ϵ_{j k}^{c b} ϵ_{j k}^{c a}} - \frac{1}{2} δ_{b c} \sum_{j k d} \frac{⟨ d i ∥ j k ⟩ ⟨ j k ∥ d a ⟩}{ϵ_{j k}^{d b} ϵ_{j k}^{d a}}$
$M_{a b, c d} = \frac{1}{2} \sum_{i j} \frac{⟨ d b ∥ i j ⟩ ⟨ i j ∥ c a ⟩}{ϵ_{i j}^{d b} ϵ_{i j}^{c b} ϵ_{i j}^{c a}} + \frac{1}{2} δ_{b c} \sum_{i j e} \frac{⟨ e d ∥ i j ⟩ ⟨ i j ∥ e a ⟩}{ϵ_{i j}^{e b} ϵ_{i j}^{e d} ϵ_{i j}^{e a}}$
$Y_{a b, i j} = (ϵ_{i} + ϵ_{j}) M_{a b, i j} + 2 \sum_{k l} M_{a b, k l} ⟨ k i ∥ l j ⟩ + 2 \sum_{c d} M_{a b, c d} ⟨ c i ∥ d j ⟩$
$Y_{a b, c d} = - \frac{1}{4} \sum_{i j} \frac{⟨ d b ∥ i j ⟩ ⟨ i j ∥ c a ⟩}{ϵ_{i j}^{c b} ϵ_{i j}^{c a}} - \frac{1}{2} (ϵ_{c} + ϵ_{d}) M_{a b, c d}$
$Y_{a b, c i} = M_{a b, c i} + \sum_{j k} M_{a b, j k} (⟨ j c ∥ k i ⟩ + ⟨ j i ∥ k c ⟩) + \sum_{d e} M_{a b, d c} (⟨ d c ∥ e i ⟩ + ⟨ d i ∥ e c ⟩)$

These equations have been derived in the original Hartree-Fock basis, but the CPFNO equation is in the FNO basis. To transform the results, one uses

D_{A^{″} A^{'}}^{(2)} = \sum_{b c} U_{b A^{″}} D_{b c}^{(2)} U_{c A^{'}} .

(41)

Differentiating this expression,

\frac{\partial D_{A^{″} A^{'}}^{(2)}}{\partial χ} = \sum_{b c} [\frac{\partial U_{b A^{″}}}{\partial χ} D_{b c}^{(2)} U_{c A^{'}} + U_{b A^{″}} \frac{\partial D_{b c}^{(2)}}{\partial χ} U_{c A^{'}} + U_{b A^{″}} D_{b c}^{(2)} \frac{\partial U_{c A^{'}}}{\partial χ}] .

(42)

The perturbed FNO density can therefore be written, using $D_{B^{″} B^{'}}^{(2)} = 0$ and the relationship expressed in Eq. (30),

\frac{\partial D_{A^{″} A^{'}}^{(2)}}{\partial χ} = - \sum_{B^{'}} D_{A^{'} B^{'}}^{(2)} V_{A^{″} B^{'}}^{χ} + \sum_{B^{″}} D_{A^{″} B^{″}}^{(2)} V_{B^{″} A^{'}}^{χ} - \frac{1}{2} \sum_{B^{'}} D_{A^{'} B^{'}}^{(2)} S_{A^{″} B^{'}}^{χ} + \frac{1}{2} \sum_{B^{″}} D_{A^{″} B^{″}}^{(2)} S_{B^{″} A^{'}}^{χ} + P_{+} (A^{″} A^{'}) [\sum_{I B} Y_{A^{″} A^{'}, B I} V_{B I}^{χ} + \sum_{I > J} (Y_{A^{″} A^{'}, I J} - Y_{A^{″} A^{'}, J I}) V_{I J}^{χ} + \sum_{B C} (Y_{A^{″} A^{'}, B C} S_{B C}^{χ} + M_{A^{″} A^{'}, B C} f_{B C}^{(χ)}) + \sum_{I J} M_{A^{″} A^{'}, I J} f_{I J}^{(χ)} - \sum_{I > J} Y_{A^{″} A^{'}, J I} S_{I J}^{χ} + \sum_{I B} M_{A^{″} A^{'}, B I} S_{B I}^{χ}] + D_{A^{″} A^{'}}^{(2), χ},

(43)

where the quantities $M$ and $Y$ have been transformed to the FNO basis. The perturbed quantities have therefore been completely separated from the CPFNO coefficients, which allows for a perturbation-independent solution of the CPFNO equations. To go further, one must choose between noncanonical and canonical perturbed FNOs.

1. Noncanonical perturbed orbitals

For the choice of noncanonical perturbed orbitals, we are free to define

V_{I J}^{χ} = - \frac{1}{2} S_{I J}^{χ},

(44a)

V_{A^{'} B^{'}}^{χ} = - \frac{1}{2} S_{A^{'} B^{'}}^{χ},

(44b)

V_{A^{″} B^{″}}^{χ} = - \frac{1}{2} S_{A^{″} B^{″}}^{χ} .

(44c)

The only terms that need to be solved are

V_{A I}^{χ}

and

V_{A^{″} B^{'}}^{χ}

⁠.

Then CPFNO equations in matrix form are

[\begin{matrix} A_{vo, vo} & 0 \\ A_{v^{″} v^{'}, vo} & A_{v^{″} v^{'}, v^{″} v^{'}} \end{matrix}] [\begin{matrix} V_{vo}^{χ} \\ V_{v^{″} v^{'}}^{χ} \end{matrix}] = [\begin{matrix} B_{vo}^{χ} \\ B_{v^{″} v^{'}}^{χ} \end{matrix}],

(45)

where

B_{A I}^{χ} = - f_{A I}^{(χ)} + S_{A I}^{χ} f_{I I} + \sum_{J K} S_{J K}^{χ} ⟨ A J ∥ I K ⟩,

(46)

B_{A^{″} A^{'}}^{χ} = P_{+} (A^{″} A^{'}) [\sum_{B C} (Y_{A^{″} A^{'}, B C} S_{B C}^{χ} + M_{A^{″} A^{'}, B C} f_{B C}^{(χ)}) + \sum_{I J} (M_{A^{″} A^{'}, I J} f_{I J}^{(χ)} - \frac{1}{2} Y_{A^{″} A^{'}, J I} S_{I J}^{χ}) + \sum_{I B} M_{A^{″} A^{'}, B I} S_{B I}^{χ}] - \frac{1}{2} \sum_{B^{'}} S_{A^{″} B^{'}}^{χ} D_{B^{'} A^{'}}^{(2)} + \frac{1}{2} \sum_{B^{″}} D_{A^{″} B^{″}}^{(2)} S_{B^{″} A^{'}}^{χ} + D_{A^{″} A^{'}}^{(2), χ},

(47)

and the elements of matrix $A$ are in Table II.

Table II.

Elements of the noncanonial CPFNO $A$ matrix.


$A_{A I, B J} = ⟨ A B ∥ I J ⟩ + ⟨ A J ∥ I B ⟩ - δ_{I J} f_{A B} + δ_{I J} δ_{A B} f_{I I}$
$A_{A^{″} A^{'}, B J} = Y_{A^{″} A^{'}, B J} + Y_{A^{'} A^{″}, B J}$
$A_{A^{″} A^{'}, B^{″}, B^{'}} = δ_{A^{″} B^{″}} D_{A^{'} B^{'}}^{(2)} - δ_{A^{'} B^{'}} D_{A^{″} B^{″}}^{(2)}$

The interchange theorem can be written

\sum_{P > Q} X_{P Q} V_{P Q}^{χ} = X^{T} V^{χ} = X^{T} A^{- 1} B^{χ} \equiv - D^{(or) T} B^{χ} = - \sum_{P > Q} D_{P Q}^{(or)} B_{P Q}^{χ} .

(48)

Therefore, one can solve the perturbation-independent equation below, instead of Eq. (45),

[\begin{matrix} A_{vo, vo}^{T} & A_{vo, v^{″} v^{'}}^{T} \\ 0 & A_{v^{″} v^{'}, v^{″}, v^{'}}^{T} \end{matrix}] [\begin{matrix} D_{vo}^{(or)} \\ D_{v^{″} v^{'}}^{(or)} \end{matrix}] = [\begin{matrix} - X_{vo} \\ - X_{v^{″} v^{'}} \end{matrix}] .

(49)

The second of these equations, which determines the orbital response of the uncorrelated molecular orbitals, can be solved independently of the first, using a standard linear equation solver. Substituting this result into the equation for the virtual-occupied block of the orbital response contribution to the density matrix,

\sum_{J B} D_{B J}^{(or)} [⟨ B A ∥ J I ⟩ + ⟨ B I ∥ J A ⟩ + (f_{I I} δ_{A B} - f_{A B}) δ_{I J}] = - {\tilde{X}}_{A I},

(50)

{\tilde{X}}_{A I} = X_{A I} + \sum_{B^{″} B^{'}} D_{B^{″} B^{'}}^{(or)} [Y_{B^{″} B^{'}, A I} + Y_{B^{'} B^{″}, A I}] .

(51)

This modified $z$ -vector equation can then be solved by the standard method.⁵² Using these orbital relaxation components, one can form the final full density matrices,

D_{I J} = γ_{I J} - \sum_{A^{″} A^{'}} (M_{A^{″} A^{'}, I J} + M_{A^{'} A^{″}, I J}) D_{A^{″} A^{'}}^{(or)},

(52a)

D_{A I} = γ_{A I} + D_{A I}^{(or)},

(52b)

D_{A B} = γ_{B C} - \sum_{C^{″} C^{'}} (M_{C^{″} C^{'}, A B} + M_{C^{'} C^{″}, A B}) D_{C^{″} C^{'}}^{(or)},

(52c)

and intermediate matrices

I_{I J} = I_{I J}^{'} - \sum_{A K} D_{A K}^{(or)} ⟨ A I ∥ K J ⟩ + \frac{1}{2} \sum_{A^{″} A^{'}} (Y_{A^{″} A^{'}, J I} + Y_{A^{'} A^{″}, J I}) D_{A^{″} A^{'}}^{(or)},

(53a)

I_{A I} = I_{A I}^{'} - f_{I I} D_{A I}^{(or)} - \sum_{B^{″} B^{'}} (M_{B^{″} B^{'}, A I} - M_{B^{'} B^{″}, A I}) D_{B^{″} B^{'}}^{(or)},

(53b)

I_{A^{'} B^{'}} = I_{A^{'} B^{'}}^{'} - \sum_{C^{″} C^{'}} (Y_{C^{″} C^{'}, A^{'} B^{'}} + Y_{C^{'} C^{″}, A^{'} B^{'}}) D_{C^{″} C^{'}}^{(or)},

(53c)

I_{A^{″} A^{'}} = I_{A^{″} A^{'}}^{'} + \frac{1}{2} \sum_{B^{'}} D_{A^{″} B^{'}}^{(or)} D_{B^{'} A^{'}}^{(2)} - \frac{1}{2} \sum_{B^{″}} D_{A^{″} B^{″}}^{(2)} D_{B^{″} A^{'}}^{(or)} - \sum_{B^{″} B^{'}} (Y_{B^{″} B^{'}, A^{″}, A^{'}} + Y_{B^{'} B^{″}, A^{″} A^{'}}) D_{B^{″} B^{'}}^{(or)},

(53d)

I_{A^{″} B^{″}} = - \sum_{C^{″} C^{'}} (Y_{C^{″} C^{'}, A^{″} B^{″}} + Y_{C^{'} C^{″}, A^{″} B^{″}}) D_{C^{″} C^{'}}^{(or)} .

(53e)

The final term left to address is

D_{A^{″} A^{'}}^{(2), χ}

from

B_{A^{″} A^{'}}^{χ}

⁠. This term is only simply expressed [as in Eq. ] when the orbitals are those from canonical Hartree-Fock. Therefore, to completely separate the perturbed orbital contribution requires additional integral transforms. First, one must back-transform the orbital response density matrix

D_{A^{″} A^{'}}^{(or)}

to the Hartree-Fock basis,

D_{a b}^{(or)} = \sum_{C^{″} C^{'}} U_{a C^{″}} D_{C^{″} C^{'}}^{(or)} U_{b C^{'}} .

(54)

This object can now be contracted with the perturbation independent pieces of $D_{a^{″} a^{'}}^{(2)}$ ⁠,

G_{a b, i j} = \sum_{c} D_{c b}^{(or)} \frac{⟨ i j ∥ a c ⟩}{ϵ_{i j}^{b c} ϵ_{i j}^{a c}} .

(55)

This term can only be formed in the Hartree-Fock basis, leading to the following final expression for noncanonical perturbed orbital gradients:

\frac{\partial E_{CC}}{\partial χ} = \sum_{P Q} D_{P Q} f_{P Q}^{(χ)} + \sum_{P Q R S} Γ_{P Q, R S} {⟨ P Q ∥ R S ⟩}^{χ} + \sum_{P Q} I_{P Q} S_{P Q}^{χ} + \sum_{a b i j} G_{a b, i j} {⟨ a b ∥ i j ⟩}^{χ} .

(56)

Because of the extra $G$ term, a separate back-transformation is necessary to write the $G_{a b, i j}$ term in the atomic orbital basis before contraction with the derivative integrals.

2. Canonical gradients

In the case of $CCSD (T)$ ⁠, it is highly advantageous to impose the condition that the perturbed orbitals remain semicanonical.^2,55 When frozen occupied or virtuals are used, the derivative is calculated using canonical perturbed orbitals as well.⁴⁸ (This condition is actually more stringent than strictly necessary; as long as mixing occurs only within the frozen and active subsets of orbitals, they do not need to be maintained canonical.) Therefore, one must formulate the FNO orbital relaxation terms in semicanonical orbitals. Supplementing the conventional Brillouin condition (and the FNO condition) are the requirements that

\frac{\partial f_{I J}}{\partial χ} = 0, \frac{\partial f_{A^{'} B^{'}}}{\partial χ} = 0 .

(57)

There is no need to impose canonicality on the uncorrelated orbitals ${A^{″}}$ because the computational advantage lies in determining the CC contribution to the density matrices, which does not involve the uncorrelated orbitals. By imposing this requirement one can no longer choose $V_{I J}^{χ} = - 1 ∕ 2 S_{I J}^{χ}$ and $V_{A^{'} B^{'}}^{χ} = - 1 ∕ 2 S_{A^{'} B^{'}}^{χ}$ ⁠. However, one can choose $V_{A^{″} B^{″}}^{χ} = - 1 ∕ 2 S_{A^{″} B^{″}}^{χ}$ ⁠, since the dropped virtuals can be noncanonical. Therefore, no iterative equations have to be solved in the uncorrelated-uncorrelated sector.

The new CPHF equations can be written in matrix form as follows:

[\begin{matrix} A_{oo, oo} & A_{oo, vo} & 0 & 0 \\ 0 & A_{vo, vo} & 0 & 0 \\ 0 & A_{v^{'} v^{'}, vo} & A_{v^{'} v^{'}, v^{'} v^{'}} & A_{v^{'} v^{'}, v^{″} v^{'}} \\ A_{v^{″} v^{'}, oo} & A_{v^{″} v^{'}, vo} & 0 & A_{v^{″} v^{'}, v^{″} v^{'}} \end{matrix}] [\begin{matrix} V_{oo}^{χ} \\ V_{vo}^{χ} \\ V_{v^{'} v^{'}}^{χ} \\ V_{v^{″} v^{'}}^{χ} \end{matrix}] = [\begin{matrix} B_{oo}^{χ} \\ B_{vo}^{χ} \\ B_{v^{'} v^{'}}^{χ} \\ B_{v^{″} v^{'}}^{χ} \end{matrix}] .

(58)

The right-hand side of this equation is given by

B_{P^{'} Q^{'}}^{χ} = - f_{P^{'} Q^{'}}^{(χ)} + S_{P^{'} Q^{'}}^{χ} f_{Q^{'} Q^{'}} + \sum_{I J} S_{I J}^{χ} ⟨ P^{'} I ∥ Q^{'} J ⟩,

(59)

B_{A^{″} A^{'}}^{χ} = P_{+} (A^{″} A^{'}) [\sum_{B C} (Y_{A^{″} A^{'}, B C} S_{B C}^{χ} + M_{A^{″} A^{'}, B C} f_{B C}^{(χ)}) + \sum_{I J} M_{A^{″} A^{'}, I J} f_{I J}^{(χ)} - \sum_{I < J} Y_{A^{″} A^{'}, J I} S_{I J}^{χ} + \sum_{I B} M_{A^{″} A^{'}, B I} S_{B I}^{χ}] - \frac{1}{2} \sum_{B^{'}} S_{A^{″} B^{'}}^{χ} D_{B^{'} A^{'}}^{(2)} + \frac{1}{2} \sum_{B^{″}} D_{A^{″} B^{″}}^{(2)} S_{B^{″} A}^{χ} + D_{A^{″} A^{'}}^{(2), χ} .

(60)

The matrix elements of $A$ are given in Table III.

Table III.

Elements of the canonical CPFNO $A$ matrix.


$A_{I J, K L} = δ_{I K} δ_{J L} (f_{I I} - f_{J J})$
$A_{I J, A K} = ⟨ I A ∥ J K ⟩ + ⟨ I K ∥ J A ⟩$
$A_{A I, B J} = δ_{I J} f_{A B} - δ_{I J} δ_{A B} f_{I I} + ⟨ A B ∥ I J ⟩ + ⟨ A J ∥ I B ⟩$
$A_{A^{'} B^{'}, C I} = ⟨ A^{'} C ∥ B^{'} I ⟩ + ⟨ A^{'} I ∥ B^{'} C ⟩$
$A_{A^{'} B^{'}, C^{'} D^{'}} = δ_{A^{'} C^{'}} δ_{B^{'} D^{'}} (f_{A^{'} A^{'}} - f_{B^{'} B^{'}})$
$A_{A^{'} B^{'}, A^{″} C^{″}} = δ_{B^{'} C^{'}} f_{A^{'} A^{″}} + δ_{A^{'} C^{'}} f_{A^{″} B^{'}}$
$A_{A^{″} A^{'}, I J} = Y_{A^{″} A^{'}, I J} + Y_{A^{″} A^{'}, I J} - Y_{A^{'} A^{″}, J I} - Y_{A^{'} A^{″}, J I}$
$A_{A^{″} A^{'}, B J} = Y_{A^{″} A^{'}, B J} + Y_{A^{'} A^{″}, B J}$
$A_{A^{″} A^{'}, B^{″} B^{'}} = δ_{A^{″} B^{″}} D_{A^{'} B^{'}}^{(2)} - δ_{A^{'} B^{'}} D_{A^{″} B^{″}}^{(2)}$

Solving the linear equation

[\begin{matrix} A_{oo, oo}^{T} & 0 & 0 & A_{oo, v^{″} v^{'}}^{T} \\ A_{vo, oo}^{T} & A_{vo, vo}^{T} & A_{vo, v^{'} v^{'}}^{T} & A_{vo, v^{″} v^{'}}^{T} \\ 0 & 0 & A_{v^{'} v^{'}, v^{'} v^{'}}^{T} & 0 \\ 0 & 0 & A_{v^{″} v^{'}, v^{'} v^{'}}^{T} & A_{v^{″} v^{'}, v^{″} v^{'}}^{T} \end{matrix}] [\begin{matrix} D_{oo}^{(or)} \\ D_{vo}^{(or)} \\ D_{v^{'} v^{'}}^{(or)} \\ D_{v^{″} v^{'}}^{(or)} \end{matrix}] = [\begin{matrix} - X_{oo} \\ - X_{vo} \\ - X_{v^{'} v^{'}} \\ - X_{v^{″} v^{'}} \end{matrix}]

(61)

yields the orbital response contribution to the overall density matrix. In this form, it is obvious that the solution for the active virtual–active virtual block does not couple to the other blocks, yielding

D_{A^{'} B^{'}}^{(or)} = \frac{- X_{A^{'} B^{'}}}{f_{A^{'} A^{'}} - f_{B^{'} B^{'}}} .

(62)

The FNO block can then be determined by inserting the new orbital relaxation terms from the active virtual block,

\sum_{B^{″} B^{'}} D_{B^{″} B^{'}}^{(or)} [δ_{A^{″} B^{″}} D_{A^{'} B^{'}}^{(2)} - δ_{A^{'} B^{'}} D_{A^{″} B^{″}}^{(2)}] = - {\tilde{X}}_{A^{″} A^{'}},

(63)

{\tilde{X}}_{A^{″} A^{'}} = X_{A^{″} A^{'}} + \sum_{B^{'}} f_{A^{″} B^{'}} D_{B^{'} A^{'}}^{(or)} .

(64)

Furthermore, after solving for the inactive-active virtual block, the orbital response for the occupied-occupied block can be solved,

\sum_{K > L} D_{K L}^{(or)} [δ_{I K} δ_{J L} (f_{I I} - f_{J J})] = - {\tilde{X}}_{I J} \Rightarrow D_{I J}^{(or)} = \frac{- {\tilde{X}}_{I J}}{f_{I I} - f_{J J}},

(65)

{\tilde{X}}_{I J} = X_{I J} + P_{-} (I J) \sum_{A^{″} A^{'}} (Y_{A^{″} A^{'}, I J} + Y_{A^{'} A^{″}, I J}) D_{A^{″} A^{'}}^{(or)} .

(66)

Finally, the response of the occupied-occupied block and the virtual-virtual blocks can be inserted into the equation for the occupied-virtual block,

\sum_{J B} D_{B J}^{(or)} [δ_{A B} δ_{I J} (f_{A A} - f_{I I}) + ⟨ B A ∥ J I ⟩ + ⟨ B I ∥ J A ⟩] = - {\tilde{X}}_{A I},

(67)

{\tilde{X}}_{A I} = X_{A I} + \frac{1}{2} \sum_{P^{'} Q^{'}} D_{P^{'} Q^{'}}^{(or)} [⟨ A P^{'} ∥ I Q^{'} ⟩ + ⟨ A Q^{'} ∥ I P^{'} ⟩] + \sum_{B^{″} B^{'}} D_{B^{″} B^{'}}^{(or)} [Y_{B^{″} B^{'}, A I} + Y_{B^{'} B^{″}, A I}] .

(68)

This equation now fits the standard form of the $z$ -vector equations.

After solving for all of the orbital response components of the density matrices, one can define the full, relaxed, density matrices via

D_{I J} = γ_{I J} + D_{I J}^{(or)} - \sum_{A^{″} A^{'}} (M_{A^{″} A^{'}, I J} + M_{A^{'} A^{″}, I J}) D_{A^{″} A^{'}}^{(or)},

(69a)

D_{A I} = γ_{A I} + D_{A I}^{(or)},

(69b)

D_{A B} = γ_{B C} + D_{A B}^{(or)} - \sum_{C^{″} C^{'}} (M_{C^{″} C^{'}, A B} + M_{C^{'} C^{″}, A B}) D_{C^{″} C^{'}}^{(or)},

(69c)

and the final intermediate matrices,

I_{I J} = I_{I J}^{'} - \sum_{P Q} ⟨ I P ∥ J Q ⟩ D_{P Q}^{(or)} - f_{J J} D_{I J}^{(or)},

(70a)

I_{A I} = I_{A I}^{'} - f_{I I} D_{A I}^{(or)} - \sum_{B^{″} B^{'}} (M_{B^{″} B^{'}, A I} - M_{B^{'} B^{″}, A I}) D_{B^{″} B^{'}}^{(or)},

(70b)

I_{A^{'} B^{'}} = I_{A^{'} B^{'}}^{'} - f_{B^{'} B^{'}} D_{A^{'} B^{'}}^{(or)} - \sum_{C^{″} C^{'}} (Y_{C^{″} C^{'}, A^{'} B^{'}} + Y_{C^{'} C^{″}, A^{'} B^{'}}) D_{C^{″} C^{'}}^{(or)},

(70c)

I_{A^{″} A^{'}} = I_{A^{″} A^{'}}^{'} - f_{A^{'} A^{'}} D_{A^{″} A^{'}}^{(or)} + \frac{1}{2} \sum_{B^{'}} D_{A^{″} B^{'}}^{(or)} D_{B^{'} A^{'}}^{(2)} - \frac{1}{2} \sum_{B^{″}} D_{A^{″} B^{″}}^{(2)} D_{B^{″} A^{'}}^{(or)} - \sum_{B^{″} B^{'}} (Y_{B^{″} B^{'}, A^{″} A^{'}} + Y_{B^{'} B^{″}, A^{″} A^{'}}) D_{B^{″} B^{'}}^{(or)},

(70d)

I_{A^{″} B^{″}} = - \sum_{C^{″} C^{'}} (Y_{C^{″} C^{'}, A^{″} B^{″}} + Y_{C^{'} C^{″}, A^{″} B^{″}}) D_{C^{″} C^{'}}^{(or)} .

(70e)

Again, it is necessary to include the

D_{A^{″} A^{'}}^{(2), χ}

term separately via a back-transformation, as is expressed in Eqs. .

A summary of the steps necessary to calculate a derivative using the FNOs is shown in Fig. 1.

FIG. 1.

View large Download slide

Steps in an FNO derivative calculation.

C. Smoothness of the potential energy surface

The FNO procedure developed here will not necessarily yield rigorously smooth potential energy surfaces (PESs). Note that the FNO truncation is performed point by point on the PES, without consideration of the connection between that point and other points on the potential energy surface. Therefore, if the structure or size of the space spanned by the correlated set of virtual orbitals changes as a function of the geometry, it is possible that the energy could change discontinuously.

To minimize the impact of discontinuities, the code recognizes orbitals that are close in occupation to the correlated orbitals. Those within a certain tolerance of the cutoff occupation are considered to be quasidegenerate and are retained. Assuming that the geometry steps are not too large, this procedure should smooth changes in the FNO structure. It should be clear that this problem is not unique to the frozen natural orbital truncation procedure, but exists for all procedures (such as localized orbitals methods) that truncate the correlation space in a geometry dependent way.^56–58 To the best of the authors’ knowledge, there is no fully satisfactory solution to this problem.

III. IMPLEMENTATION

The FNO-CC gradients have been implemented within the ACES II program system.⁵⁹ It takes advantage of real Abelian point group symmetry, and all equations are fully spin-summed and applicable to closed- or open-shells using single determinant relativistic Hartree-Fock (RHF) or spin-polarized unrestricted Hartree-Fock (UHF) references. In a FNO energy calculation, a partial integral transformation is performed before the FNO truncation, and then a full integral transformation is performed in the resultant truncated basis. This computational advantage is unachievable for gradients; instead, one must perform a full integral transformation for both the truncated and full basis sets, requiring the storage of more integrals. The correlated calculations are then performed within the truncated basis. The formation of the density matrices proceeds in two parts: First, the correlated contributions are formed within the truncated basis, then these density matrices are expanded to the full basis, and the orbital relaxation terms are calculated and included. One calculates the $M$ and $Y$ in the HF basis and then stores them in the truncated FNO basis, so that they can be added to the orbital relaxation equations. The back-transformation of the FNO density matrix and the $G_{a b, i j}$ term are performed separately. These terms are then summed before contraction with derivative integrals.

Compared to gradient calculations that do not use FNOs, the largest added expense is the necessity of calculating and storing several new intermediates of a dimension similar to that of the two-electron integrals. The computational cost is far less than the cost of the CC procedure, though, and does not change the overall scaling of the coupled-cluster, but, instead reduces its cost in applications. However, the additional storage costs could be problematic for some combinations of computer and molecule.

All gradient calculations were verified by comparing the analytical gradient expression to those obtained by numerical differentiation of the energy.

IV. RESULTS AND DISCUSSION

A. Calibration

To determine the capability of FNO truncated gradients to reproduce structures, we applied FNO $CCSD (T)$ and $Λ CCSD (T)$ to the set of well-characterized molecules from Bak et al.⁶⁰ Comparative statistics are shown in Tables IV–VII. We have chosen to show the dependence of geometrical properties versus the percentage of the virtual space retained in the truncated calculations. In some ways, this way of choosing a truncation is unsatisfying; it would be better if one were able to examine the MBPT(2) occupation numbers, and then chose proper cutoffs based on these values. However, while we have looked into this issue, we have not been able to determine any consistent truncation criterion: The occupation numbers go smoothly from high to low occupation, without any sharp changes that would indicate a place to truncate. Because the goal of the method is to reduce the computational cost of the calculation, at this point it seems better to use a truncation scheme where the speed-up can be predicted, even if it is less satisfying theoretically.

Table IV.

Comparison of optimized equilibrium bond lengths for different correlation-consistent basis sets Refs. 65 and 66 for multiple FNO truncations for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The percentage indicates what percent of the virtual space of each molecule was active. For truncated basis sets (20%–80%) errors are relative to the untruncated basis set result; for 100%, errors are relative to experiment. Averages were calculated over the set of molecules from Ref. 60. Only valence electrons were correlated. $\bar{Δ}$ is the (signed) mean error, ${\bar{Δ}}_{abs}$ is the mean absolute error, $Δ_{\max}$ is the maximum absolute error, and $Δ_{std}$ is the standard deviation. All numbers are in units of pm.

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pVDZ
20	0.02	0.66	2.30	0.80	0.05	0.63	1.93	0.74
40	$- 0.17$	0.54	1.73	0.70	$- 0.14$	0.54	1.85	0.71
60	$- 0.53$	0.54	2.43	0.63	$- 0.51$	0.52	2.3	0.60
80	$- 0.26$	0.30	3.45	0.51	$- 0.24$	0.29	3.22	0.65
100^a	1.72	1.72	4.51	0.82	1.69	1.69	4.12	0.76
cc-pVTZ
20	0.30	0.64	5.50	1.21	0.32	0.63	5.47	1.20
40	0.15	0.21	0.87	0.25	0.16	0.22	0.88	0.24
60	$- 0.03$	0.14	0.57	0.19	$- 0.02$	0.14	0.53	0.19
80	$- 0.03$	0.09	0.45	0.13	$- 0.03$	0.09	0.41	0.13
100^a	0.05	0.22	0.90	0.29	0.02	0.22	0.71	0.27
cc-pVQZ
20	0.05	0.18	0.68	0.24	0.07	0.18	0.77	0.23
40	$- 0.11$	0.18	1.11	0.27	$- 0.11$	0.17	1.01	0.26
60	0.00	0.05	0.26	0.07	0.02	0.05	0.26	0.07
80	0.00	0.03	0.10	0.04	0.00	0.02	0.10	0.03
100^a	$- 0.06$	0.13	0.71	0.19	$- 0.10$	0.14	0.71	0.19

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pVDZ
20	0.02	0.66	2.30	0.80	0.05	0.63	1.93	0.74
40	$- 0.17$	0.54	1.73	0.70	$- 0.14$	0.54	1.85	0.71
60	$- 0.53$	0.54	2.43	0.63	$- 0.51$	0.52	2.3	0.60
80	$- 0.26$	0.30	3.45	0.51	$- 0.24$	0.29	3.22	0.65
100^a	1.72	1.72	4.51	0.82	1.69	1.69	4.12	0.76
cc-pVTZ
20	0.30	0.64	5.50	1.21	0.32	0.63	5.47	1.20
40	0.15	0.21	0.87	0.25	0.16	0.22	0.88	0.24
60	$- 0.03$	0.14	0.57	0.19	$- 0.02$	0.14	0.53	0.19
80	$- 0.03$	0.09	0.45	0.13	$- 0.03$	0.09	0.41	0.13
100^a	0.05	0.22	0.90	0.29	0.02	0.22	0.71	0.27
cc-pVQZ
20	0.05	0.18	0.68	0.24	0.07	0.18	0.77	0.23
40	$- 0.11$	0.18	1.11	0.27	$- 0.11$	0.17	1.01	0.26
60	0.00	0.05	0.26	0.07	0.02	0.05	0.26	0.07
80	0.00	0.03	0.10	0.04	0.00	0.02	0.10	0.03
100^a	$- 0.06$	0.13	0.71	0.19	$- 0.10$	0.14	0.71	0.19

^a

Relative to experiment.

Table V.

Comparison of optimized equilibrium bond lengths for different correlation-consistent basis sets Refs. 65 and 66 for multiple FNO truncations for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The percentage indicates what percent of the virtual space of each molecule was active. For truncated basis sets (20%–80%) errors are relative to the untruncated basis set result; for 100%, errors are relative to experiment. Averages were calculated over the set of molecules from Ref. 60. All electrons were correlated. $\bar{Δ}$ is the (signed) mean error, ${\bar{Δ}}_{abs}$ is the mean absolute error, $Δ_{\max}$ is the maximum absolute error, and $Δ_{std}$ is the standard deviation. All numbers are in units of pm.

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pVDZ
20	$- 0.09$	0.66	2.21	0.81	$- 0.05$	0.62	1.84	0.74
40	$- 0.58$	0.63	2.96	0.73	$- 0.56$	0.60	2.86	0.69
60	$- 0.30$	0.33	2.34	0.49	$- 0.29$	0.31	2.14	0.45
80	$- 0.14$	0.17	0.99	0.24	$- 0.14$	0.17	0.98	0.24
100^a	1.66	1.66	4.42	0.80	1.63	1.63	4.03	0.74
cc-pVTZ
20	$- 0.11$	0.26	2.37	0.50	$- 0.09$	0.25	2.15	0.47
40	$- 0.24$	0.32	0.92	0.34	$- 0.23$	0.31	0.91	0.33
60	$- 0.29$	0.33	1.34	0.43	$- 0.28$	0.32	1.34	0.42
80	$- 0.09$	0.11	0.33	0.12	$- 0.09$	0.11	0.33	0.12
100^a	0.19	0.26	1.04	0.28	0.15	0.25	0.86	0.27
cc-pVQZ
20	$- 0.21$	0.25	0.71	0.23	$- 0.21$	0.25	0.70	0.22
40	$- 0.06$	0.07	0.33	0.10	$- 0.06$	0.08	0.30	0.09
60	$- 0.05$	0.09	0.52	0.13	$- 0.04$	0.09	0.49	0.12
80	0.00	0.01	0.05	0.02	0.00	0.01	0.05	0.02
100^a	$- 0.01$	0.09	0.63	0.17	$- 0.06$	0.10	0.64	0.18

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pVDZ
20	$- 0.09$	0.66	2.21	0.81	$- 0.05$	0.62	1.84	0.74
40	$- 0.58$	0.63	2.96	0.73	$- 0.56$	0.60	2.86	0.69
60	$- 0.30$	0.33	2.34	0.49	$- 0.29$	0.31	2.14	0.45
80	$- 0.14$	0.17	0.99	0.24	$- 0.14$	0.17	0.98	0.24
100^a	1.66	1.66	4.42	0.80	1.63	1.63	4.03	0.74
cc-pVTZ
20	$- 0.11$	0.26	2.37	0.50	$- 0.09$	0.25	2.15	0.47
40	$- 0.24$	0.32	0.92	0.34	$- 0.23$	0.31	0.91	0.33
60	$- 0.29$	0.33	1.34	0.43	$- 0.28$	0.32	1.34	0.42
80	$- 0.09$	0.11	0.33	0.12	$- 0.09$	0.11	0.33	0.12
100^a	0.19	0.26	1.04	0.28	0.15	0.25	0.86	0.27
cc-pVQZ
20	$- 0.21$	0.25	0.71	0.23	$- 0.21$	0.25	0.70	0.22
40	$- 0.06$	0.07	0.33	0.10	$- 0.06$	0.08	0.30	0.09
60	$- 0.05$	0.09	0.52	0.13	$- 0.04$	0.09	0.49	0.12
80	0.00	0.01	0.05	0.02	0.00	0.01	0.05	0.02
100^a	$- 0.01$	0.09	0.63	0.17	$- 0.06$	0.10	0.64	0.18

^a

Relative to experiment.

Table VI.

Comparison of optimized equilibrium bond angles for different correlation-consistent basis sets (Refs. 65 and 66) for multiple FNO truncations for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The percentage indicates what percent of the virtual space of each molecule was active. For truncated basis sets (20%–80%) errors are relative to the untruncated basis set result; for 100%, errors are relative to experiment. Averages were calculated over the set of molecules from Ref. 60. Only valence electrons were correlated. $\bar{Δ}$ is the (signed) mean error, ${\bar{Δ}}_{abs}$ is the mean absolute error, $Δ_{\max}$ is the maximum absolute error, and $Δ_{std}$ is the standard deviation. All numbers are in units of degrees.

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pVDZ
20	0.18	0.34	0.66	0.36	0.20	0.29	0.64	0.30
40	0.64	0.77	1.36	0.61	0.62	0.76	1.36	0.63
60	0.53	0.60	1.00	0.46	0.51	0.58	0.96	0.46
80	0.13	0.17	0.49	0.21	0.12	0.16	0.47	0.21
100^a	$- 1.99$	1.99	4.97	1.53	$- 1.97$	1.97	4.92	1.52
cc-pVTZ
20	0.10	0.36	0.89	0.46	0.09	0.34	0.85	0.45
40	0.01	0.19	0.47	0.27	0.01	0.18	0.47	0.27
60	0.12	0.19	0.41	0.22	0.12	0.19	0.41	0.22
80	0.07	0.10	0.27	0.11	0.07	0.10	0.28	0.11
100^a	$- 0.90$	0.91	4.26	1.31	$- 0.88$	0.89	4.20	1.29
cc-pVQZ
20	$- 0.03$	0.31	0.68	0.40	$- 0.03$	0.30	0.65	0.40
40	0.08	0.18	0.40	0.23	0.08	0.18	0.41	0.23
60	0.04	0.07	0.14	0.08	0.04	0.07	0.13	0.08
80	0.01	0.02	0.09	0.03	0.01	0.02	0.09	0.03
100^a	$- 0.69$	0.69	3.90	1.23	$- 0.68$	0.68	3.83	1.21

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pVDZ
20	0.18	0.34	0.66	0.36	0.20	0.29	0.64	0.30
40	0.64	0.77	1.36	0.61	0.62	0.76	1.36	0.63
60	0.53	0.60	1.00	0.46	0.51	0.58	0.96	0.46
80	0.13	0.17	0.49	0.21	0.12	0.16	0.47	0.21
100^a	$- 1.99$	1.99	4.97	1.53	$- 1.97$	1.97	4.92	1.52
cc-pVTZ
20	0.10	0.36	0.89	0.46	0.09	0.34	0.85	0.45
40	0.01	0.19	0.47	0.27	0.01	0.18	0.47	0.27
60	0.12	0.19	0.41	0.22	0.12	0.19	0.41	0.22
80	0.07	0.10	0.27	0.11	0.07	0.10	0.28	0.11
100^a	$- 0.90$	0.91	4.26	1.31	$- 0.88$	0.89	4.20	1.29
cc-pVQZ
20	$- 0.03$	0.31	0.68	0.40	$- 0.03$	0.30	0.65	0.40
40	0.08	0.18	0.40	0.23	0.08	0.18	0.41	0.23
60	0.04	0.07	0.14	0.08	0.04	0.07	0.13	0.08
80	0.01	0.02	0.09	0.03	0.01	0.02	0.09	0.03
100^a	$- 0.69$	0.69	3.90	1.23	$- 0.68$	0.68	3.83	1.21

^a

Relative to experiment.

Table VII.

Comparison of optimized equilibrium bond angles for different correlation-consistent basis sets (Refs. 65 and 66) for multiple FNO truncations for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The percentage indicates what percent of the virtual space of each molecule was active. For truncated basis sets (20%–80%) errors are relative to the untruncated basis set result; for 100%, errors are relative to experiment. Averages were calculated over the set of molecules from Ref. 60. All electrons were correlated. $\bar{Δ}$ is the (signed) mean error, ${\bar{Δ}}_{abs}$ is the mean absolute error, $Δ_{\max}$ is the maximum absolute error, and $Δ_{std}$ is the standard deviation. All numbers are in units of degrees.

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pCVDZ
20	0.39	0.40	1.06	0.41	0.38	0.38	1.01	0.39
40	0.43	0.61	1.26	0.65	0.42	0.60	1.25	0.65
60	0.37	0.37	0.65	0.18	0.37	0.37	0.66	0.18
80	0.10	0.24	0.99	0.41	0.10	0.24	0.99	0.41
100^a	$- 1.99$	1.99	4.97	1.53	$- 1.98$	1.98	4.93	1.52
cc-pCVTZ
20	0.17	0.29	0.72	0.34	0.16	0.29	0.73	0.34
40	0.38	0.49	1.02	0.46	0.38	0.50	1.02	0.47
60	0.32	0.37	1.28	0.44	0.32	0.37	1.28	0.44
80	$- 0.03$	0.20	0.53	0.26	$- 0.02$	0.20	0.52	0.26
100^a	$- 1.01$	1.01	4.26	1.28	$- 1.00$	1.00	4.20	1.26
cc-pCVQZ
20	0.02	0.28	0.63	0.37	0.01	0.27	0.62	0.36
40	0.09	0.16	0.39	0.21	0.09	0.16	0.39	0.21
60	0.07	0.11	0.28	0.12	0.07	0.10	0.27	0.12
80	$- 0.01$	0.01	0.04	0.02	$- 0.02$	0.03	0.10	0.04
100^a	$- 0.70$	0.70	3.92	1.23	$- 0.69$	0.69	3.84	1.21

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
cc-pCVDZ
20	0.39	0.40	1.06	0.41	0.38	0.38	1.01	0.39
40	0.43	0.61	1.26	0.65	0.42	0.60	1.25	0.65
60	0.37	0.37	0.65	0.18	0.37	0.37	0.66	0.18
80	0.10	0.24	0.99	0.41	0.10	0.24	0.99	0.41
100^a	$- 1.99$	1.99	4.97	1.53	$- 1.98$	1.98	4.93	1.52
cc-pCVTZ
20	0.17	0.29	0.72	0.34	0.16	0.29	0.73	0.34
40	0.38	0.49	1.02	0.46	0.38	0.50	1.02	0.47
60	0.32	0.37	1.28	0.44	0.32	0.37	1.28	0.44
80	$- 0.03$	0.20	0.53	0.26	$- 0.02$	0.20	0.52	0.26
100^a	$- 1.01$	1.01	4.26	1.28	$- 1.00$	1.00	4.20	1.26
cc-pCVQZ
20	0.02	0.28	0.63	0.37	0.01	0.27	0.62	0.36
40	0.09	0.16	0.39	0.21	0.09	0.16	0.39	0.21
60	0.07	0.11	0.28	0.12	0.07	0.10	0.27	0.12
80	$- 0.01$	0.01	0.04	0.02	$- 0.02$	0.03	0.10	0.04
100^a	$- 0.70$	0.70	3.92	1.23	$- 0.69$	0.69	3.84	1.21

^a

Relative to experiment.

One immediate conclusion is that the FNO convergence behavior is identical for both $CCSD (T)$ and $Λ CCSD (T)$ ⁠. Mean absolute errors (probably the best single measure of the results) are almost identical, especially for larger basis sets. The convergence with respect to truncation of the FNO geometries is not monotonic; while there is generally a trend that less truncation leads to better ${\bar{Δ}}_{abs}$ ⁠, there are exceptions. Even more dramatic are the maximum errors, which do not show a clear convergence behavior. These results are not necessarily surprising. Especially for the double- $ζ$ basis sets, the truncated basis sets can become so small that one cannot consider them meaningful points for extrapolation of the convergence behavior. Unlike the convergence of the energy, the convergence of geometric properties will tend to be less clear-cut: Optimized geometries are dependent not just on the energy at a point, but rather the relative energy at a point to the points around it. There is, therefore, a delicate balance to the best choice of basis and method for geometry prediction, leading to more complicated convergence behavior.

Examining the tables of bond angles (Tables VI and VII), it is clear that both methods underestimate bond angles, even with all the different truncations. Importantly, the FNO truncations do not significantly affect the standard deviations of the geometries.

To better understand the convergence behavior of the FNOs, the mean absolute errors from experiment for this data set are plotted in Figs. 2–5. In these plots, each point represents retaining 20%, 40%, 60%, 80%, or 100% of the virtual space of the corresponding basis. On the horizontal axis is a measure of the relative size of the truncated basis set as compared to the largest basis in the calculation (100% quadruple- $ζ$ ⁠). For bond lengths, all choices of basis sets 60% triple- $ζ$ or larger perform similarly. The picture in the bond angle plots is more mixed, with full convergence not achieved until 40% of the quadruple- $ζ$ basis, though the 60% triple- $ζ$ basis performs quite well. Double- $ζ$ basis sets are inadequate at every truncation. These plots provide a guide for the choice of an optimal basis set of a given size. For example, 20% of a cc-pVQZ basis or 40% of a cc-pVTZ basis yield results that are approximately the same for bond lengths (as shown in Fig. 2) and have the same cost (at the correlated level) as the inferior untruncated cc-pVDZ basis set.

FIG. 2.

View large Download slide

Mean deviation from experiment for bond lengths (in pm) for the equilibrium geometries of the set of molecules from Ref. 60 as a function of correlation-consistent valence basis set and FNO truncation for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The horizontal axis is the average number of virtual basis functions as a percentage of the virtual space of the largest basis cc-pVQZ. Only valence orbitals were correlated.

FIG. 3.

View large Download slide

Mean deviation from experiment for bond lengths (in pm) for the equilibrium geometries of the set of molecules from Ref. 60 as a function of correlation-consistent core-valence basis set and FNO truncation for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The horizontal axis is the average number of virtual basis functions as a percentage of the virtual space of the largest basis cc-pCVQZ. All electrons were correlated.

FIG. 4.

View large Download slide

Mean deviation from experiment for bond angles (in degrees) for the equilibrium geometries of the set of molecules from Ref. 60 as a function of correlation-consistent valence basis set and FNO truncation for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The horizontal axis is the average number of virtual basis functions as a percentage of the virtual space of the largest basis cc-pVQZ. Only valence orbitals were correlated.

FIG. 5.

View large Download slide

Mean deviation from experiment for bond angles (in degrees) for the equilibrium geometries of the set of molecules from Ref. 60 as a function of correlation-consistent core-valence basis set and FNO truncation for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The horizontal axis is the average number of virtual basis functions as a percentage of the virtual space of the largest basis cc-pCVQZ. All electrons were correlated.

The tables provide the required information to rezero the results to those in the complete basis, to isolate the FNO effect from the experimental values. Such plots would reach zero deviation much more rapidly as the tables show.

Even more sensitive to the electron structure method than geometries are vibrational frequencies. In Tables VIII and IX the data for vibrational frequencies as compared to the untruncated basis set results for aug-cc-pVDZ and aug-cc-pVTZ basis sets are shown. For the closed-shell molecules used to calculate the averages in Table VIII, the mean absolute errors are acceptable for basis set truncations of 40% or more, with mean errors of $10 {cm}^{- 1}$ or less. This stands in stark contrast to the open-shell results in Table IX, where 80% of the given basis sets are required to reproduce the untruncated results. For the open-shell molecules, we use UHF reference functions because we have not yet implemented FNOs for restricted open-shell Hartree-Fock (ROHF) reference functions. For cyanide radical, it is known that a ROHF reference function provides significantly better results than UHF (Ref. 61) for perturbation theory, which may be skewing the averages. However, even the $N H_{2}$ results show more dependence than the closed-shell molecules. It is possible that the UHF reference function is the source of this discrepancy, but further work is necessary to verify that conjecture.

Table VIII.

Comparison of vibrational frequencies for the selected closed-shell molecules $H_{2} O$ (Ref. 67), $N H_{3}$ (Ref. 68), $H_{2} C O$ (Refs. 69 and 70), and $C_{2} H_{4}$ (Ref. 71) at equilibrium with different augmented correlation-consistent basis sets (Refs. 65, 66, and 72) for multiple FNO truncations for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The percentage indicates what percent of the virtual space of each molecule was active. For truncated basis sets (20%–80%) errors are relative to the untruncated basis set result; for 100%, errors are relative to experiment. Only valence electrons were correlated. $\bar{Δ}$ is the (signed) mean error, ${\bar{Δ}}_{abs}$ is the mean absolute error, $Δ_{\max}$ is the maximum absolute error, and $Δ_{std}$ is the standard deviation. All numbers are in units of ${cm}^{- 1}$ ⁠.

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
aug-cc-pVDZ
20	$- 11$	27	86	37	$- 10$	27	86	37
40	$- 8$	10	37	10	$- 8$	10	36	10
60	$- 6$	7	34	8	$- 5$	7	34	8
80	$- 1$	1	4	1	$- 1$	1	4	1
100^a	18	26	70	27	17	25	70	27
aug-cc-pVTZ
20	$- 11$	18	54	22	$- 10$	18	53	22
40	$- 6$	8	20	8	$- 6$	8	20	7
60	$- 4$	9	48	15	$- 4$	9	48	15
80	$\approx 0$	1	5	1	$\approx 0$	1	5	1
100^a	4	18	49	25	3	17	48	25

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
aug-cc-pVDZ
20	$- 11$	27	86	37	$- 10$	27	86	37
40	$- 8$	10	37	10	$- 8$	10	36	10
60	$- 6$	7	34	8	$- 5$	7	34	8
80	$- 1$	1	4	1	$- 1$	1	4	1
100^a	18	26	70	27	17	25	70	27
aug-cc-pVTZ
20	$- 11$	18	54	22	$- 10$	18	53	22
40	$- 6$	8	20	8	$- 6$	8	20	7
60	$- 4$	9	48	15	$- 4$	9	48	15
80	$\approx 0$	1	5	1	$\approx 0$	1	5	1
100^a	4	18	49	25	3	17	48	25

^a

Relative to experiment.

Table IX.

Comparison of vibrational frequencies for the selected open-shell radicals CN (Ref. 73) and $N H_{2}$ (Ref. 74) at equilibrium with different augmented correlation-consistent basis sets (Refs. 65, 66, and 72) for multiple FNO truncations for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The percentage indicates what percent of the virtual space of each molecule was active. For truncated basis sets (20%–80%) errors are relative to the untruncated basis set result; for 100%, errors are relative to experiment. Only valence electrons were correlated. $\bar{Δ}$ is the (signed) mean error, ${\bar{Δ}}_{abs}$ is the mean absolute error, $Δ_{\max}$ is the maximum absolute error, and $Δ_{std}$ is the standard deviation. All numbers are in units of ${cm}^{- 1}$ ⁠.

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
aug-cc-pVDZ
20	$- 99$	99	280	122	$- 105$	105	309	137
40	$- 30$	30	101	48	$- 27$	27	91	43
60	$- 16$	16	31	10	$- 17$	17	34	12
80	$- 3$	3	8	3	$- 1$	1	2	1
100^a	$- 69$	69	121	52	$- 76$	76	122	44
aug-cc-pVTZ
20	$- 36$	47	126	63	$- 22$	34	76	41
40	$- 17$	17	45	20	$- 19$	20	58	27
60	$- 12$	12	30	12	$- 13$	13	33	15
80	$\approx 0$	3	6	4	$- 8$	8	30	15
100^a	$- 105$	105	171	52	$- 114$	114	173	51

Basis set(%)	$CCSD (T)$				$Λ CCSD (T)$
Basis set(%)	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$	$\bar{Δ}$	${\bar{Δ}}_{abs}$	$Δ_{\max}$	$Δ_{std}$
aug-cc-pVDZ
20	$- 99$	99	280	122	$- 105$	105	309	137
40	$- 30$	30	101	48	$- 27$	27	91	43
60	$- 16$	16	31	10	$- 17$	17	34	12
80	$- 3$	3	8	3	$- 1$	1	2	1
100^a	$- 69$	69	121	52	$- 76$	76	122	44
aug-cc-pVTZ
20	$- 36$	47	126	63	$- 22$	34	76	41
40	$- 17$	17	45	20	$- 19$	20	58	27
60	$- 12$	12	30	12	$- 13$	13	33	15
80	$\approx 0$	3	6	4	$- 8$	8	30	15
100^a	$- 105$	105	171	52	$- 114$	114	173	51

^a

Relative to experiment.

Deviations from experiment for these sets of molecules are shown in Figs. 6 and 7. More so than the geometries, the deviations in the vibrational frequencies are nonuniform, with different percentages exhibiting radically different agreements with experiment. In the open-shell set, what is immediately clear is that the results agree much more poorly (at all basis sets sizes) with experiment than the closed-shell set. One surprising feature of the open-shell figure is that the augmented double- $ζ$ basis set results are significantly better than the triple- $ζ$ results. This behavior holds for all FNO truncations maintaining more than 40% of the basis set.

FIG. 6.

View large Download slide

Mean deviation from experiment for vibrational frequencies (in ${cm}^{- 1}$ ⁠) for the equilibrium geometries of the closed-shell molecules $H_{2} O$ ⁠, $N H_{3}$ ⁠, $H_{2} C O$ ⁠, and $C_{2} H_{4}$ as a function of augmented correlation-consistent valence basis set and FNO truncation for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The horizontal axis is the average number of virtual basis functions as a percentage of the virtual space of the largest basis aug-cc-pVTZ. Only valence orbitals are correlated.

FIG. 7.

View large Download slide

Mean deviation from experiment for vibrational frequencies (in ${cm}^{- 1}$ ⁠) for the equilibrium geometries of the open-shell molecules CN and $N H_{2}$ as a function of augmented correlation-consistent valence basis set and FNO truncation for $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The horizontal axis is the average number of virtual basis functions as a percentage of the virtual space of the largest basis aug-cc-pVTZ. Only valence orbitals are correlated.

These results are more difficult to interpret than those for geometries and energies. Especially when compared to experiment, the results are much less uniform and show more dependence on the degree of FNO truncation than other properties. This fact should not be surprising; a hessian depends more strongly on the energy differences around the equilibrium structure than does a (first) derivative. A note of caution: It is possible for the FNO procedure to show discontinuities in the vibrational frequencies. We did not see this appear in the results for the set of molecules used here, but in other cases small changes in the truncation level can lead to larger changes in the vibrational frequencies. This dependence illustrates the problem of local smoothness around any given point once the FNO procedure has been applied. We are currently looking more closely at these issues in an attempt to provide less truncation dependent vibrational frequencies for all molecules. There is also the issue of vibrational frequencies from atomic natural orbital basis sets^62,63 compared to those from correlation-consistent ones, which will be considered in future work.

B. Nitroethane

The decomposition of nitroethane can occur via several different pathways.⁴² Schematics of the possible reaction paths are shown in Figs. 8–10. The numbering of transition states and intermediates corresponds to that used in Ref. 42. To sort out the relative importance of each of the individual pathways, all of the relevant species are optimized using FNO $CCSD (T)$ and FNO $Λ CCSD (T)$ ⁠. In the main pathways, there is one reactant (nitroethane), five intermediates, ten transition states, and a total of twelve products. Each of these 28 critical points are fully optimized in a cc-pVTZ basis set with 60% of the virtual space kept using FNOs. This basis sets both performed well in the calibration tests and are small enough to allow the calculations to be completed using our computational resources. Calculations were performed both locally, on our SG1 Altix, as well as at Department of Defense Major Shared Resource Centers. The core occupied orbitals and corresponding core virtual orbitals are dropped as well. For nitroethane and its isomers, this yielded a total of 15 active occupied orbitals and 117 active virtual orbitals. The expected savings per geometry optimization step of each critical point, as compared to a full basis set calculation, is approximately 75%. RHF references are used for closed-shell species, and for open-shell species UHF references are used. At the optimized critical points, finite-difference Hessians are calculated to verify that the geometries did, in fact, correspond to either minima or first-order transition states, as well as to determine the vibrational frequencies, allowing zero-point energy (ZPE) corrections to be included.

FIG. 8.

View large Download slide

Schematic of the one-step HONO elimination and direct fission pathways for decomposition of nitroethane. The vertical axis measures the ZPE corrected energies (in kcal/mol) relative to nitroethane from $Λ CCSD (T)$ calculations with 60% of the virtual space of the cc-pVTZ basis set retained via FNOs.

FIG. 9.

View large Download slide

Schematic of the decomposition pathway of nitroethane through isomerization to ethylnitrite. The vertical axis measures the ZPE corrected energies (in kcal/mol) relative to nitroethane from $Λ CCSD (T)$ calculations with 60% of the virtual space of the cc-pVTZ basis set retained via FNOs.

FIG. 10.

View large Download slide

Schematic of the decomposition pathway of nitroethane through isomerization to ethyl hydroxy nitroxide. The vertical axis measures the ZPE corrected energies (in kcal/mol) relative to nitroethane from $Λ CCSD (T)$ calculations with 60% of the virtual space of the cc-pVTZ basis set retained via FNOs.

The decomposition of nitroethane can be broken into four main classes of pathways: Direct fission of nitroethane to form ethyl radical and nitrogen dioxide (Fig. 8), single-step elimination of HONO (Fig. 8), isomerization to ethylnitrite (denoted INT3) and subsequent decomposition (Fig. 9), and isomerization to ethyl hydroxy nitroxide [CH3CHN(OH)O] (denoted INT5) and then further decomposition (Fig. 10). In the figures mentioned, we have used the notation from Ref. 42 for the intermediates, transition states, and some products (P4 and P8 are two cyclic isomers of nitroethane). When compared to nitromethane, analogies of each of these pathways exists—except for the HONO elimination. For the set of pathways beginning with isomerization to ethylnitrite, we focus on the mechanism that yields the lowest energy products $C H_{3} C H O + H N O$ ⁠. For isomerization through ethyl hydroxy nitroxide, we choose to focus on the thermodynamically minimum set of products $C H_{3} C N O + H_{2} O$ ⁠, elimination of water. To provide an estimate of the importance of these different paths, in Fig. 11 we plot a qualitative picture of their relative energies.

FIG. 11.

View large Download slide

Schematic of the most important pathways for each possible isomerization for the decomposition of nitroethane. The vertical axis measures the ZPE corrected energies (in kcal/mol) relative to nitroethane from $Λ CCSD (T)$ calculations with 60% of the virtual space of the cc-pVTZ basis set retained via FNOs.

Table X compares the results from B3LYP in a $6 - 311 + G (3 d f, 2 p)$ basis and the 60% FNO calculations with $CCSD (T)$ and $Λ CCSD (T)$ at their respective optimized geometries. Focusing first on the B3LYP results from Ref. 42, the energy differences between the different pathways are relatively small. To appropriately model the kinetics of the decomposition of these reactions, it is important that the stationary point energies are converged with respect to electronic structure—small changes in barrier heights can lead to large differences in kinetics.

Table X.

Relative energies of important stationary points for the decomposition of nitroethane in kcal/mol. All species are at their appropriately optimized structures and energies are relative to that of nitroethane including zero-point energy corrections. The B3LYP DFT results are from Ref. 42 and use the $6 - 311 + G (3 d f, 2 p)$ basis. Results for $CCSD (T)$ and $Λ CCSD (T)$ are from this work using a cc-pVTZ basis set with 60% of the virtual orbitals kept by the FNO procedure. Only valence electrons were correlated. Species labels correspond to those in Figs. 8–10.

Species	B3LYP	$CCSD (T)$	$Λ CCSD (T)$
Transition states
TS5	42.11	48.29	48.32
TS6	59.40	64.83	64.94
TS8	35.35	38.07	38.46
TS9	52.50	57.62	57.72
TS10	63.08	60.15	60.81
TS11	64.74	67.27	67.42
TS12	55.93	67.10	68.32
TS13	61.41	60.48	63.30
TS14	31.83	31.07	31.31
TS15	70.41	66.20	68.26
Intermediates
INT3	1.60	$- 0.12$	$- 0.27$
INT5	9.64	14.66	14.61
INT7	$- 5.54$	$- 10.32$	$- 10.54$
Products
$C H_{3} C H_{2} + N O_{2}$	52.32	57.12	56.95
$C H_{2} C H_{2} + H O N O$	15.62	18.35	18.26
$C H_{3} C H O + H N O$	8.62	5.26	4.98
$C H_{3} C N O + H_{2} O$	$- 7.72$	1.57	1.34
$C H_{3} C H_{2} O + N O$	36.22	34.53	34.00
P4	54.40	53.88	54.61
P8	24.47	21.99	21.90

Species	B3LYP	$CCSD (T)$	$Λ CCSD (T)$
Transition states
TS5	42.11	48.29	48.32
TS6	59.40	64.83	64.94
TS8	35.35	38.07	38.46
TS9	52.50	57.62	57.72
TS10	63.08	60.15	60.81
TS11	64.74	67.27	67.42
TS12	55.93	67.10	68.32
TS13	61.41	60.48	63.30
TS14	31.83	31.07	31.31
TS15	70.41	66.20	68.26
Intermediates
INT3	1.60	$- 0.12$	$- 0.27$
INT5	9.64	14.66	14.61
INT7	$- 5.54$	$- 10.32$	$- 10.54$
Products
$C H_{3} C H_{2} + N O_{2}$	52.32	57.12	56.95
$C H_{2} C H_{2} + H O N O$	15.62	18.35	18.26
$C H_{3} C H O + H N O$	8.62	5.26	4.98
$C H_{3} C N O + H_{2} O$	$- 7.72$	1.57	1.34
$C H_{3} C H_{2} O + N O$	36.22	34.53	34.00
P4	54.40	53.88	54.61
P8	24.47	21.99	21.90

Before considering the differences between the coupled-cluster results and those from DFT, note that the results for $CCSD (T)$ and $Λ CCSD (T)$ agree closely, with minimal changes in energy ordering to the different species, despite the fact that $Λ CCSD (T)$ does much better for RHF-based CC bond breaking. Because of this similarity, we will simply refer to the CC results when comparing against B3LYP rather than choosing one or another. Qualitatively, the results from CC and DFT seem to agree quite well; products and intermediates are ordered the same in CC and DFT, and transition states are not radically rearranged. As is noted in Ref. 42, B3LYP tends to underestimate energy barriers; our coupled-cluster results support this conclusion, as the majority of the transition states were determined to be higher in energy than predicted by DFT. The shifts are not uniform, however, leading to a reordering of several of the high-lying transition states.

The lower-lying transition states were left unchanged in order, leading to the same conclusions about the kinetically favored channel. The transition state for the elimination of HONO via a concerted reaction has the lowest barrier by $10 kcal ∕ mol$ in B3LYP and by $8 kcal ∕ mol$ for both $CCSD (T)$ and $Λ CCSD (T)$ ⁠. The concerted nature of this transition state might raise concern about the applicability of the perturbative $CCSD (T)$ method, which fails for RHF-based bond breaking, but recent work³⁹ shows, surprisingly, that $CCSD (T)$ and $Λ CCSD (T)$ (which ameliorates the RHF failure) tend to reproduce transition states with equal accuracy.

From the B3LYP calculations, the elimination of water is the most thermodynamically stable product by more than $16 kcal ∕ mol$ ⁠. On the other hand, the coupled-cluster calculations predict an energy gap between the elimination of water and the elimination of HNO of only $6.8 kcal ∕ mol$ $[CCSD (T)]$ or $6.3 kcal ∕ mol$ $[Λ CCSD (T)]$ ⁠. The elimination of water is exoenergetic in B3LYP by more than $7.5 kcal ∕ mol$ ⁠, while it is endoenergetic by $1.5 kcal ∕ mol$ by both CC methods. When comparing to the energies of the intermediates, the global minimum on the CC potential energy surface is now 1,1-nitrosoethanol (INT7) and ethylnitrite (INT3) is slightly lower in energy than nitroethane. The coupled-cluster calculations also suggest that the elimination of HNO is less favorable kinetically, as the barriers along the reaction pathway are higher relative to those from B3LYP.

V. CONCLUSION

The application of methods that reduce basis set size will always be limited unless analytical gradients are available. For methods such as FNO-CC, where the basis set reduction is based on an auxiliary calculation for the molecule at a particular geometry, the inclusion of orbital relaxation terms is substantially more complicated than it is for more traditional methods that simply modify the orbital eigenvalue equations. In our case, because of the dependence on a MBPT(2) density matrix, there is an orbital relaxation contribution to the two-particle density matrix that is new. Because of the one- and two-particle natures of all the interactions in the Hamiltonian, the most general such truncation procedure should only contribute orbital relaxation effects to both density matrices.

Despite the complexity of the orbital relaxation terms, we are able to show that just as in the case for Hartree-Fock orbitals one can separate the perturbation-dependent integral derivatives from the perturbation-independent orbital relaxation. Therefore, one needs to solve the CPFNO equations (or equivalently, the $z$ -vector equations) once instead of for each perturbation. Then the CC results follow with substantial savings in time that can approach an order of magnitude, depending upon the level of CC correlation. Unfortunately, the price paid for this computational saving is the need to store several quantities of the dimension of two-electron integrals. Proper combination of the terms in an integral-direct formalism may be able to circumvent that complication.

The FNO procedure initiates the optimized virtual space (OVOS) method,^16,17,21,22 which imposes the additional constraint of trying to obtain the lowest MBPT(2) energy¹⁶ or maximizes the overlap between the truncated and untruncated MBPT(2) wavefunctions.^22,64 This constraint can easily be added to the analytical FNO gradient procedure presented here to enable OVOS structures and hessians to be obtained analytically. In fact, the OVOS method is an example of the general issue of imposing additional conditions on a virtual space to fulfill a desired objective.

The application of the FNO truncation methods to the test set of molecules showed that while a cc-pVDZ or cc-pCVDZ basis is inadequate to be predictive for geometries a truncated cc-pVTZ basis of the same number of active orbitals is substantially better. It is always preferable to use the largest possible basis set and then reduce its effective virtual orbital space dimension via the FNO method than to compromise on the size of the underlying basis set. Results for vibrational frequencies are more mixed, without the clear preference for FNO truncations over untruncated smaller basis sets. This conclusion may partly be due to the limited set of molecules studied, but it also suggests that some further developments of the proper treatment of vibrational frequencies within FNO-CC may be needed.

Our results support the general conclusions reached by Denis et al.⁴² about the decomposition of nitroethane. The one-step elimination of HONO appears to be favored kinetically, with the barrier for that reaction which is $8 kcal ∕ mol$ lower than that for the direct bond fission. However, there are important differences in the energetics, with the energies along the pathways initiated by the isomerization to ethylnitrite being most affected. The gap between the thermodynamically favored products $(C H_{3} C N O + H_{2} O)$ and $C H_{3} C H O + H N O$ is reduced to roughly $6 kcal ∕ mol$ versus $16 kcal ∕ mol$ from B3LYP calculations.

ACKNOWLEDGMENTS

The authors would like to recognize support for this work from the Army Research Office through a MURI grant. Computational facilities for our group were supported by the Air Force Office of Scientific Research through DURIP funding. One of the authors (A.G.T.) would also like to acknowledge support from a Department of Defense graduate fellowship.

REFERENCES

1.

K.

Raghavachari

,

G. W.

Trucks

,

J. A.

Pople

, and

M.

Head–Gordon

,

Chem. Phys. Lett.

https://doi.org/10.1016/S0009-2614(89)87395-6

157

,

479

(

1989

).

Google Scholar

Crossref

2.

J. D.

Watts

,

J.

Gauss

, and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.464480

98

,

8718

(

1993

).

Google Scholar

Crossref

3.

R. J.

Bartlett

and

M.

Musial

,

Rev. Mod. Phys.

https://doi.org/10.1103/RevModPhys.79.291

79

,

291

(

2007

).

Google Scholar

Crossref

4.

W. J.

Hunt

and

W. A.

Goddard

III,

Chem. Phys. Lett.

https://doi.org/10.1016/S0009-2614(99)00340-1

3

,

414

(

1969

).

Google Scholar

Crossref

5.

E. R.

Davidson

,

J. Chem. Phys.

https://doi.org/10.1063/1.1678521

57

,

1999

(

1972

).

Google Scholar

Crossref

6.

D.

Feller

and

E. R.

Davidson

,

J. Chem. Phys.

https://doi.org/10.1063/1.441577

74

,

3977

(

1981

).

Google Scholar

Crossref

7.

I. L.

Cooper

and

C. N. M.

Pounder

,

J. Chem. Phys.

https://doi.org/10.1063/1.443677

77

,

5045

(

1982

).

Google Scholar

Crossref

8.

S.

Bird

and

T. A.

Claxton

,

J. Chem. Phys.

https://doi.org/10.1063/1.447151

80

,

3742

(

1984

).

Google Scholar

Crossref

9.

P.

Fantucci

,

V.

Bonačić–Koutechký

, and

J.

Koutechký

,

J. Comput. Chem.

https://doi.org/10.1002/jcc.540060515

6

,

462

(

1985

).

Google Scholar

Crossref

10.

W. L.

Luken

and

B. A. B.

Seiders

,

Chem. Phys.

https://doi.org/10.1016/0301-0104(85)85019-9

92

,

235

(

1985

).

Google Scholar

Crossref

11.

F.

Illas

,

M.

Merchan

,

M.

Pelissier

, and

J. P.

Malrieu

,

Chem. Phys.

https://doi.org/10.1016/0301-0104(86)85014-5

107

,

361

(

1986

).

Google Scholar

Crossref

12.

H. P.

Kelly

,

Perturbation Theory and its Application in Quantum Mechanics

(

Wiley

,

New York

,

1966

), pp.

215

–

241

.

Google Scholar

13.

D. M.

Silver

and

R. J.

Bartlett

,

Phys. Rev. A

https://doi.org/10.1103/PhysRevA.13.1

13

,

1

(

1976

).

Google Scholar

Crossref

14.

D. M.

Silver

,

S.

Wilson

, and

R. J.

Bartlett

,

Phys. Rev. A

https://doi.org/10.1103/PhysRevA.16.477

16

,

477

(

1977

).

Google Scholar

Crossref

15.

P. R.

Taylor

,

G. B.

Bacskay

,

N. S.

Hush

, and

A. C.

Hurley

,

Chem. Phys. Lett.

https://doi.org/10.1016/0009-2614(76)85390-0

41

,

444

(

1976

).

Google Scholar

Crossref

16.

L.

Adamowicz

and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.452468

86

,

6314

(

1987

).

Google Scholar

Crossref

17.

L.

Adamowicz

,

R. J.

Bartlett

, and

A.

Sadlej

,

J. Chem. Phys.

https://doi.org/10.1063/1.454721

88

,

5749

(

1988

).

Google Scholar

Crossref

18.

C.

Sosa

,

J.

Geertsen

,

G. W.

Trucks

,

R. J.

Bartlett

, and

J. A.

Franz

,

Chem. Phys. Lett.

https://doi.org/10.1016/0009-2614(89)87399-3

159

,

148

(

1989

).

Google Scholar

Crossref

19.

K. J.

Wilson

, Ph.D. thesis,

University of Florida

,

2002

.

20.

A. G.

Taube

and

R. J.

Bartlett

,

Collect. Czech. Chem. Commun.

https://doi.org/10.1135/cccc20050837

70

,

837

(

2005

).

Google Scholar

Crossref

21.

M.

Pitoňák

,

P.

Neogrády

,

V.

Kellö

, and

M.

Urban

,

Mol. Phys.

104

,

2277

(

2006

).

Google Scholar

22.

M.

Pitoňák

,

F.

Holka

,

P.

Neogrády

, and

M.

Urban

,

J. Mol. Struct.: THEOCHEM

https://doi.org/10.1016/j.theochem.2006.05.018

768

,

79

(

2006

).

Google Scholar

Crossref

23.

R.

Ahlrichs

and

W.

Kutzelnigg

,

Theor. Chim. Acta

https://doi.org/10.1007/BF00528769

10

,

377

(

1968

).

Google Scholar

Crossref

24.

R.

Ahlrichs

,

H.

Lischka

,

B.

Zurawski

, and

W.

Kutzelnigg

,

J. Chem. Phys.

https://doi.org/10.1063/1.431254

63

,

4685

(

1975

).

Google Scholar

Crossref

25.

C.

Edmiston

and

M.

Krauss

,

J. Chem. Phys.

https://doi.org/10.1063/1.1727841

45

,

1833

(

1966

).

Google Scholar

Crossref

26.

M.

Gelus

,

R.

Ahlrichs

,

V.

Staemmler

, and

W.

Kutzelnigg

,

Chem. Phys. Lett.

https://doi.org/10.1016/0009-2614(70)80159-2

7

,

503

(

1970

).

Google Scholar

Crossref

27.

T.

Barr

and

E. R.

Davidson

,

Phys. Rev. A

https://doi.org/10.1103/PhysRevA.1.644

1

,

644

(

1970

).

Google Scholar

Crossref

28.

H.

Jørgen

,

A.

Jensen

,

P.

Jørgensen

,

H.

Ågren

, and

J.

Olsen

,

J. Chem. Phys.

https://doi.org/10.1063/1.453884

88

,

3834

(

1988

).

Google Scholar

Crossref

29.

P. J.

Hay

,

J. Chem. Phys.

https://doi.org/10.1063/1.1680359

59

,

2468

(

1973

).

Google Scholar

Crossref

30.

W.

Meyer

,

J. Chem. Phys.

https://doi.org/10.1063/1.1679283

58

,

1017

(

1973

).

Google Scholar

Crossref

31.

B.

Mintz

,

K.

Lennox

, and

A.

Wilson

,

J. Chem. Phys.

https://doi.org/10.1063/1.1785145

121

,

5629

(

2004

).

Google Scholar

Crossref

PubMed

32.

B.

Mintz

and

A.

Wilson

,

J. Chem. Phys.

https://doi.org/10.1063/1.1867412

122

,

134106

(

2005

).

Google Scholar

Crossref

PubMed

33.

B.

Prascher

,

B.

Wilson

, and

A.

Wilson

,

J. Chem. Phys.

https://doi.org/10.1063/1.2768602

127

,

124110

(

2007

).

Google Scholar

Crossref

PubMed

34.

P.

Politzer

and

J. S.

Murray

,

Energetic Materials, Part 1. Decomposition, Crystal and Molecular Properties

,

Theoretical and Computational Chemistry

, Vol.

12

(

Elsevier

,

Amsterdam, the Netherlands

,

2003

).

Google Scholar

Crossref

35.

T. B.

Brill

,

T. P.

Russell

,

W. C.

Tao

, and

R. B.

Wardle

, in

Decomposition, Combustion and Detonation Chemistry of Energetic Materials

(

Materials Research Society

,

Pittsburgh, PA

,

1996

), Vol.

418

.

Google Scholar

36.

M. L. J.

McKee

,

J. Am. Chem. Soc.

https://doi.org/10.1021/ja00279a020

108

,

5784

(

1986

).

Google Scholar

Crossref

PubMed

37.

M. T.

Nguyen

,

H. T.

Le

,

B.

Hajgató

,

T.

Veszprémi

, and

M. C.

Lin

,

J. Phys. Chem. A

https://doi.org/10.1021/jp027532h

107

,

4286

(

2003

).

Google Scholar

Crossref

38.

W.-F.

Hu

,

T.-J.

He

,

D.-M.

Chen

, and

F.-C.

Liu

,

J. Phys. Chem. A

https://doi.org/10.1021/jp020070n

106

,

7294

(

2002

).

Google Scholar

Crossref

39.

A. G.

Taube

and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.2830237

128

,

044111

(

2008

).

Google Scholar

Crossref

PubMed

40.

P.

Gray

,

A. D.

Yoffe

, and

L.

Roselaar

,

Trans. Faraday Soc.

https://doi.org/10.1039/tf9555101489

51

,

1489

(

1955

).

Google Scholar

Crossref

41.

G.

Spokes

and

S.

Benson

,

J. Am. Chem. Soc.

https://doi.org/10.1021/ja01000a007

89

,

6030

(

1967

).

Google Scholar

Crossref

42.

P.

Denis

,

O.

Ventura

,

H.

Le

, and

M. T.

Nguyen

,

Phys. Chem. Chem. Phys.

https://doi.org/10.1039/b300275f

5

,

1730

(

2003

).

Google Scholar

Crossref

43.

S. A.

Kucharski

and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.475961

108

,

5243

(

1998

).

Google Scholar

Crossref

44.

T. D.

Crawford

and

J. F.

Stanton

,

Int. J. Quantum Chem.

https://doi.org/10.1002/(SICI)1097-461X(1998)70:4/5<601::AID-QUA6>3.0.CO;2-Z

70

,

601

(

1998

).

Google Scholar

Crossref

45.

A. G.

Taube

and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.2830236

128

,

044110

(

2008

).

Google Scholar

Crossref

PubMed

46.

E. A.

Salter

,

G. W.

Trucks

, and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.456069

90

,

1752

(

1989

).

Google Scholar

Crossref

47.

J.

Gauss

,

J. F.

Stanton

, and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.460916

95

,

2639

(

1991

).

Google Scholar

Crossref

48.

K. K.

Baeck

,

J. D.

Watts

, and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.474744

107

,

3853

(

1997

).

Google Scholar

Crossref

49.

J.

Gerratt

and

I.

Mills

,

J. Chem. Phys.

https://doi.org/10.1063/1.1670299

49

,

1719

(

1968

).

Google Scholar

Crossref

50.

A.

Dalgarno

and

A.

Stewart

,

Proc. R. Soc. London, Ser. A

238

,

269

(

1956

).

Google Scholar

51.

R. J.

Bartlett

,

C. E.

Dykstra

, and

J.

Paldus

, in

Advanced Theories and Computational Approaches to Electronic Structure of Molecules

, edited by

C. E.

Dykstra

(

Reidel

,

Dordrecht

,

1984

), pp.

127

–

159

.

Google Scholar

Crossref

52.

N. C.

Handy

and

H. F.

Schaefer

III,

J. Chem. Phys.

https://doi.org/10.1063/1.447489

81

,

5031

(

1984

).

Google Scholar

Crossref

53.

R. J.

Bartlett

,

Molecular Quantum Mechanics: Selected Papers of N. C. Handy

(

Taylor & Francis

,

London

,

2004

), pp.

127

–

130

.

Google Scholar

54.

R. J.

Harrison

,

G. B.

Fitzgerald

,

W. D.

Laidig

, and

R. J.

Bartlett

,

Chem. Phys. Lett.

https://doi.org/10.1016/0009-2614(86)87048-8

124

,

291

(

1986

).

Google Scholar

Crossref

55.

T. J.

Lee

and

A. P.

Rendell

,

J. Chem. Phys.

https://doi.org/10.1063/1.460411

94

,

6229

(

1991

).

Google Scholar

Crossref

56.

K.

Toyota

,

M.

Ehara

, and

H.

Nakatsuji

,

Chem. Phys. Lett.

https://doi.org/10.1016/S0009-2614(02)00250-6

356

,

1

(

2002

).

Google Scholar

Crossref

57.

J. E.

Subotnik

and

M.

Head–Gordon

,

J. Chem. Phys.

https://doi.org/10.1063/1.2000252

123

,

064108

(

2005

).

Google Scholar

Crossref

58.

J. E.

Subotnik

,

A. D.

Dutoi

, and

M.

Head–Gordon

,

J. Chem. Phys.

https://doi.org/10.1063/1.2033687

123

,

114108

(

2005

).

Google Scholar

Crossref

PubMed

59.

J. F.

Stanton

,

J.

Gauss

,

J. D.

Watts

,

M.

Nooijen

,

N.

Oliphant

,

S. A.

Perera

,

P. G.

Szalay

,

W. J.

Lauderdale

,

S. A.

Kucharski

,

S. R.

Gwaltney

,

S.

Beck

,

A.

Balková

,

D. E.

Bernholdt

,

K. K.

Baeck

,

P.

Rozyczko

,

H.

Sekino

,

C.

Hober

, and

R. J.

Bartlett

, ACES II, a program product of the Quantum Theory Project (University of Florida, Gainesville, FL). Integral packages included are VMOL (J. Almlöf and P. R. Taylor), VPROPS (P. Taylor), ABACUS (T. Helgaker, H. J. Aa. Jensen, P. Jørgensen, J. Olsen, and P. R. Taylor).

60.

K. L.

Bak

,

J.

Gauss

,

P.

Jørgensen

,

J.

Olsen

,

T.

Helgaker

, and

J. F.

Stanton

,

J. Chem. Phys.

https://doi.org/10.1063/1.1357225

114

,

6548

(

2001

).

Google Scholar

Crossref

61.

W. J.

Lauderdale

,

J. F.

Stanton

,

J.

Gauss

,

J. D.

Watts

, and

R. J.

Bartlett

,

J. Chem. Phys.

https://doi.org/10.1063/1.463664

97

,

6606

(

1992

).

Google Scholar

Crossref

62.

P.-O.

Widmark

,

P. A.

Malmqvist

, and

B.

Roos

,

Theor. Chim. Acta

https://doi.org/10.1007/BF01120130

77

,

291

(

1990

).

Google Scholar

Crossref

63.

P.-O.

Widmark

,

B. J.

Persson

, and

B.

Roos

,

Theor. Chim. Acta

https://doi.org/10.1007/BF01112569

79

,

419

(

1991

).

Google Scholar

Crossref

64.

P.

Neogrády

,

M.

Pitoňák

, and

M.

Urban

,

Mol. Phys.

https://doi.org/10.1080/00268970500096251

103

,

2141

(

2005

).

Google Scholar

Crossref

65.

T. H.

Dunning

, Jr.,

J. Chem. Phys.

https://doi.org/10.1063/1.456153

90

,

1007

(

1989

).

Google Scholar

Crossref

66.

D. E.

Woo

and

T. H.

Dunning

, Jr.,

J. Chem. Phys.

https://doi.org/10.1063/1.470645

103

,

4572

(

1995

).

Google Scholar

Crossref

67.

A. R.

Hoy

,

I. M.

Mills

, and

G.

Strey

,

Mol. Phys.

https://doi.org/10.1080/00268977200102361

24

,

1265

(

1972

).

Google Scholar

Crossref

68.

J. L.

Duncan

and

I. M.

Mills

,

Spectrochim. Acta

https://doi.org/10.1016/0371-1951(64)80047-3

20

,

523

(

1964

).

Google Scholar

Crossref

69.

D. E.

Reisner

,

R. W.

Field

,

J. L.

Kinsey

, and

H. L.

Dai

,

J. Chem. Phys.

https://doi.org/10.1063/1.446677

80

,

5968

(

1984

).

Google Scholar

Crossref

70.

M. W.

Wohar

and

P. W.

Jadgodzinski

,

J. Mol. Spectrosc.

https://doi.org/10.1016/0022-2852(91)90030-E

148

,

13

(

1991

).

Google Scholar

Crossref

71.

J. L.

Duncan

,

D. C.

McKean

, and

P. D.

Mallinson

,

J. Mol. Spectrosc.

https://doi.org/10.1016/0022-2852(73)90154-9

45

,

221

(

1973

).

Google Scholar

Crossref

72.

R. A.

Kendall

,

T. H.

Dunning

, Jr., and

R. J.

Harrison

,

J. Chem. Phys.

https://doi.org/10.1063/1.462569

96

,

6796

(

1992

).

Google Scholar

Crossref

73.

K. P.

Huber

and

G.

Herzberg

,

Constants of Diatomic Molecules

,

Molecular Spectra and Molecular Structure

, Vol.

4

(

Van Nostrand and Reinhold

,

New York

,

1979

).

Google Scholar

74.

K.

Dressler

and

D. A.

Ramsay

,

Philos. Trans. R. Soc. London, Ser. A

https://doi.org/10.1098/rsta.1959.0011

251

,

553

(

1959

).

Google Scholar

Crossref

2008

American Institute of Physics

Frozen natural orbital coupled-cluster theory: Forces and application to decomposition of nitroethane

I. INTRODUCTION

II. THEORY

A. Frozen natural orbital CC

B. Gradients

1. Noncanonical perturbed orbitals

2. Canonical gradients

C. Smoothness of the potential energy surface

III. IMPLEMENTATION

IV. RESULTS AND DISCUSSION

A. Calibration

B. Nitroethane

V. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Frozen natural orbital coupled-cluster theory: Forces and application to decomposition of nitroethane Open Access

I. INTRODUCTION

II. THEORY

A. Frozen natural orbital CC

B. Gradients

1. Noncanonical perturbed orbitals

2. Canonical gradients

C. Smoothness of the potential energy surface

III. IMPLEMENTATION

IV. RESULTS AND DISCUSSION

A. Calibration

B. Nitroethane

V. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Frozen natural orbital coupled-cluster theory: Forces and application to decomposition of nitroethane