Communication: Automatic code generation enables nuclear gradient computations for fully internally contracted multireference theory

Substantial effort has been devoted to implementing complex formulas in quantum chemistry, resulting in accurate and useful computational tools for chemical applications. However, some equations are too complicated to be handled manually. Among such examples is the analytical nuclear gradient theory of fully internally contracted complete active space second-order perturbation theories (FIC-CASPT2 or simply CASPT2)^1–3 and its variants,^4–6 whose complexity has hampered their implementation for more than two decades since FIC-CASPT2 was developed.¹ In this study, we address this problem by employing an automatic code generation approach to help enable many chemical applications. For instance, such implementations can be used for the geometry optimization⁷ of strongly correlated molecules as complex as those investigated by respective single-point calculations. They also have the potential to replace mean-field models (e.g., complete active space self-consistent field or CASSCF) often used in photodynamics simulations involving the ground and excited states of molecules.⁸ 

The challenge in implementing the nuclear gradients for the FIC-CASPT2 method can be easily recognized as follows. The state-specific CASPT2 energy functional for the nth state is

in which $f ˆ$ is the standard state-specific Fock operator, and $E 0 ( n ) = 〈 Φ 0 | f ˆ | Φ 0 〉$⁠. The reference and correlated wave functions are defined as (using the reference CI coefficients $c I ( n )$⁠)

where I labels Slater determinants, and $E ˆ Ω$ are external and semi-external excitation operators (such as $E ˆ a i , b j$ and $E ˆ r s , a t$⁠).¹ In this article, a and b label virtual orbitals, i and j label closed orbitals, r, s, t, and u label active orbitals, and x, y, z, and w label any orbitals. The energy functional $E$ is minimized with respect to the T_Ω amplitudes to give the CASPT2 energy. The nuclear energy gradients are the total derivative of the energy with respect to nuclear displacements R,

in which C is the molecular-orbital coefficient matrix parameterized by an anti-Hermitian matrix κ,

The first term on the right hand side of Eq. (4) is the Hellmann–Feynman force,⁹ and the second term is zero in the simplest state-specific CASPT2 case because $E$ is stationary with respect to variation of T_Ω. The third term also appears in the standard single-reference algorithms. The complexity of the equations for FIC-CASPT2 nuclear gradients stems from the last term, which is associated with the reference-coefficient derivatives; since the first-order wave functions are expanded in a basis that is dependent on $c I ( n )$ [Eq. (3)], every single term in $E$ contributes to $∂E/∂ c I ( n )$ in a nontrivial way.

The use of partially internally contracted or uncontracted basis functions¹⁰ (referred to as WK-CASPT2 in the following) for first-order wave functions greatly simplifies the equations, making it tractable to manually implement nuclear gradients for such variants.^11–15 This is because, in these methods, parts (or all) of the first-order wave function are expanded in terms of excited Slater determinants,

which are not dependent on $c I ( n )$ (note that in practice only distinct determinants are included in the sum). There is also an implementation of nuclear gradients of uncontracted multireference configuration interaction (MRCI).¹⁶ However, the formal scaling of the size of first-order wave functions in these methods is factorial with respect to the number of active orbitals, rendering them sub-optimal for large calculations.⁴ 

Over the last decade, various automatic code generation approaches have been developed to replace tedious, error-prone manual implementation processes in quantum chemistry.¹⁷ The automated higher-order coupled-cluster (CC) implementations by Kállay and Surján¹⁸ and by Hirata^19,20 were the first to demonstrate that the automation strategy can produce programs that are competitive with hand-optimized codes in serial and massively parallel environments, respectively. Hanrath and co-workers have realized an arbitrary-order CC code generator that is almost optimal.^21,22 Recently, this strategy has been extended to various methods such as CC-F12,^23–26 relativistic CC,²⁷ local CC,^28,29 open-shell CC,^30,31 and excited-state CC methods.^20,32 Note that these automation techniques are complementary to meta-language approaches (such as sial,³³ itf,³⁴ libtensor,³⁵ tiledarray,³⁶ to name a few), because code generators can be used to synthesize computer codes in any meta-language as well.

Furthermore, the automation strategy has been applied to development of multireference electron correlation theories. Neuscamman et al. used an automated scheme for the canonical transformation theory;³⁷ Parkhill et al. developed local active-space methods (e.g., an active-space perfect quadruples model³⁸) using a sparse framework;³⁹ Hanauer and Köhn extended their string-based code to implement an internally contracted MRCC method.⁴⁰ Saitow et al. reported a fully internally contracted MRCI method based on density-matrix-renormalization-group reference functions.⁴¹ 

In this work, we extend the automatic code generation approach to realize FIC-CASPT2 nuclear gradients that have been sought for a long time. Following the standard approach,^9,13 we use the CASPT2 Lagrangian and the so-called Z-vector equation,⁴² instead of directly evaluating Eq. (4), to avoid computation of geometry derivatives of wave function parameters (such as $d c I ( m ) /dR$⁠). The state-specific CASPT2 Lagrangian is¹³ 

Each term on the right hand side (other than $E$⁠) corresponds to a constraint arising from the CASSCF reference calculation. A is the orbital gradient of CASSCF, and Z is its Lagrange multiplier. S is the overlap matrix, and X is the Lagrange multiplier for orbital orthogonality. The next two terms are related to the stationary condition in the full configuration interaction in the active space performed in CASSCF. Here, m labels states averaged in CASSCF and W_m is the weight in the state averaging. The final term originates from the use of the frozen core approximation in CASPT2. See details in Ref. 13. When the stationary conditions on $L$ with respect to all the parameters and multipliers are met, the nuclear gradients can be computed as

since only molecular integrals have explicit dependence on the nuclear displacement R. We will consider level shifts⁴³ in future work, which requires the additional implementation of the so-called λ equation (as do multistate variants¹⁴).

We have developed a new automated code generator, named smith3, that derives equations on the basis of the spin-free version of Wick’s theorem in which the normal ordering is defined with respect to the closed-Fock vacuum. In addition, smith3 expresses the terms that involve active-orbital indices as a sum of canonical terms so that they can be computed from canonical density matrices and its derivatives, e.g.,

with σ and ρ labeling spins. Note that $( Γ 0 ) Λ = 〈 Φ 0 | E ˆ Λ | Φ 0 〉$ and $( Γ 0 ) Λ I = 〈 I | E ˆ Λ | Φ 0 〉$⁠, where $E ˆ Λ$ is a general operator. Here, smith3 is used to implement the following expressions:

in which $|Ψ〉= ∑ I t I |I〉$ is any multi-configuration reference function, and $G ˆ = ∑ Λ G Λ E ˆ Λ$ and $R ˆ = ∑ Ω R Ω E ˆ Ω$ are general and excitation operators, respectively. Note that the determinant index I is treated analogously to the orbital indices in the generated code; for instance, $( Γ 0 ) t u I$ is viewed as a three-index tensor whose size is $n det n act 2$ (n_det and n_act are the numbers of the determinants and the active orbitals, respectively).

Using this machinery, we automate the nuclear-gradient implementation as follows. First, to optimize T_Ω, the program for computing the CASPT2 residual vectors is generated using Eq. (10a), i.e.,

Second, the Z-CASSCF equation is solved, which is a set of coupled equations defined as

For Eq. (12a), the reference-coefficient derivatives of the CASPT2 energy,

are implemented by smith3 using Eq. (10c). Next, for Eq. (12b) the MO-coefficient derivatives of the CASPT2 energy,

are calculated from the one- and two-body density matrices (see Refs. 13 and 15 for explicit formulas). The density matrices are defined as

and are implemented using Eq. (10b). Given $y I ( n )$ and Y_rs, solutions of Z-CASSCF [Z, X, $z I ( m )$⁠, and z_ij in Eq. (7)] can be obtained as detailed in Ref. 13. Using these wave function parameters and the Lagrange multipliers, effective density matrices are formed, which are then contracted to two-index and three-index gradient integrals.^13,15

The generated code uses a tile-based data layout that is similar to that used in tce ¹⁹ and in the earlier version of smith.^23,24 All the codes implemented in the bagel package⁴⁴ and the code generator smith3 ⁴⁵ are openly available under the GNU General Public License. The technical details on the implementation, working equations, and source code of the smith3 program are also found in supplementary material.⁴⁶ CASSCF and Z-CASSCF were manually implemented in bagel using density fitting (DF) for efficiency as reported in Ref. 15. In CASPT2, four-index two-external integrals were explicitly constructed from DF integrals. The smith3 program generated ca. 1150 tasks, the majority of which are tensor contractions. Each task is expressed as a node of a tree-like directed acyclic graph, which we traverse at runtime. This infrastructure should assist in interfacing smith3 code to parallel-runtime libraries in the future.

First, to show the numerical impact of full internal contraction on geometrical parameters, we optimized the ground-state geometry of the N,N ′-diiminato-copper-dioxygen complex [(H₅C₃N₂)CuO₂] in its side-on coordination configuration. The ground state is singlet. We used the (14e, 9o) active space consisting of an in-plane d orbital of copper and eight valence orbitals of dioxygen as suggested in earlier work.^15,34 The aug-cc-pVDZ^47,48 and def2-QZVPP/JKFIT⁴⁹ basis sets were used for orbital and auxiliary functions, respectively. Table I compiles optimized Cu–O and O–O bond lengths. It is evident that neither the degrees of internal contraction (i.e., CASPT2 and WK-CASPT2) nor the DF approximation has impact on the bond lengths. Our program did not take advantage of spatial symmetry. One optimization step took about 13 min using 2 Xeon E5-2650 CPUs (2.0 GHz). More than half the time is spent for evaluation of Eq. (13). Note, however, that the code generated by smith3 has not been threaded efficiently, and there is room for further improvement.

Next, we calculated the vertical and adiabatic ionization potentials (IPs) of the porphin molecule (C₂₀H₁₄N₄) using the optimized geometries computed by CASPT2. The cc-pVDZ basis set⁴⁷ was used together with the corresponding JKFIT basis set⁵¹ for DF. The (4e, 4o) and (3e, 4o) active spaces were used,⁴⁶ which consist of the four frontier orbitals of Gouterman’s model.⁵² The numbers of (correlated) inactive and virtual orbitals were 55 and 323, respectively. One optimization step took about 30 min on the same hardware. For comparison, we also computed the IPs from the PBE functional⁵³ and MP2 (with an unrestricted variant for the radical cation) using turbomole.⁵⁴ The PBE calculations were performed using the def2-SVP basis set.⁵⁵ Geometry optimization using all methods including CASPT2 was performed without imposing spatial symmetry; the geometry of the neutral porphin was found to belong to the D_2h symmetry group, whereas that of the radical cation was found to be C_2h (even when an initial geometry was set to a D_2h structure, optimization converged to this minimum). Similar symmetry breaking of metalloporphyrin cation radicals due to the pseudo-Jahn–Teller effect has been reported in the literature.⁵⁶ The optimized geometry of the radical cation from unrestricted PBE was found to be D_2h. We were not able to optimize the geometry of the radical cation using unrestricted MP2 due to wave function instability. The geometrical parameters are compiled in supplementary material.⁴⁶ The IPs are shown in Table II. The vertical IP computed by CASPT2 (6.84 eV) is in good agreement with the experimental value (6.9 eV).⁵⁷ The difference between the vertical and adiabatic IPs computed at 0.18 eV is an order of magnitude larger than that computed using PBE. Furthermore, we computed the IPs with CASPT2 using the PBE-optimized geometries. As expected, the vertical IP is almost identical to that computed at the CASPT2 geometry; however, the adiabatic IP is larger than the vertical IP, which attests to the importance of geometry optimization at the CASPT2 level.

In summary, we have used automatic code generation to realize analytical CASPT2 nuclear gradients with full internal contraction. Our implementation has been applied to the N,N ′-diiminato-copper-dioxygen complex to show that errors due to full internal contraction are marginal. We have also computed the vertical and adiabatic IPs of the porphin molecule to illustrate the capability of our implementation. There is, however, room for improvement in our program in terms of efficiency and storage requirement; currently, application of our program is limited by the storage of

and the T_Ω amplitudes of $n occ 2 n bas 2$ size (n_occ and n_bas are the numbers of occupied orbitals and basis functions, respectively). Further optimization and distributed-memory parallelization of our program are warranted and will be performed in the future. We will also consider the level shift,⁴³ other zeroth-order Hamiltonians,^6,58 and multistate extensions.^14,59,60

The debugging of parts of the code in bagel, including Z-CASSCF, was facilitated by existing implementations in molpro.⁵⁰ This work has been supported by Department of Energy, Basic Energy Sciences (Grant No. DE-FG02-13ER16398), and the Air Force Office of Scientific Research Young Investigator Program (Grant No. FA9550-15-1-0031).