Benchmark of dynamic electron correlation models for seniority-zero wavefunctions and their application to thermochemistry

TL;DR

This paper develops and benchmarks new perturbation theory models combined with AP1roG wavefunctions to accurately describe both static and dynamic electron correlation effects, improving thermochemistry predictions.

Contribution

The authors extend perturbation theory models with AP1roG to include dynamic correlation, and benchmark their performance against high-level reference data for reaction energies and bond-breaking processes.

Findings

Linearized coupled cluster correction performs best among tested models.

Second-order perturbation theory with specific Hamiltonian choices yields reliable results.

Models accurately describe systems with both static and dynamic electron correlation.

Abstract

Wavefunctions restricted to electron-pair states are promising models to describe static/nondynamic electron correlation effects encountered, for instance, in bond-dissociation processes and transition-metal and actinide chemistry. To reach spectroscopic accuracy, however, the missing dynamic electron correlation effects that cannot be described by electron-pair states need to be included \textit{a posteriori}. In this article, we extend the previously presented perturbation theory models with an Antisymmetric Product of 1-reference orbital Geminal (AP1roG) reference function that allow us to describe both static/nondynamic and dynamic electron correlation effects. Specifically, our perturbation theory models combine a diagonal and off-diagonal zero-order Hamiltonian, a single-reference and multi-reference dual state, and different excitation operators used to construct the projection manifold. We benchmark all proposed models as well as an \textit{a posteriori} linearized coupled cluster correction on top of AP1roG against CR-CCSD(T) reference data for reaction energies of several closed-shell molecules that are extrapolated to the basis set limit. Moreover, we test the performance of our new methods for multiple bond breaking processes in the N$_2$, C$_2$, and BN dimers against MRCI-SD and MRCI-SD+Q reference data. Our numerical results indicate that the best performance is obtained from a linearized coupled cluster correction as well as second-order perturbation theory corrections employing a diagonal and off-diagonal zero-order Hamiltonian and a single-determinant dual state. These dynamic corrections on top of AP1roG allow us to reliably model molecular systems dominated by static/nondynamic as well as dynamic electron correlation.

Tables

Table 1. Table 1: Summary of PT models with zero-order Hamiltonian H ^ 0 subscript ^ 𝐻 0 \hat{H}_{0} , perturbation V ^ ^ 𝑉 \hat{V} , dual ⟨ Ψ ~ | bra ~ Ψ \langle\tilde{\Psi}| , and excitation operator T ^ ^ 𝑇 \hat{T} . All operators are defined in the text. The computational scaling is given in the last column.

Model	${\hat{H}}_0$	$\hat{V}$	$⟨\tilde{\Psi }|$	$\hat{T}$	scaling
PT2SDd	${\hat{F}}_N^d$	$F_N^o+{\hat{W}}_N^{\prime }$	$⟨0|$	${\hat{T}}_2^{\prime }$, ${\hat{T}}_1+{\hat{T}}_2^{\prime }$	$𝒪\left(o^2{v}^2\right)$
PT2MDd	${\hat{F}}_N^d$	$F_N^o+{\hat{W}}_N^{\prime }$	$⟨\mathrm{AP1roG}|$	${\hat{T}}_2^{\prime }$, ${\hat{T}}_1+{\hat{T}}_2^{\prime }$	$𝒪\left(o^2{v}^2\right)$
PT2SDo	${\hat{F}}_N^d+{\hat{F}}_N^o$	${\hat{W}}_N^{\prime }$	$⟨0|$	${\hat{T}}_2^{\prime }$, ${\hat{T}}_1+{\hat{T}}_2^{\prime }$	$𝒪\left(o^2{v}^3\right)$
PT2MDo	${\hat{F}}_N^d+{\hat{F}}_N^o$	${\hat{F}}_N-{\overline{F}}_N+{\hat{W}}_N^{\prime }$	$⟨\mathrm{AP1roG}|$	${\hat{T}}_2^{\prime }$, ${\hat{T}}_1+{\hat{T}}_2^{\prime }$	$𝒪\left(o^2{v}^3\right)$
PT2b	${\hat{F}}_N^d+{\hat{F}}_N^o$	${\hat{H}}_N^{\prime }$	$⟨\mathrm{AP1roG}|$	${\hat{T}}_2$, ${\hat{T}}_2^{\prime }$, ${\hat{T}}_1+{\hat{T}}_2$, ${\hat{T}}_1+{\hat{T}}_2^{\prime }$	$𝒪\left(o^2{v}^3\right)$

Table 2. Table 2: Spectroscopic constants for the dissociation of the C 2 , N 2 , and BN molecule for different quantum chemistry methods. The differences are with respect to MRCI-SD reference data  Peterson  et al. ( 1993 ) for C 2 and N 2 , and with respect to MRCI-SD+Q reference data  Peterson ( 1995 ) for BN. E e subscript 𝐸 𝑒 E_{e} : ground state energy at r e subscript 𝑟 𝑒 r_{e} .

	Method	$E_e$ [$E_h$]	$r_e$ [Å]	$D_e$ [$\frac{\mathrm{kcal}}{\mathrm{mol}}$]		${\omega }_e$ [cm-1]
C2	AP1roG	$-$75.58569	1.227($-$0.025)	132.9(	$-$8.2)	1780(	$-$56)
	AP1roG-PT2SDd(d)	$-$75.80222	1.251($+$0.001)	160.4(	$+$19.3)	1889(	$+$53)
	AP1roG-PT2SDd(sd)	$-$75.81041	1.249($-$0.003)	154.6(	$+$13.5)	1915(	$+$79)
	AP1roG-PT2MDd(d)	$-$75.77832	1.239($-$0.013)	121.6(	$-$19.5)∗	1940(	$+$104)
	AP1roG-PT2MDd(sd)	$-$75.78630	1.238($-$0.014)	115.6(	$-$25.5)∗	1963(	$+$127)
	AP1roG-PT2SDo(d)	$-$75.81778	1.242($-$0.010)	156.5(	$+$15.4)	1919(	$+$83)
	AP1roG-PT2SDo(sd)	$-$75.81783	1.242($-$0.010)	156.4(	$+$15.3)	1919(	$+$83)
	AP1roG-PT2MDo(d)	$-$75.79032	1.231($-$0.021)	125.9(	$-$15.2)∗	2010(	$+$174)
	AP1roG-PT2MDo(sd)	$-$75.79108	1.230($-$0.022)	119.9(	$-$21.2)∗	2019(	$+$183)
	AP1roG-PT2b(d)	$-$75.78350	1.235($-$0.017)	127.2(	$-$13.9)∗	1938(	$+$102)
	AP1roG-PT2b(sd)	$-$75.79400	1.228($-$0.024)	113.0(	$-$28.1)∗	2049(	$+$213)
	AP1roG-PT2b(d\p)	$-$75.78381	1.231($-$0.021)	123.9(	$-$17.2)∗	2016(	$+$180)
	AP1roG-PT2b(sd\p)	$-$75.79370	1.228($-$0.024)	112.9(	$-$28.2)∗	2048(	$+$212)
	AP1roG-LCCD Boguslawski and Ayers (2015)	$-$75.81125	1.240($-$0.012)	139.3(	$-$1.8)	1916(	$+$80)
	AP1roG-LCCSD Boguslawski and Ayers (2015)	$-$75.81257	1.240($-$0.012)	143.0(	$+$1.9)	1926(	$+$90)
	NEVPT2 Tecmer et al. (2014)	$-$75.78829	1.244($-$0.008)	148.0(	$+$6.9)	1886(	$+$50)
	CR-CCSD(T)	$-$75.80484	1.242($-$0.010)	152.1(	$+$9.0)	1989(	$+$153)
	MRCI-SD Peterson et al. (1993)	$-$75.78079	1.252	141.1		1836
N2	AP1roG	$-$109.12686	1.087($+$0.083)	255.9(	$-$38.0)	2435(	$+$94)
	AP1roG-PT2SDd(d)	$-$109.37326	1.094($-$0.010)	249.8(	$-$31.9)	2301(	$-$40)
	AP1roG-PT2SDd(sd)	$-$109.37383	1.093($-$0.011)	233.4(	$+$15.5)	2299(	$-$42)
	AP1roG-PT2MDd(d)	$-$109.35918	1.104($+$0.000)	210.3(	$-$7.6)∗	2744(	$+$403)
	AP1roG-PT2MDd(sd)	$-$109.35988	1.105($+$0.001)	202.0(	$-$15.1)∗	2779(	$+$438)
	AP1roG-PT2SDo(d)	$-$109.39838	1.110($+$0.006)	246.9(	$+$29.0)	2278(	$-$63)
	AP1roG-PT2SDo(sd)	$-$109.39838	1.108($+$0.004)	248.8(	$+$30.9)	2254(	$-$87)
	AP1roG-PT2MDo(d)	$-$109.38106	1.089($-$0.015)	220.8(	$+$2.9)∗	2200(	$-$141)
	AP1roG-PT2MDo(sd)	$-$109.38132	1.090($-$0.014)	216.7(	$-$1.0)∗	2203(	$-$138)
	AP1roG-PT2b(d)	$-$109.37490	1.092($-$0.012)	219.2(	$+$1.3)∗	2246(	$-$95)
	AP1roG-PT2b(sd)	$-$109.37582	1.095($-$0.009)	207.1(	$-$10.8)∗	2272(	$-$69)
	AP1roG-PT2b(d\p)	$-$109.37431	1.092($-$0.012)	218.9(	$+$1.0)∗	2240(	$-$101)
	AP1roG-PT2b(sd\p)	$-$109.37525	1.094($-$0.010)	206.8(	$-$11.1)∗	2263(	$-$78)
	AP1roG-fLCCD	$-$109.39671	1.102($-$0.002)	206.7(	$-$11.2)	2334(	$-$7)
	AP1roG-fLCCSD	$-$109.39871	1.103($-$0.001)	211.6(	$-$6.3)	2337(	$-$4)
	CASSCF Peterson et al. (1993)	$-$109.13190	1.106($+$0.002)	211.6(	$-$6.3)	2340(	$+$1)
	CR-CCSD(T)	$-$109.39632	1.101($-$0.003)	229.7(	$+$11.8)∗	2397(	$+$56)
	MRCI-SD Peterson et al. (1993)	$-$109.36162	1.104	217.9		2341
BN	AP1roG	$-$79.07582	1.304($+$0.019)	114.3(	$-$40.1)	1630(	$-$52)
	AP1roG-PT2SDd(d)	$-$79.28333	1.283($-$0.002)	181.3(	$+$26.9)	1751(	$+$68)
	AP1roG-PT2SDd(sd)	$-$79.29473	1.271($-$0.014)	166.1(	$+$11.7)∗	1784(	$+$102)
	AP1roG-PT2MDd(d)	$-$79.26001	1.273($-$0.012)	160.7(	$+$6.3)	1778(	$+$96)
	AP1roG-PT2MDd(sd)	$-$79.27571	1.264($-$0.021)	144.5(	$-$9.9)	1760(	$+$78)
	AP1roG-PT2SDo(d)	$-$79.29432	1.282($-$0.003)	169.9(	$+$15.5)	1718(	$+$36)
	AP1roG-PT2SDo(sd)	$-$79.29427	1.281($-$0.004)	169.7(	$+$15.3)∗	1719(	$+$37)
	AP1roG-PT2MDo(d)	$-$79.26814	1.274($-$0.011)	158.5(	$+$4.1)∗	1743(	$+$61)
	AP1roG-PT2MDo(sd)	$-$79.27584	1.269($-$0.016)	139.0(	$-$15.4)∗	1725(	$+$43)
	AP1roG-PT2b(d)	$-$79.25960	1.273($-$0.012)	152.6(	$-$1.8)∗	1739(	$+$57)
	AP1roG-PT2b(sd)	$-$79.27961	1.260($+$0.025)	141.9(	$-$12.5)∗	1746(	$+$64)
	AP1roG-PT2b(d\p)	$-$79.25938	1.273($+$0.012)	156.6(	$+$2.2)∗	1741(	$+$59)
	AP1roG-PT2b(sd\p)	$-$79.27934	1.261($-$0.024)	142.3(	$-$12.1)∗	1746(	$+$64)
	AP1roG-LCCD Boguslawski and Ayers (2015)	$-$79.28205	1.277($-$0.008)	174.1(	$-$19.8)	1734(	$+$52)
	AP1roG-LCCSD Boguslawski and Ayers (2015)	$-$79.28975	1.275($-$0.010)	178.6(	$+$24.2)	1749(	$+$67)
	CCSD(T) Li and Paldus (2006)	—	1.244($-$0.008)	—		1886(	$+$50)
	CR-CCSD(T)	$-$79.28080	1.275($-$0.010)	164.9(	$+$10.5)	1726(	$+$44)
	RMR CCSD(T) Li and Paldus (2006)	—	1.284($-$0.001)	—		1681(	$+$1)
	MRCI-SD+Q Peterson (1995)	$-$79.26794	1.285	154.4		1682

Table 3. Table 3: List of reactions containing closed-shell molecules.

\ceF2 + H2 -¿ 2HF	{1}
\ceF2O + H2 -¿ F2 + H2O	{2}
\ceH2O2 + H2 -¿ 2H2O	{3}
\ceN2 + 3H2 -¿ 2NH3	{4}
\ceN2O + H2 -¿ N2 + H2O	{5}
\ceC2H2 + H2 -¿ C2H4	{6}
\ceBH3 + 3HF -¿ BF3 + 3 H2	{7}
\ceCO + H2O -¿ CO2 + H2	{8}
\ceCO + 3H2 -¿ CH4 + H2O	{9}
\ce2BH3 -¿ B2H6	{10}
\ce2H2O -¿ (H2O)2	{11}
\ce2HF -¿ (HF)2	{12}
\ceHCOOH -¿ CO2 + H2	{13}
\ceCO + CH4 -¿ CH3CHO	{14}
\ce2NH3 -¿ (NH3)2	{15}

Table 4. Table 4: Reaction energies of { 3 }–{ 3 } obtained from different theoretical models using a cc-pVQZ basis set (in kcal/mol). All dynamic energy corrections on top of AP1roG are determined for the optimized natural AP1roG basis. Only the CR-CCSD(T) reference energy and the energy difference with respect to CR-CCSD(T) ( E p CR − CCSD ​ ( T ) − E p method ) subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method (E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}}) is given in the Table. ME: mean error ( ME = ∑ p N ( E p CR − CCSD ​ ( T ) − E p method ) / N ME superscript subscript 𝑝 𝑁 subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method 𝑁 \mathrm{ME}=\sum_{p}^{N}(E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}})/N ); RMSE: root mean square error ( RMSE = ∑ p N ( E p CR − CCSD ​ ( T ) − E p method ) 2 / N RMSE superscript subscript 𝑝 𝑁 superscript subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method 2 𝑁 \mathrm{RMSE}=\sqrt{\sum_{p}^{N}(E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}})^{2}/N} ); MAE: mean absolute error ( MAE = ∑ p N | E p CR − CCSD ​ ( T ) − E p method | / N MAE superscript subscript 𝑝 𝑁 subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method 𝑁 \mathrm{MAE}=\sum_{p}^{N}|E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}}|/N ); max AE: maximum absolute error ( max AE = max ( { | E p CR − CCSD ​ ( T ) − E p method | } \mathrm{max\,AE}=\max{(\{|E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}}|}\} ).

	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	ME	RMSE	MAE	max AE
CR-CCSD(T)	$-135.8$	$-70.9$	$-88.1$	$-40.5$	$-83.2$	$-50.0$	$-92.5$	$-6.0$	$-66.5$	$-44.0$	$-4.7$	$-4.4$	$2.4$	$3.8$	$-2.7$	$-$	$-$	$-$	$-$
HF	$9.1$	$5.2$	$4.9$	$-6.6$	$16.7$	$0.9$	$0.3$	$-5.1$	$-9.0$	$-22.1$	$-1.0$	$-0.7$	$1.9$	$-8.2$	$-1.1$	$-1.0$	$8.7$	$6.2$	$22.1$
AP1roG	$-19.5$	$-5.9$	$-14.5$	$-26.8$	$23.7$	$0.2$	$-3.1$	$-10.9$	$-21.7$	$-31.9$	$1.5$	$-2.1$	$-3.3$	$-6.7$	$-2.2$	$-8.2$	$15.5$	$11.6$	$31.9$
AP1roG-MP2	$15.9$	$14.6$	$17.3$	$-10.5$	$28.3$	$-12.3$	$-20.1$	$-29.6$	$-3.7$	$-16.8$	$-10.3$	$-2.7$	$-6.4$	$-14.8$	$-2.3$	$-3.6$	$15.9$	$13.7$	$29.6$
PT2SDd(d)	$15.7$	$11.5$	$16.9$	$-14.7$	$21.3$	$-11.9$	$-8.5$	$-23.1$	$-9.4$	$-14.1$	$-11.6$	$-1.9$	$-2.4$	$-13.5$	$-1.7$	$-3.2$	$13.4$	$11.9$	$23.1$
PT2SDd(sd)	$15.0$	$11.4$	$16.8$	$-14.6$	$21.2$	$-11.8$	$-8.5$	$-23.7$	$-9.9$	$-13.9$	$-11.6$	$-1.9$	$-2.8$	$-13.4$	$-1.7$	$-3.3$	$13.4$	$11.9$	$23.7$
PT2MDd(d)	$15.6$	$11.1$	$16.5$	$-9.4$	$18.3$	$-8.5$	$-6.4$	$-19.4$	$-5.9$	$-12.9$	$-11.6$	$-1.5$	$-1.5$	$-11.8$	$-1.4$	$-1.9$	$11.6$	$10.1$	$19.4$
PT2MDd(sd)	$15.7$	$11.2$	$17.5$	$-8.5$	$18.4$	$-8.4$	$-6.7$	$-21.8$	$-8.0$	$-13.5$	$-12.1$	$-1.7$	$-2.3$	$-12.9$	$-1.6$	$-2.3$	$12.3$	$10.7$	$21.8$
PT2SDo(d)	$3.0$	$4.3$	$0.7$	$1.9$	$3.8$	$0.0$	$-9.9$	$-1.4$	$4.3$	$-9.0$	$1.7$	$-1.0$	$-2.7$	$2.2$	$-0.4$	$-0.2$	$4.2$	$3.1$	$9.9$
PT2SDo(sd)	$2.9$	$4.3$	$0.7$	$1.8$	$3.7$	$0.0$	$-9.8$	$-1.4$	$4.3$	$-9.0$	$1.7$	$-1.0$	$-2.6$	$2.3$	$-0.4$	$-0.2$	$4.2$	$3.1$	$9.8$
PT2MDo(d)	$7.0$	$6.7$	$3.4$	$4.2$	$3.3$	$0.7$	$-10.3$	$-0.2$	$4.9$	$-10.6$	$0.1$	$-1.9$	$-0.8$	$0.4$	$-1.7$	$0.3$	$5.1$	$3.7$	$10.6$
PT2MDo(sd)	$8.0$	$7.2$	$4.6$	$5.5$	$3.6$	$0.7$	$-10.8$	$-1.9$	$3.9$	$-11.6$	$-0.4$	$-2.1$	$-1.1$	$-0.7$	$-1.9$	$0.2$	$5.5$	$4.3$	$11.6$
PT2b(d)	$-1.4$	$-0.8$	$-1.9$	$-4.7$	$-4.7$	$-0.4$	$5.6$	$6.7$	$-3.3$	$-3.2$	$0.8$	$0.4$	$3.7$	$3.1$	$0.6$	$0.0$	$3.4$	$2.7$	$6.7$
PT2b(sd)	$-1.3$	$-0.6$	$-0.9$	$-3.5$	$-4.4$	$-0.2$	$5.4$	$3.9$	$-5.8$	$-3.9$	$0.2$	$0.2$	$2.6$	$2.1$	$0.3$	$-0.4$	$3.0$	$2.4$	$5.8$
PT2b(d\p)	$-1.3$	$-0.8$	$-1.9$	$-4.8$	$-4.4$	$-0.5$	$5.5$	$6.5$	$-3.1$	$-3.5$	$0.8$	$0.4$	$3.3$	$3.2$	$0.6$	$0.0$	$3.3$	$2.7$	$6.5$
PT2b(sd\p)	$-1.2$	$-0.6$	$-0.9$	$-3.5$	$-4.2$	$-0.3$	$5.3$	$3.7$	$-5.6$	$-4.2$	$0.2$	$0.1$	$2.2$	$2.2$	$0.3$	$-0.4$	$3.0$	$2.3$	$5.6$
AP1roG-LCCD	$-2.1$	$1.6$	$-0.8$	$-4.6$	$7.1$	$-2.2$	$-1.9$	$-4.6$	$-1.3$	$-3.9$	$0.4$	$-0.1$	$-3.5$	$1.2$	$0.3$	$-1.0$	$3.1$	$2.4$	$7.1$
AP1roG-LCCSD	$-1.8$	$-1.3$	$-0.8$	$-2.7$	$-2.0$	$-1.0$	$-0.1$	$1.0$	$0.4$	$-1.5$	$0.7$	$0.1$	$-2.0$	$2.8$	$0.4$	$-0.5$	$1.5$	$1.2$	$2.8$
MP2	$5.9$	$0.9$	$2.6$	$-2.7$	$-7.1$	$-2.7$	$2.4$	$4.7$	$-0.8$	$0.6$	$0.7$	$0.4$	$2.9$	$1.7$	$0.8$	$0.7$	$3.1$	$2.5$	$7.1$
BCCD	$1.4$	$0.8$	$0.4$	$-0.7$	$4.0$	$0.4$	$0.0$	$-2.0$	$-0.8$	$-3.7$	$0.3$	$0.2$	$-1.2$	$-0.2$	$0.4$	$0.0$	$1.6$	$1.1$	$4.0$
CCSD	$1.2$	$0.2$	$0.3$	$-1.0$	$3.2$	$0.2$	$0.4$	$-1.6$	$-1.2$	$-3.6$	$0.3$	$0.2$	$-0.9$	$-0.4$	$0.4$	$-0.1$	$1.4$	$1.0$	$3.6$
CCSD(T)	$-1.7$	$-2.5$	$-1.7$	$-2.7$	$-2.5$	$-0.6$	$1.7$	$0.7$	$-2.8$	$-0.8$	$0.6$	$0.4$	$0.4$	$0.3$	$0.7$	$-0.7$	$1.6$	$1.3$	$2.8$
LCCD	$-2.9$	$0.3$	$-1.1$	$0.0$	$2.2$	$0.8$	$-0.8$	$-2.1$	$1.0$	$-1.2$	$0.4$	$0.3$	$-1.8$	$0.9$	$0.5$	$-0.2$	$1.4$	$1.1$	$2.9$
LCCSD	$-5.8$	$-4.6$	$-2.9$	$-2.1$	$-3.9$	$0.2$	$1.4$	$0.5$	$-2.0$	$-0.4$	$0.5$	$0.4$	$-0.6$	$0.9$	$0.6$	$-1.2$	$2.5$	$1.8$	$5.8$

Table 5. Table 5: Reaction energies of { 3 }–{ 3 } obtained from different theoretical models using a cc-pVQZ basis set (in kcal/mol). All dynamic energy corrections on top of AP1roG are determined for canonical Hartree–Fock orbitals. Only the CR-CCSD(T) reference energy and the energy difference with respect to CR-CCSD(T) ( E p CR − CCSD ​ ( T ) − E p method ) subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method (E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}}) is given in the Table. ME: mean error ( ME = ∑ p N ( E p CR − CCSD ​ ( T ) − E p method ) / N ME superscript subscript 𝑝 𝑁 subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method 𝑁 \mathrm{ME}=\sum_{p}^{N}(E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}})/N ); RMSE: root mean square error ( RMSE = ∑ p N ( E p CR − CCSD ​ ( T ) − E p method ) 2 / N RMSE superscript subscript 𝑝 𝑁 superscript subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method 2 𝑁 \mathrm{RMSE}=\sqrt{\sum_{p}^{N}(E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}})^{2}/N} ); MAE: mean absolute error ( MAE = ∑ p N | E p CR − CCSD ​ ( T ) − E p method | / N MAE superscript subscript 𝑝 𝑁 subscript superscript 𝐸 CR CCSD T 𝑝 superscript subscript 𝐸 𝑝 method 𝑁 \mathrm{MAE}=\sum_{p}^{N}|E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}}|/N ); max AE: maximum absolute error ( max AE = max ( { | E p CR − CCSD ​ ( T ) − E p method | } \mathrm{max\,AE}=\max{(\{|E^{\mathrm{CR-CCSD(T)}}_{p}-E_{p}^{\mathrm{method}}|}\} ).

	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	{3}	ME	RMSE	MAE	max AE
CR-CCSD(T)	$-135.8$	$-70.9$	$-88.1$	$-40.5$	$-83.2$	$-50.0$	$-92.5$	$-6.0$	$-66.5$	$-44.0$	$-4.7$	$-4.4$	$2.4$	$3.8$	$-2.7$	$-$	$-$	$-$	$-$
HF	$9.1$	$5.2$	$4.9$	$-6.6$	$16.7$	$0.9$	$0.3$	$-5.1$	$-9.0$	$-22.1$	$-1.0$	$-0.7$	$1.9$	$-8.2$	$-1.1$	$-1.0$	$8.7$	$6.2$	$22.1$
AP1roG	$0.7$	$27.0$	$17.1$	$-47.6$	$21.2$	$-15.4$	$-41.8$	$-15.4$	$-41.1$	$-35.6$	$-23.9$	$-25.9$	$18.6$	$-29.6$	$-17.5$	$-14.0$	$27.9$	$25.2$	$47.6$
AP1roG-MP2	$-1.2$	$8.7$	$3.3$	$-9.6$	$-5.7$	$-3.4$	$2.4$	$4.9$	$-3.2$	$1.7$	$-0.8$	$-0.8$	$4.4$	$1.9$	$0.1$	$0.2$	$4.4$	$3.5$	$9.6$
PT2SDd(d)	$1.5$	$2.6$	$1.4$	$-5.6$	$-5.9$	$-2.8$	$5.4$	$4.6$	$-3.6$	$0.9$	$0.8$	$0.3$	$3.0$	$1.5$	$0.5$	$0.3$	$3.3$	$2.7$	$5.9$
PT2MDd(d)	$2.3$	$-2.4$	$-0.6$	$-3.4$	$-6.5$	$-2.5$	$8.1$	$4.5$	$-4.4$	$0.0$	$2.1$	$1.2$	$2.3$	$1.1$	$0.6$	$0.2$	$3.6$	$2.8$	$8.1$
PT2MDd(sd)	$2.3$	$-2.1$	$0.0$	$-3.7$	$-6.3$	$-2.9$	$8.2$	$4.2$	$-4.8$	$-0.1$	$1.3$	$1.2$	$2.5$	$0.9$	$0.2$	$0.1$	$3.6$	$2.7$	$8.2$
AP1roG-LCCD	$0.0$	$4.0$	$1.5$	$1.3$	$2.9$	$1.3$	$-2.9$	$-2.0$	$3.7$	$0.4$	$-0.1$	$-0.2$	$-2.2$	$2.2$	$0.9$	$0.7$	$2.1$	$1.7$	$4.0$
AP1roG-LCCSD	$-1.7$	$0.1$	$-0.2$	$-0.4$	$-1.5$	$0.6$	$-0.6$	$-0.3$	$1.1$	$1.1$	$0.1$	$0.0$	$-1.3$	$1.8$	$1.0$	$0.0$	$1.0$	$0.8$	$1.8$
MP2	$5.9$	$0.9$	$2.6$	$-2.7$	$-7.1$	$-2.7$	$2.4$	$4.7$	$-0.8$	$0.6$	$0.7$	$0.4$	$2.9$	$1.7$	$0.8$	$0.7$	$3.1$	$2.5$	$7.1$
BCCD	$1.4$	$0.8$	$0.4$	$-0.7$	$4.0$	$0.4$	$0.0$	$-2.0$	$-0.8$	$-3.7$	$0.3$	$0.2$	$-1.2$	$-0.2$	$0.4$	$0.0$	$1.6$	$1.1$	$4.0$
CCSD	$1.2$	$0.2$	$0.3$	$-1.0$	$3.2$	$0.2$	$0.4$	$-1.6$	$-1.2$	$-3.6$	$0.3$	$0.2$	$-0.9$	$-0.4$	$0.4$	$-0.1$	$1.4$	$1.0$	$3.6$
CCSD(T)	$-1.7$	$-2.5$	$-1.7$	$-2.7$	$-2.5$	$-0.6$	$1.7$	$0.7$	$-2.8$	$-0.8$	$0.6$	$0.4$	$0.4$	$0.3$	$0.7$	$-0.7$	$1.6$	$1.3$	$2.8$
LCCD	$-2.9$	$0.3$	$-1.1$	$0.0$	$2.2$	$0.8$	$-0.8$	$-2.1$	$1.0$	$-1.2$	$0.4$	$0.3$	$-1.8$	$0.9$	$0.5$	$-0.2$	$1.4$	$1.1$	$2.9$
LCCSD	$-5.8$	$-4.6$	$-2.9$	$-2.1$	$-3.9$	$0.2$	$1.4$	$0.5$	$-2.0$	$-0.4$	$0.5$	$0.4$	$-0.6$	$0.9$	$0.6$	$-1.2$	$2.5$	$1.8$	$5.8$

Equations

$$|AP1roG⟩=\exp \left(\sum _{i=1}^P\sum _{a=P+1}^K{c}_i^a{a}_a^{†}{a}_{\bar{a}}^{†}{a}_{\bar{i}}{a}_i\right)|0⟩,$$

$${\hat{H}}_N^{\prime }$$

$$=({\hat{H}}_0-⟨0|{\hat{H}}_0|0⟩)+(\hat{V}-⟨0|\hat{V}|0⟩-E_{\mathrm{corr}}^{(0)})$$

$$=({\hat{H}}_0{)}_N+{\hat{V}}_N^{\prime },$$

$$\hat{H}={\hat{H}}_1+\hat{W}=\sum _{pq}{h}_{pq}{a}_p^{†}{a}_q+\frac{1}{2}\sum _{pqrs}⟨pq|rs⟩a_p^{†}{a}_q^{†}{a}_s{a}_r.$$

$$|{\Psi }^{(1)}⟩=\sum _p{t}_p|{\Phi }_p⟩,$$

$$⟨{\Psi }^{(1)}|AP1roG⟩=0.$$

$$({\hat{H}}_0{)}_N|{\Psi }^{(1)}⟩+{\hat{V}}_N^{\prime }|AP1roG⟩$$

$$\sum _p{t}_p({\hat{H}}_0{)}_N|{\Phi }_p⟩+{\hat{V}}_N^{\prime }|AP1roG⟩$$

$$E^{(0)}+E^{(1)}=\frac{1}{⟨\tilde{\Psi }|AP1roG⟩}⟨\tilde{\Psi }|{\hat{V}}_N^{\prime }|AP1roG⟩=0,$$

$$E^{(2)}=\frac{⟨\tilde{\Psi }|{\hat{V}}_N^{\prime }|{\Psi }^{(1)}⟩}{⟨\tilde{\Psi }|AP1roG⟩}.$$

$$|{\Psi }^{(1)}⟩=\hat{T}|0⟩,$$

$${\hat{T}}_2^{\prime }=\frac{1}{2}\sum _{ij}^{\mathrm{occ}}\sum _{ab}^{\mathrm{virt}}{}^{\prime }{t}_{ij}^{ab}{\hat{E}}_{ai}{\hat{E}}_{bj},$$

$${\hat{T}}_1=\sum _i^{\mathrm{occ}}\sum _a^{\mathrm{virt}}{t}_i^a{\hat{E}}_{ai}$$

$$⟨\overline{{}_{ij}^{ab}}|=\frac{1}{3}{⟨}_{ij}^{ab}|+\frac{1}{6}{⟨}_{ji}^{ab}|,$$

$$⟨\overline{{}_i^a}|=\frac{1}{2}{⟨}_i^a|,$$

$$⟨\overline{{}_{ij}^{ab}}{|}_{kl}^{cd}⟩={\delta }_{iajb,kcld}+{\delta }_{jbia,kcld}$$

$$⟨\overline{{}_i^a}{|}_j^b⟩={\delta }_{ia,jb}.$$

$$\sum _p{t}_p⟨\overline{{\Phi }_q}|({\hat{H}}_0{)}_N|{\Phi }_p⟩+⟨\overline{{\Phi }_q}|{\hat{V}}_N^{\prime }|AP1roG⟩=0,$$

$$\hat{F}$$

$$=\sum _{pq}{f}_{pq}{a}_p^{†}{a}_q,$$

$$({\hat{H}}_0{)}_N={\hat{F}}_N^{\mathrm{d}}=\sum _p{f}_{pp}\{a_p^{†}{a}_p\},$$

$${\hat{V}}_N^{\prime }$$

$$=\sum _{p\ne q}{f}_{pq}\{a_p^{†}{a}_q\}+\frac{1}{2}\sum _{pqrs}{}^{\prime }⟨pq|rs⟩\{a_p^{†}{a}_q^{†}{a}_s{a}_r\}$$

$$E^{(2)}=⟨0|{\hat{V}}_N^{\prime }|{\Psi }^{(1)}⟩.$$

$$E_{\mathrm{d}}^{(2)}=\sum _{iajb}{t}_{ij}^{ab}(⟨ij||ab⟩+⟨ij|ab⟩),$$

$$E_{\mathrm{sd}}^{(2)}=2\sum _{ia}{f}_{ia}{t}_i^a+\sum _{iajb}{t}_{ij}^{ab}(⟨ij||ab⟩+⟨ij|ab⟩).$$

$$({\hat{H}}_0{)}_N$$

$$=\sum _p{\bar{f}}_{pp}\{a_p^{†}{a}_p\}.$$

$${\hat{V}}_N^{\prime }$$

$$=\sum _{p,q}(f_{pq}-{\bar{f}}_{pp}{\delta }_{pq})\{a_p^{†}{a}_q\}+\frac{1}{2}\sum _{pqrs}{}^{\prime }⟨pq|rs⟩\{a_p^{†}{a}_q^{†}{a}_s{a}_r\}.$$

$$|\overline{{\Psi }^{(1)}}⟩=\sum _p\frac{t_p}{⟨AP1roG|AP1roG⟩}|{\Phi }_p⟩=\sum _p{\bar{t}}_p|{\Phi }_p⟩.$$

$$E^{(2)}=\sum _p{\bar{t}}_p⟨AP1roG|{\hat{V}}_N^{\prime }|{\Phi }_p⟩.$$

$$({\hat{H}}_0{)}_N={\hat{F}}_N^{\mathrm{d}}+{\hat{F}}_N^{\mathrm{o}}=\sum _{p,q}{f}_{pq}\{a_p^{†}{a}_q\}.$$

Full text

I Introduction

Accurate and reliable theoretical predictions of thermodynamic properties of chemical reactions remain an active field of research Peterson et al. (2012). Computational studies are particularly desirable especially when experimental manipulations are limited or impossible. Conventional “ab initio” quantum chemistry tools used in thermochemistry are commonly based on different flavours of many-body perturbation theory and coupled cluster approaches Coester (1958); Cizek (1966); Cizek and Paldus (1971); Paldus et al. (1972); Bartlett (1981); Helgaker et al. (2000); Shavitt and Bartlett (2009); Bartlett and Musiał (2007). Over the years, the CCSD(T) (Coupled Cluster Singles, Doubles and perturbative Triples) method has unfolded as one of the most reliable and accurate models for systems that are well-represented by a single electron configuration. Therefore, it is often referred to as the “gold standard” of quantum chemistry. Although CCSD(T) performs extraordinary well for systems dominated by dynamic correlation, the model often fails for multi-reference systems. For molecular systems with a substantial multi-reference character, multi-reference methods Szalay et al. (2012) are commonly used, which are computationally very expensive. Furthermore, the most popular multi-reference approaches used in quantum chemistry do not a priori account for the dynamic part of the correlation energy. Exceptions are multi-reference coupled cluster models, some variants of the multi-reference configuration interaction method, and the perturb-then-diagonalize multi-reference perturbation theory approaches that use effective Hamiltonians. To remedy this problem, a posteriori corrections have been developed that are usually based on perturbation theory Andersson et al. (1990, 1992); Angeli et al. (2001); Pulay (2011); Malmqvist et al. (2008), canonical transformation theory Neuscamman et al. (2010) or coupled cluster methods Lyakh et al. (2012).

Although multi-reference methods are frequently applied to model strongly-correlated systems, like bond-breaking processes and heavy-element chemistry, some flavours of single-reference coupled cluster theory pose promising alternatives to describe static/nondynamic electron correlation. Among the most representative examples are the CR-CC (Completely Renormalized Coupled Cluster) and CR-EOMCC (Completely Renormalized Equation of Motion Coupled Cluster) approaches Piecuch et al. (2002, 2004); Piecuch and Wloch (2005); Wloch et al. (2007); Kowalski and Piecuch (2004), the active-space CC and EOMCC methods Kowalski and Piecuch (2000); Kowalski et al. (2005, 2010); Piecuch (2010), the EA/IP- (Electron Affinity/Ionization Potential) and DEA/DIP-(Double EA/IP) EOMCC models Shen and Piecuch (2013), and the spin-flip CC/EOMCC formalism Krylov (2006); Levchenko and Krylov (2004); Krylov (2008). A different, computationally feasible approach suitable for strongly-correlated systems uses seniority-zero wavefunctions to describe the static/nondynamic part of the electron correlation energy. Boguslawski et al. (2014a); Tecmer et al. (2014); Boguslawski et al. (2014b, c); Tecmer et al. (2015); Rassolov (2002); Surján et al. (2012); Neuscamman (2012); Limacher et al. (2013); Ellis et al. (2013); Limacher et al. (2014a); Stein et al. (2014); Henderson et al. (2014a); Neuscamman (2016) The missing dynamic electron correlation effects are included a posteriori in these ansätze using, for instance, many-body perturbation theory Jeszenszki et al. (2014); Limacher et al. (2014b); Limacher (2015), coupled-cluster theory Zoboki et al. (2013); Small et al. (2014); Henderson et al. (2014b); Boguslawski and Ayers (2015); Veis et al. (2016), extended random phase approximation Pernal (2014), and density functional theory (DFT) corrections Garza et al. (2015a, b).

In this article, we seek computationally cheap and reliable dynamic electron correlation models on top of the Antisymmetric Product of 1-reference orbital Geminal (AP1roG) of different levels of approximation. The AP1roG wavefunction can be written using an exponential ansatz,

[TABLE]

where the sum runs over all electron pairs $P$ (equal to the number of occupied orbitals) and virtual orbitals $K-P$, while $|0\rangle$ is some reference determinant. $a_{p}^{\dagger}$ and $a_{p}$ are the fermionic creation and annihilation operators for spin-up $p$ and spin-down $\bar{p}$ electrons, while $c_{i}^{a}$ are the AP1roG amplitudes, also called geminal coefficients. The geminal coefficients $c_{i}^{a}$ are obtained by solving the projected Schrödinger equation, as in coupled cluster theory. The AP1roG model allows us to describe static/nondynamic correlation and to a minor extent some dynamic correlation (see also ref. Boguslawski et al. (2016)). To include the missing part of the dynamic electron correlation effects a posteriori that go beyond the AP1roG ansatz, that is, beyond electron pair states, we will use perturbation theory as well as a linearized coupled cluster corrections. Boguslawski and Ayers (2015) Specifically, our perturbation theory models will combine a diagonal and off-diagonal zero-order Hamiltonian, a single- and multi-determinant wavefunction as dual, and different excitation operators used to construct the projection manifold. The performance of all perturbation theory models as well as the linearized coupled cluster correction, as presented in ref. 51, will be assessed against spectroscopic constants of homo- and heteronuclear dimers and reaction energies of closed shell systems.

This article is organized as follows. In section II, we present different theoretical models to correct for the missing dynamic electron correlation in the AP1roG ansatz. Computational details are described in section III. Numerical results and comparison to standard electron correlation methods for molecular systems dominated by static/nondynamic (bond dissociation processes) and dynamic electron correlation (reaction energies of main group compounds) are presented in section IV. Finally, we conclude in section V.

II Theoretical models for dynamic correlation

A common way to account for dynamic correlation effects in quantum many-body systems is to use perturbation theory (PT). Specifically for the AP1roG wavefunction, two different PT models have been proposed that allow us to describe electron correlation effects beyond electron pairs. Limacher et al. (2014b) While these approaches provide reliable spectroscopic constants for some first-row diatomic molecules, Tecmer et al. (2014); Boguslawski and Ayers (2015) their performance deteriorates when moving to heavy-element containing compounds like actinide species. Boguslawski and Ayers (2015) To remedy this problem, a linearized coupled cluster (LCC) corrections on top of AP1roG can be used that allows us to accurately describe closed-shell molecules with (quasi-)degenerate $f$-, $d$-, $p$-, and $s$-shells, as encountered in actinide-containing compounds. Boguslawski and Ayers (2015) In the following, we will briefly summarize some theoretical models for dynamic correlation with an AP1roG reference function. Specifically, we will extend the existing PT models and relate the proposed models to the recently presented PTa and PTb methods. Limacher et al. (2014b)

II.1 2nd-order Perturbation Theory

One drawback of multi-reference perturbation theory is the arbitrariness of the theoretical model. For example, there are different choices for the zero-order Hamiltonian $\hat{H}_{0}$, the dual state $\langle\tilde{\Psi}|$ in the projector, and the choice of the projection space. Pulay (2011) Poor choices may lead to technical difficulties and unphysical solutions. In this work, we will investigate different selections for the zero-order Hamiltonian $\hat{H}_{0}$, the dual state $\langle\tilde{\Psi}|$, and the projection space, while the zero-order wavefunction is restricted to the AP1roG reference function of eq. (1), $|\Psi^{(0)}\rangle=|\mathrm{AP1roG}\rangle$.

Furthermore, in the derivation of all PT models, we will use the quantum chemical Hamiltonian in its normal product from shifted by the correlation energy of AP1roG $E_{\mathrm{corr}}^{(0)}$ (the shift in energy is indicated by “′”),

[TABLE]

where the quantum chemical Hamiltonian $\hat{H}$ is divided into a zero-order contribution $\hat{H}_{0}$ and a perturbation $\hat{V}$. It is convenient to rewrite $\hat{H}$ into a sum of a one- ($\hat{H}_{1}$) and a two-electron ($\hat{H}_{2}$, here indicated by $\hat{W}$) part,

[TABLE]

In the above equation, $h_{pq}$ and $\langle pq|rs\rangle$ are the one- and two-electron (written in physicists’ notation) integrals, respectively, determined for the one-particle basis functions $p,q,r,s$. We should note that, in this work, we will restrict $\hat{H}_{0}$ to be a one-body operator so that the normal-product form of the perturbation $\hat{V}=\hat{H}-\hat{H}_{0}$ can be written as an operator shifted by the AP1roG correlation energy, $\hat{V}_{N}^{\prime}=\hat{V}_{N}-E_{\mathrm{corr}}^{(0)}$ (again indicated by “′”), as we have $E_{\mathrm{corr}}^{(0)}=\langle 0|\hat{W}_{N}|\mathrm{AP1roG}\rangle$. Introducing a shifted perturbation operator $\hat{V}^{\prime}_{N}$ is equivalent to neglecting contractions (or diagrams) that correspond to the AP1roG correlation energy in the PT equations, which will be indicated by the ′ in the sum of the $\hat{W}_{N}^{\prime}$ operator.

As in conventional Rayleigh–Schrödinger PT, the exact wavefunction can be written as an order-by-order expansion, $|\Psi\rangle=|\mathrm{AP1roG}\rangle+\lambda|\Psi^{(1)}\rangle+\ldots$, where $\lambda$ is the order parameter. The first-order correction to the wavefunction is expanded in a set of Slater determinants $\Phi_{p}$,

[TABLE]

and forced to be orthogonal to the zero-order wavefunction, here, $|\mathrm{AP1roG}\rangle$,

[TABLE]

This orthogonality constraint restricts the choice of the projection space used for the expansion of $|{\Psi^{(1)}}\rangle$ (and higher orders). By construction, all pair-excitations with respect to $|0\rangle$ have to be excluded as they don’t satisfy $\langle\Phi_{p}|\mathrm{AP1roG}\rangle=0$. Similarly, introducing an order parameter $\lambda$ in the Hamiltonian $\hat{H}=\hat{H}_{0}+\lambda\hat{V}$ and equating coefficients of powers of $\lambda$, we obtain the zero-, first-, and higher-order PT equations. Specifically for the first-order correction to the wavefunction, we have to solve

[TABLE]

Since we have introduced a shifted normal-product Hamiltonian, the zero- and first-order energy corrections vanish,

[TABLE]

where $\langle\tilde{\Psi}|$ is the dual of the unperturbed state $|\mathrm{AP1roG}\rangle$. Specific choices for $\langle\tilde{\Psi}|$ will be considered below. The first non-zero correction to the energy is of second-order and can be calculated from the first-order wavefunction and the shifted normal-product perturbation Hamiltonian,

[TABLE]

Before we focus on possible choices of $(\hat{H}_{0})_{N}$ and $\hat{V}_{N}^{\prime}$ as well as $\langle\tilde{\Psi}|$, we will define our projection space used in the expansion of $|\Psi^{(1)}\rangle$ in eq. (4). Following previous PT models as well as linear CC corrections with an AP1roG reference function, the projection space will contain all possible excitations with respect to a reference determinant. This reference determinant is, however, not arbitrary, but restricted to the reference determined of AP1roG, $|0\rangle$ of eq. (1). We should emphasize that $|0\rangle$ is not equivalent to the Hartree–Fock determinant as in conventional CC theory, but adjusted during the optimization of the AP1roG wavefunction. Choosing $|0\rangle$ as reference determinant, the first-order correction can be written as

[TABLE]

where $\hat{T}$ is some excitation operator that substitutes electrons from the occupied to the virtual space with respect to $|0\rangle$. Furthermore, we will restrict $\hat{T}$ to contain double excitations without electron pairs, $\hat{T}=\hat{T}_{2}^{\prime}$, as well as single and double excitations, $\hat{T}=\hat{T}_{1}+\hat{T}_{2}^{\prime}$. If only double excitations are included, the excitation operator is specified as

[TABLE]

where $\hat{E}_{ai}=a^{\dagger}_{a}a_{i}+a^{\dagger}_{\bar{a}}a_{\bar{i}}$ is the singlet excitation operator and the perturbation amplitudes $t_{ij}^{ab}$ are symmetric with respect to pair-exchange, i.e., $t_{ij}^{ab}=t_{ji}^{ba}$. Similar to our notation of the shifted normal-product Hamiltonian, the prime in the above summation indicates that pair-excited determinants are excluded in the excitation operator, i.e., $t_{i\bar{i}}^{a\bar{a}}=0$. Exclusion of pair-excitations fulfils the orthogonality condition eq. (5). If the excitation operator contains both single and double excitations, the single excitations can be accounted for by adding

[TABLE]

to the double excitation operator. Furthermore, as basis for the bra states of the projection manifold, we follow refs. 46; 51 and use the convenient choice

[TABLE]

for doubly excited determinants where $\langle^{ab}_{ij}|=\langle 0|\hat{E}_{jb}\hat{E}_{ia}$ and

[TABLE]

for singly excited determinants with $\langle^{a}_{i}|=\langle 0|\hat{E}_{ia}$. The bra basis then forms a biorthogonal basis that satisfies the normalization conditions

[TABLE]

for doubly-excited determinants, while for the singly-excited determinants we have

[TABLE]

Choosing eq. (4) as the first-order correction to the wavefunction with excitation operators as defined in eqs. (10) and (11), the corresponding perturbation amplitudes are then obtained by solving

[TABLE]

Note that the above equations depend on the partition of the Hamiltonian $\hat{H}$ into the zero-order and the perturbation part as well as the projection manifold, but not on the choice of the dual state $\langle\tilde{\Psi}|$. In the following, we will consider different partitionings of $\hat{H}$ as well as two choices for $\langle\tilde{\Psi}|$.

II.1.1 Diagonal one-body zero-order Hamiltonian

First, we will restrict $(\hat{H}_{0})_{N}$ to be a diagonal one-electron operator. Specifically, we will choose $(\hat{H}_{0})_{N}$ to be the diagonal of the inactive Fock operator

[TABLE]

where $\langle pi||qi\rangle$ are the two-electron integrals in physicists’ notation containing the Coulomb $\langle pi|qi\rangle$ and exchange $\langle pi|iq\rangle$ terms. Using normal-product operators, the zero-order Hamiltonian reads

[TABLE]

while the perturbation becomes (cf. eq. (II.1))

[TABLE]

Eqs. (18) and (II.1.1) are then substituted into eq. (16) to solve for the PT amplitudes. Note that, in the case of a diagonal zero-order Hamiltonian, the PT amplitudes are obtained from a set of uncoupled equations.

Restricting the dual $\langle\tilde{\Psi}|$ to $\langle 0|$.

If $\langle\tilde{\Psi}|$ is restricted to the AP1roG reference determinant $\langle 0|$, we can straightforwardly evaluate the overlap $\langle\Phi_{0}|\mathrm{AP1roG}\rangle$, which equals 1 due to intermediate normalization of the AP1roG wavefunction. The expression for the second-order energy correction $E^{(2)}$ given in eq. (8) thus simplifies to

[TABLE]

Specifically, the energy correction for $\hat{T}=\hat{T}^{\prime}_{2}$ is given as

[TABLE]

while for single and double excitations we have

[TABLE]

We should emphasize that this PT model is equivalent to the PTa model presented in ref. 46. Note, however, that pair excitations are not excluded in ref. 46 and that the full Hamiltonian $\hat{H}$ is taken as the perturbation Hamiltonian so that $E^{(2)}=\langle 0|\hat{H}|\Psi^{(1)}\rangle$. Since the PTa amplitude equations for the pair excitations vanish (as they equal the AP1roG amplitude equations), the PTa energy corrections is nonetheless determined from eq. (20). For reasons of consistency, we will abbreviate these PT models using a single determinant (SD) as dual and a diagonal (d) zero-order Hamiltonian as PT2SDd. The choice of the excitation operator will be indicated in parentheses, that is, PT2SDd(d) for double excitations and PT2SDd(sd) for both single and double excitations.

Choosing $\langle\mathrm{AP1roG}|$ as dual $\langle\tilde{\Psi}|$.

If $\langle\mathrm{AP1roG}|$ is taken as the dual state $\langle\tilde{\Psi}|$, we have to evaluate terms as $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ in the energy expression, which becomes computationally intractable for large systems. In order to arrive at a computationally feasible model, we will follow ref. 46 to, at least partially, eliminate the overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ in the PT equations and energy expression. For that purpose, we redefine the zero-order Hamiltonian of eq. (18) by introducing the inverse of the overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ as a scaling factor,

[TABLE]

By changing the zero-order Hamiltonian of eq. (18), we also have to adjust the corresponding perturbation part given in eq. (II.1.1),

[TABLE]

To fully avoid the evaluation of the overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ in the PT amplitude equations eq. (16), the inverse of the AP1roG wavefunction overlap will be absorbed in the PT amplitudes. Thus, the first-order wavefunction contains scaled amplitudes,

[TABLE]

By substituting eq. (II.1.1) into eq. (16) and introducing the scaled PT amplitudes from eq. (25), we can eliminate the wavefunction overlap from the PT working equations. Furthermore, scaling the PT amplitudes by the inverse of the overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ restores the zero-order Hamiltonian of eq. (18) and we get $(\hat{H}_{0})_{N}=\hat{{F}}_{N}^{\mathrm{d}}$. Note that the wavefunction overlap is still present in the perturbation part $\hat{V}_{N}^{\prime}$. Due to the structure of the zero-order wavefunction and the choice of the projection manifold, the diagonal part of the modified Fock operator in eq. (II.1.1) does not contribute to the PT amplitude equations and the resulting perturbation reduces to $\hat{V}_{N}^{\prime}=\hat{{F}}_{N}^{\mathrm{o}}+\hat{W}_{N}^{\prime}$. Since we have chosen $\langle\mathrm{AP1roG}|$ as dual, the second-order energy correction is determined from $E^{(2)}=\langle\mathrm{AP1roG}|\hat{V}_{N}^{\prime}|\overline{\Psi^{(1)}}\rangle$ (cf. eq. (8)), where the first-order correction to the wavefunction is calculated from the scaled PT amplitudes (eq. (25)),

[TABLE]

The sum in the above equation runs over all determinants in the projection manifold (doubly excited or singly- and doubly-excited determinants). We should note that although we can exactly evaluate the energy correction and the PT equations, the zero- and first-order energy corrections eq. (7) do not vanish. However, we can neglect the weights of the (AP1roG) amplitudes beyond single pair excitations and assume that $E^{(0)}+E^{(1)}\approx 0$. Limacher et al. (2014b) We will label this PT model as PT2MDd as it uses a multi-determinant (MD) wavefunction as dual and a diagonal (d) zero-order Hamiltonian. Furthermore, PT2MDd(d) indicates that the excitation operator contains only double excitations (without pairs), while in PT2MDd(sd) both single and double excitations are included in $\hat{T}$.

II.1.2 Off-diagonal one-body zero-order Hamiltonian

Next, we will consider an off-diagonal one-electron operator as the zero-order Hamiltonian. Similar to the PT models employing a diagonal zero-order Hamiltonian, we choose the Fock operator of eq. (II.1.1) as our $\hat{H}_{0}$ Hamiltonian,

[TABLE]

Then, the perturbation part $\hat{V}=\hat{H}-\hat{H}_{0}$ in its (shifted) normal-product form reads

[TABLE]

To obtain the working equation for the PT amplitudes, we substitute eqs. (27) and (28) into eq. (16). Note that, in contrast to the PT methods with a diagonal zero-order Hamiltonian, the PT amplitudes are now obtained from a set of coupled equations.

Restricting the dual $\langle\tilde{\Psi}|$ to $\langle 0|$.

If the dual state is restricted to the reference determinant of the AP1roG wavefunction, we benefit from the intermediate normalization when evaluating the overlap $\langle 0|\mathrm{AP1roG}\rangle$. Analogous to PT2SDd-type methods, the second order energy can be evaluated from eq. (20). Note, however, that only the doubly-excited determinants directly contribute to the energy correction. Since the perturbation Hamiltonian is a two-electron operator, the second-order energy correction of the single excitations vanishes. Single excitations contribute indirectly through coupling to the double excitation manifold in the PT amplitude equations. For both including and excluding the singles projection manifold, the second-order energy $E^{(2)}$ is thus calculated from eq. (21). We will abbreviate the PT corrections using an off-diagonal (o) zero-order Hamiltonian and a single determinant for its dual state as PT2SDo, while the projection manifold will be indicated in parentheses (d for doubles, sd for singles and doubles, respectively).

Choosing $\langle\mathrm{AP1roG}|$ as dual $\langle\tilde{\Psi}|$.

Similar to the PT methods with a diagonal $\hat{H}_{0}$ Hamiltonian, choosing $\langle\mathrm{AP1roG}|$ as dual forces us to evaluate the wavefunction overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$, which becomes prohibitive for large systems. In order to (partially) remove the wavefunction overlap from the working equations, we follow the procedure from above and introduce a scaled zero-order Hamiltonian, where we have to modify both the diagonal ${{F}}_{N}^{\mathrm{d}}$ and off-diagonal ${{F}}_{N}^{\mathrm{o}}$ Fock operator,

[TABLE]

Since we use a modified Fock operator as zero-order Hamiltonian, we have to account for it in the definition of the perturbation part,

[TABLE]

Analogous to PT2MDd-type methods, the scaling factor in the $(\hat{H}_{0})_{N}$ Hamiltonian is absorbed in the PT amplitudes (cf.eq. (25)) so that the zero-order Hamiltonian in the PT amplitude equations eq. (16) contains only the unscaled Fock operator, $(\hat{H}_{0})_{N}=\hat{{F}}_{N}^{\mathrm{d}}+\hat{{F}}_{N}^{\mathrm{o}}$, while the PT amplitudes $t_{p}$ are replaced by the scaled amplitudes $\bar{t}_{p}$. Although we eliminated the wavefunction overlap in the second-order energy expression and in the $(\hat{H}_{0})_{N}$ part of the PT amplitude equations, $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ still remains in the perturbation part (cf. eq. (II.1.2)). In contrast to the PT2MDd-type methods discussed above, the perturbation $\hat{V}_{N}^{\prime}$ contains also modified off-diagonal elements in the one-electron Fock operator that do not vanish in the PT equations. Therefore, we have to evaluate the wavefunction overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle$ before we can determine the PT amplitudes. To obtain a computationally feasible model, we will approximate the overlap $\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle\approx 1+\sum_{ia}|c_{i}^{a}|^{2}$, keeping only the quadratic terms in the AP1roG amplitudes. This is usually a good approximation as the AP1roG wavefunction amplitudes are typically much smaller than 1 ($c_{i}^{a}\ll 1$). Thus, the scaled off-diagonal Fock matrix elements can be approximated by $\bar{f}_{pq}\approx\frac{f_{pq}}{1+\sum_{ia}|c_{i}^{a}|^{2}}$. As for PT2MDd, the second-order energy $E^{(2)}$ can be determined from eq. (26) with $\hat{V}_{N}^{\prime}$ defined in eq. (II.1.2). Note that in contrast to PT2SDo, the single excitations directly contribute to the energy correction through both the one- and two-electron operators in the perturbation Hamiltonian. The PT models with an off-diagonal (o) $\hat{H}_{0}$ Hamiltonian and a multi-determinant (MD) wavefunction ($\langle\mathrm{AP1roG}|$) as dual will be abbreviated as PT2MDo, while the projection manifold will be again indicated in parentheses (d for doubles, sd for singles and doubles, respectively). Note that pair-excitations are excluded in the projection manifold due to the orthogonality constraint.

Relation to PTb theory.

The PT2MDo model is similar, but not equivalent, to the recently presented PTb approach. Limacher et al. (2014b) In the PTb method, the wavefunction overlap in the perturbation Hamiltonian is neglected, that is, we assume $\frac{1}{\langle\mathrm{AP1roG}|\mathrm{AP1roG}\rangle}\rightarrow 0$ so that $\bar{f}_{pq}\rightarrow 0$. By neglecting the scaled components of the Fock operator, the perturbation Hamiltonian reduces to the full quantum chemical Hamiltonian $\hat{V}_{N}^{\prime}=\hat{H}_{N}^{\prime}$. The second-order energy correction is then given as

[TABLE]

Furthermore, in PTb theory, pair excitations are not excluded in the projection manifold and the first-order correction of the wavefunction contains all double excitations. These pair excitations do not contribute to the energy corrections eq. (31). Their contribution vanishes as the corresponding terms in eq. (31) equal the AP1roG amplitude equations. However, pair excitations indirectly enter the energy correction by coupling to the remaining PT amplitudes in the PT equations as well as to the pair excitations of the AP1roG model. In this work, we have extended the original PTb model by including also single excitations in the projection manifold. Furthermore, we also investigate how the pair excitations in the projection manifold influence the PTb energy correction. For that purpose we have excluded the pair excitations in the projection manifold when optimizing the PTb amplitudes. To emphasize the order of the energy correction in PTb, these PT models are abbreviated as PT2b, while the projection manifold is indicated in parentheses (d for doubles, sd for singles and doubles, d\p for doubles without pairs, sd\p for singles and doubles excluding pair excitations).

Comment on computational scaling of the PT models.

The computational scaling of all PT corrections discussed above is determined by the first term of eq. (16). Note that we can introduce suitable intermediates for the second term in eq. (16) so that summations are performed only once in the beginning of the calculation. If $(\hat{H}_{0})_{N}$ is restricted to be a diagonal one-body operator, the computational cost of the corresponding PT models scales as $\mathcal{O}(o^{2}v^{2})$, where $o$ is the number of occupied and $v$ the number of virtual orbitals, respectively. Choosing an off-diagonal one-body zero-order Hamiltonian $(\hat{H}_{0})_{N}$ increases the computational cost to $\mathcal{O}(o^{2}v^{3})$ (see also Table 1). Since we now have to solve a coupled set of linear equations iteratively, we have to consider an additional prefactor. However, this prefactor is typically much smaller than $v$. Thus, PT2SDd and PT2MDd scale similar to AP1roG (or conventional electronic structure methods like MP2), while in PT2SDo, PT2MDo, and PT2b the computational cost increases by a factor of $v$. All PT approaches presented in this work are summarized in Table 1 for comparison.

II.2 Combining AP1roG with MP2

In our last PT model, we will combine conventional MP2 theory with AP1roG. Specifically, all wavefunction amplitudes corresponding to singly and doubly excited determinants with respect to $|0\rangle$ are determined from the MP2 amplitude equations. To prevent double counting of the correlation contribution associated with electron pairs, we omit the amplitudes of all pair-excited determinants in the MP2 equations. In MP2, the PT amplitudes are determined from a set of uncoupled equations,

[TABLE]

where $|0\rangle$ is the AP1roG reference determinant (cf. eq. (16)). The second-order energy corrections then accounts for correlation contributions beyond electron-pairs and can be determined from eq. (22), which contains both single and double excitations with respect to $|0\rangle$. We will label this methods as AP1roG-MP2.

II.3 Linearized Coupled Cluster Corrections

Besides PT-type methods, dynamic correlation effects can be built in the AP1roG wavefunction a posteriori using an exponential coupled cluster ansatz, Boguslawski and Ayers (2015)

[TABLE]

where $\hat{T}=\sum_{\nu}t_{\nu}\hat{\tau}_{\nu}$ is a general cluster operator. The corresponding time-independent Schrödinger equation then reads

[TABLE]

Instead of considering the full expansion of the cluster operator in the exponential function, we can include terms only linear in $\hat{T}$. If we truncate the Baker–Campbell–Hausdorff expansion after the second term,

[TABLE]

we obtain the linearized coupled cluster Schrödinger equation

[TABLE]

The cluster amplitudes $t_{\nu}$ are determined by solving a linear set of coupled equations obtained from multiplying the above equation from left by determinants of the projection manifold $\langle\nu|$,

[TABLE]

The energy can be calculated by projecting against the reference determinant of $|{\rm AP1roG}\rangle$,

[TABLE]

Specifically, in the AP1roG-LCCD approach, the cluster operator is restricted to double excitations with respect to $|0\rangle$ as defined in eq. (10), while the AP1roG-LCCSD method also single excitations are included and hence $\hat{T}=\hat{T}_{1}+\hat{T}_{2}^{\prime}$ (cf. eqs.(10) and (11)). For more information concerning the choice of the cluster operator and the projection manifold, we refer the reader to ref. 51. Note that the computational scaling of AP1roG-LCCD/LCCSD is similar to the computational scaling of conventional LCCD/LCCSD methods, that is, $\mathcal{O}(o^{2}v^{4})$.

III Computational details

III.1 Geometry optimization and single point calculations

All molecular structures used in the calculations of reaction energies were fully optimized in the  TURBOMOLE7.0 software package Furche et al. (2014), employing the B3LYP exchange–correlation functional Becke (1993); Stephens et al. (1994) and a polarized valence triple-$\zeta$ basis set Eichkorn et al. (1997); Schafer et al. (1994). The resulting geometries in xyz format are included in the Supporting Information.

The single-point MP2, BCCD, and CCSD(T) calculations have been carried out in the Molpro2012.1 software package Werner et al. (2012, 2012), while the NWChem software suite (in combination with the tensor contraction engine Hirata (2003, 2004); Hirata et al. (2004); Kowalski et al. (2011)) was used for LCCD, LCCSD, and CR-CCSD(T) calculations. All AP1roG calculations including a posteriori corrections of the dynamic correlation energy were performed in our in-house quantum chemistry code. The single-point energies of all investigated reactions were evaluated using the cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets of Dunning Dunning (1989) and extrapolated to the basis set limit. Two sets of calculations were performed, one with the canonical Hartree–Fock orbitals and another with the AP1roG optimized orbital basis.

The dissociation curve of the C2 (ground state) molecule was determined for Dunning’s aug-cc-pVTZ basis set, while for N2 (ground state) and BN (first excited state) the cc-pVTZ basis set was employed. All the dynamic energy corrections on top of AP1roG were calculated in the AP1roG optimized orbital basis. The points of the potential energy curves of all investigated diatomics were used for a subsequent generalized Morse function Coxon (1992) fit to obtain the equilibrium bond lengths ($r_{\rm e}$) and potential energy depths ($D_{\rm e}$). The harmonic vibrational frequencies ($\omega_{\rm e}$) were calculated numerically using the five-point finite difference stencil Abramowitz and Stegun (1970).

III.2 Extrapolation to the basis set limit

The basis set limit of the Hartree–Fock energy was obtained by fitting an exponential function of the form Helgaker et al. (1997)

[TABLE]

to the Hartree–Fock energies obtained in the cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets. In the above equation $X$ indicates the cardinal number of the basis set (2 for D, 3 for T, etc.). For all correlation calculations, the basis set limit of the correlation energy was obtained by a two-point fit using the fit function

[TABLE]

as suggested in refs. Helgaker et al. (1997); Halkier et al. (1998). In the above equation, $E^{\mathrm{corr}}({X})$ indicates the correlation energy of a given method defined as $E^{\mathrm{corr}}({X})=E^{\mathrm{tot}}({X})-E^{\mathrm{SCF}}({X})$. In the case of dynamic correction on top of AP1roG, the correlation energy is thus a sum of the correlation energy of AP1roG and of the dynamic correlation model. Note, however, that extrapolating the sum of correlation energies or each contribution separately yields similar basis set limits of the (total) correlation energy $E^{\mathrm{corr}}_{\infty}$. For all correlation calculations, only the correlation energies of the cc-pVTZ ($X=3$) and cc-pVQZ ($X=4$) basis sets were employed in the fitting procedure.

IV Results and discussion

The previously discussed dynamic correlation models are benchmarked against spectroscopic constants of multiply bonded diatomics and thermochemical data for compounds containing main-group elements. Note that the former test systems are dominated by static/nondynamic electron correlation effects, while dynamic electron correlation becomes important for the latter. Both test sets allow us to assess the accuracy of the proposed dynamic corrections for strongly- and weakly correlated molecular systems.

IV.1 Molecular systems dominated by static/nondynamic correlation

We will first briefly discuss the performance and trends of the proposed PT corrections in modeling dissociation processes of three diatomic molecules, C2, N2, and BN. Their electronic structures have been extensively discussed in the literature and we refer the interested reader to, for instances, refs. Peterson et al. (1993); Peterson (1995); Abrams and Sherrill (2004); Boguslawski et al. (2013); Xu and Dunning (2014) for more details. Table 2 summarizes the spectroscopic constants for the C2, N2, and BN molecules. The corresponding potential energy surfaces are shown in the Supplementary Information. For C2 and N2, MRCI-SD data by Peterson et al. Peterson et al. (1993) is taken as reference, while MRCI-SD+Q results are used as reference for the BN molecule Peterson (1995).

In general, PT2SDd and PT2SDo provide equilibrium bond lengths and vibrational frequencies that agree well with the corresponding reference values for all investigated diatomics. Furthermore, addition of single excitations in the excitation operator $\hat{T}$ does not increase the accuracy of PT2SDd(d) or PT2SDo(d). In contrast to PT2SDd and PT2SDo, the performance of PT models using a multi-determinant dual state (PT2MDd, PT2MDo, and PT2b) is less satisfying. Specifically, these methods provide equilibrium bond distances and vibrational frequencies that deviate most from MRCI-SD/MRCI-SD+Q reference data. Similar to PT2SDd and PT2SDo, addition of single excitations worsens spectroscopic constants and hence the excitation operator should be restricted to double excitations only. Furthermore, including pair excitations in $\hat{T}$ (PT2b-type methods) does not significantly affect spectroscopic constants and both models (with and without pair excitations) yield similar values for $r_{e}$ and $\omega_{e}$.

In contrast to equilibrium bond distances and vibrational frequencies, the accurate prediction of dissociation energies $D_{e}$ is more challenging. In general, none of the proposed PT models can reliably predict potential energy well depths, which differ by approximately 15 to 30 kcal/mol from MRCI-SD/MRCI-SD+Q reference data. Note that most PT corrections (PT2MD and PT2b) diverge in the vicinity of dissociation and hence only an estimated dissociation energy is given in Table 2 (indicated by the $*$ in the Table).

The performance of the LCCD/LCCSD correction on top of AP1roG is more robust in predicting spectroscopic constants for all investigated diatomics. Although AP1roG-LCCSD yields $r_{e}$ and $\omega_{e}$ that deviate more from MRCI-SD/MRCI-SD+Q reference values compared to PT2SDd and PT2SDo results, an LCC correction allows us to reliably model dissociation energies for C2 and N2, while the errors in $D_{e}$ increase for BN (differences are 2 kcal/mol for C2, 6 kcal/mol for N2, and 24 kcal/mol for BN). Note that for the N2 molecule the coupled cluster doubles amplitudes for excitations within the valence shell (vs) have been kept frozen in the vicinity of dissociation ($r>2.0$ Å with $t_{r>2.0}^{\textrm{vs}}=t_{r=2.0}^{\rm vs}$) to prevent divergencies in the dissociation limit and to obtain a smooth potential energy surface. This is indicate by fLCCD/fLCCSD in Table 2. The corresponding AP1roG-fLCCSD spectroscopic constants agree well with CASSCF results by Peterson et al. Peterson et al. (1993) and are in general closest to MRCI-SD reference values.

IV.2 Molecular systems dominated by dynamic correlation

Our second test case includes 15 reactions containing main-group elements summarized in Table 3. For all studied systems, the reaction energies obtained from CR-CCSD(T) calculations are taken as reference values. Furthermore, we will focus on two different molecular orbital basis sets used in our AP1roG calculations and in all a posteriori corrections for dynamic correlation with an AP1roG reference function; these are optimized natural AP1roG orbitals and canonical Hartree–Fock orbitals.

IV.2.1 Optimized natural AP1roG orbitals

Figure 1 shows the reaction energies and root mean square errors (RMSEs) with respect to CR-CCSD(T) reference data (grey bars) for all investigated methods in the optimized AP1roG basis (see Table 4 for the definition of all error measures). Compared to Hartree–Fock data, AP1roG predicts reaction energies that differ more from the CR-CCSD(T) reference. In general, AP1roG underestimates reaction energies and yields an RMSE of approximately 15 kcal/mol (see also Table 4 for cc-pVQZ). This discrepancy can be attributed to electron correlation effects beyond electron pairs that are not accounted for in the AP1roG method. Inclusion of open-shell determinants in the electronic wavefunction improves reaction energies. However, the performance of the PT models strongly depends on the choice of the zero-order Hamilton and the perturbation. Specifically, all PT methods with a diagonal one-electron Hamiltonian (PT2SDd and PT2MDd) provide slightly improved reaction energies with an RME of around 13 kcal/mol (see also Table 4 for cc-pVQZ). Yet, these methods are strongly basis-set dependent if optimized AP1roG orbitals are used as orbital basis: the RMSE gradually increases with the basis set size (see Figure 1). An a posteriori MP2 correction in the AP1roG basis predicts the largest RMSE of approximately 16 kcal/mol and therefore does not allow us to accurately model reaction energies. Choosing an off-diagonal one-electron zero-order Hamiltonian (PT2SDo and PT2MDo) in the PT model further reduces the error with respect to the CCSD(T) reference with an RMSE of approximately 4-5 kcal/mol. The error in reaction energies predicted by PT corrections can be minimized if the perturbation Hamiltonian of PT2MDo (see Table 1) is replaced by the full quantum chemical Hamiltonian $\hat{H}_{N}$ as in all PT2b-type methods ($\mathrm{RMSE}\approx 3.0$ kcal/mol), where the accuracy is similar to MP2 theory (see Table 4). Note that inclusion of single excitations in the excitation operator ($\hat{T}=\hat{T}_{1}+\hat{T}_{2}^{\prime}$) does not affect the accuracy of the PT corrections in predicting reaction energies (see also Table 4). Only the performance of PT2b is slightly improved.

For all investigated basis sets, the best agreement with CR-CCSD(T) reference data is obtained for the linearized coupled cluster correction including singles and doubles (see Figure 1) with an RMSE of 1.5 kcal/mol and a mean absolute error of 1.2 kcal/mol. Note that the errors of AP1roG-LCCSD is similar to those of conventional CC methods like BCC, CCSD, CCSD(T), and LCCD. Only LCCSD performs worse with an RMSE of 2.5 kcal/mol and a mean absolute error of 1.8 kcal/mol. We should emphasize that, unlike the LCC method with an Hartree–Fock reference function, single excitations in the LCC corrections on top of AP1roG are particularly important and allow us to approach chemical accuracy (approximately 1 kcal/mol) when modeling thermochemistry of main-group compounds.

IV.2.2 Canonical Hartree–Fock orbitals

If canonical Hartree–Fock orbitals are chosen as orbital basis, we can use the diagonal structure of the Fock matrix and simplify the PT amplitude equations. Specifically for all PT models, the PT amplitude equations are transformed into a set of uncoupled linear equations as $\hat{F}_{N}^{\mathrm{o}}=0$ (cf. Table 1). Furthermore, single excitations do not contribute (either directly in the energy expression or indirectly through coupling to the doubles manifold) if the AP1roG reference determinant is taken as dual $\langle\tilde{\Psi}|$. Thus, the PT2SDd and PT2SDo methods become equivalent (with and without single excitations) and only the PT2SDd(d) results are shown in Table 5 and Figure 2. Similarly, the PT2MDd and PT2MDo models are equivalent to each other. In contrast to PT2SD-type approaches, single excitations do contribute directly to the second-order energy correction as their contribution is determined from eq. (26), which also contains one-electron terms. Finally, in PT2b-type methods, the coupling to pair-excited determinants vanishes in the PT amplitude equations as they are coupled through off-diagonal Fock matrix elements. Furthermore, the perturbation Hamiltonian $\hat{H}_{N}^{\prime}$ of PT2b reduces to $\hat{V}_{N}^{\prime}$ and both the amplitude equations and second-order energy expression are equivalent to those in PT2MD-type methods. If $\langle\mathrm{AP1roG}|$ is chosen as dual, all PT methods coincide in the canonical Hartree–Fock basis. Therefore, Table 5 and Figure 2 show only the PT2MDd(d) and PT2MDd(sd) results. The canonical Hartree–Fock basis thus allows us to directly assess how the choice of the dual state affects the performance of the PT models.

In contrast to the optimized natural AP1roG basis, all PT models yield similar reaction energies with a mean absolute error of approximately 3 kcal/mol and a RMSE of about 3.5 kcal/mol, irrespective of the excitation operator $\hat{T}$ (see Table 5 and Figure 2). Only the combination of MP2 (for open-shell configurations) and AP1roG (for closed-shell configurations) performs worse resulting in an RMSE of 4.4 kcal/mol. As observed above, the LCC correction outperforms all PT models if both single and double excitations are included in the cluster operator (cf. Table 5 and Figure 2), where the RMSE reduces to 1.0 kcal/mol. Furthermore, the accuracy of the PT methods does not deteriorate for increasing sizes of the basis set (see Figure 2).

Finally, we should note that the performance of all dynamic correlation corrections depends on the molecular orbital basis. The errors with respect to CR-CCSD(T) reference data significantly decrease for PT methods with a single determinant as dual $\langle 0|$ when canonical Hartree–Fock orbitals are used as molecular orbital basis. In cases of $\langle\mathrm{AP1roG}|$, the dependence on the molecular orbital basis is less severe. While reaction energies obtained by PT2MDd/PT2MDo in the canonical Hartree–Fock basis slightly improve reducing the error by a factor of 2, the mean absolute error and RMSE of PT2b remain almost unchanged compared to the natural AP1roG basis (cf. a RMSE of 3.0 kcal/mol for AP1roG orbitals and of 3.6 kcal/mol for canonical Hartree–Fock orbitals). Similar to PT2MD-type methods, the LCC corrections are less sensitive to the choice of the molecular orbital basis and yield reaction energies with slightly smaller error measures (cf. a RMSE of 1.5 kcal/mol for AP1roG orbitals and LCCSD and of 1.0 kcal/mol for canonical Hartree–Fock orbitals and LCCSD). We should emphasize that AP1roG-LCCSD in the canonical Hartree–Fock basis outperforms all conventional CC methods investigated in this work (BCC, CCSD, CCSD(T), LCCD, and LCCSD) and provides reaction energies that are within chemical accuracy compared to CR-CCSD(T) reference data.

V Conclusions

Wavefunctions restricted to electron-pair states allow us to reliably model static/nondynamic electron correlation effects. Boguslawski et al. (2016) However, in order to predict spectroscopic constants and thermochemistry within chemical accuracy, we need to include dynamic correlation effects that go beyond the simple electron-pair model. This can be achieved a posteriori using PT approaches, Limacher et al. (2014b); Limacher (2015) coupled cluster corrections, Henderson et al. (2014b); Boguslawski and Ayers (2015) or DFT-type methods. Garza et al. (2015a, b) In this work, we have extended the previously presented PT models with an AP1roG reference function and benchmarked those models against spectroscopic constants for multiply bonded diatomics and thermochemical data extrapolated to the basis set limit. Most importantly, combining AP1roG with the investigated corrections allows us to reliably model molecular systems dominated by both static/nondynamic and dynamic electron correlation.

Specifically, our PT extensions combine a diagonal and off-diagonal one-electron zero-order Hamiltonian, a single-determinant and multi-determinant dual state, and a projection manifold restricted to double as well as single and double excitations. In general, the performance of all PT methods can be divided in three different groups: (i) those with a diagonal zero-order Hamiltonian (PT2SDd/PT2MDd), (ii) those with an off-diagonal zero-order Hamiltonian (PT2SDo/PT2MDo), and (iii) those with an off-diagonal zero-order Hamiltonian and the full quantum-chemical Hamiltonian as perturbation operator (PT2b-type methods). For the dissociation of multiply bonded diatomics, the PT corrections using a single-determinant dual state outperform all other investigated PT models. In particular, the PT2SDd and PT2SDo methods (employing a diagonal and an off-diagonal zero-order Hamiltonian) provide accurate equilibrium bond lengths and vibrational frequencies compared to MRCI-SD/MRCI-SD+Q reference data, followed by AP1roG-LCCSD that is similar in accuracy, though equilibrium bond distances deviate more from reference values. The reliable prediction of dissociation energies, however, remains challenging for all proposed PT models with errors between 15 to 30 kcal/mol. While PT2SDd and PT2SDo provide smooth potential energy surfaces, all remaining PT methods diverge in the vicinity of dissociation. In order to model the dissociation of multiply bonded diatomics, an LCCSD correction has to be applied.

In case of reaction energies, the accuracy of the PT corrections with respect to CR-CCSD(T) reference data increases when going from PT methods (i) to (iii) reducing the RMSE from approximately 14 kcal/mol in PT2SDd/PT2MDd, to about 5 kcal/mol in PT2SDo/PT2MDo, to 3 kcal/mol in PT2b-type methods. The choice of the dual state and the inclusion of single excitations in the excitation operator do not significantly affect the accuracy of the PT methods (mean error, root mean square error, mean absolute error, maximum absolute error). Furthermore, excluding pair-excited determinants from the projection manifold in PT2b-type methods improves the accuracy of PT2b only marginally. Since pair-excitations are already described in the AP1roG reference function, it might, however, be advantageous to exclude pair excitations from the excitation operator $\hat{T}$ and hence eliminate the coupling to pair excitations modeled in the AP1roG reference function and pair excitations of the PT method, which both couple to the remaining PT amplitudes in the PT amplitude equations.

If the optimized natural AP1roG orbitals are replaced by canonical Hartree–Fock orbitals, only two distinct PT models persist, namely, PT2SDd (with double excitations) and PT2MDd (with double as well as single and double excitations). In contrast to the natural AP1roG orbitals basis, all PT methods yield similar error measures in the canonical Hartree–Fock basis with a standard error of 3.6 kcal/mol. Therefore, the optimization of the molecular orbital basis and the AP1roG reference determinant might be unnecessary if the molecular system is dominated by dynamic correlation and molecular properties around the equilibrium geometry are considered, provided dynamic correlation effects are accounted for in the AP1roG model. If the optimal natural AP1roG orbitals are used in calculations, PT2b-type methods result in the smallest error measures (around 3.0 kcal/mol) and thus outperform all other PT models. If the orbital optimization step is omitted, PT2SDd/PT2SDo provide the smallest errors that are similar to PT2b-type methods in the optimized AP1roG basis (RMSE around 3.3 kcal/mol).

Finally, a linearized coupled cluster correction with an AP1roG reference function, as presented in ref. 51, predicts reaction energies that deviate least from CR-CCSD(T) reference data reducing the RMSE to 1.5 kcal/mol. To minimize the error in AP1roG-LCC, single excitations are indispensable and have to be included in the cluster operator, both using optimized AP1roG orbitals and canonical Hartree–Fock orbitals. Most importantly, the AP1roG orbital basis does not need to be optimized if chemical accuracy (approximately 1 kcal/mol) is desired for predicting equilibrium properties of weakly-correlated systems. To conclude, AP1roG-LCCSD provides the most accurate reaction energies with respect to CR-CCSD(T) reference data, outperforming all investigated PT models as well as conventional electronic structure methods like MP2, BCC, CCSD, CCSD(T), LCCD, and LCCSD. In order to describe equilibrium properties of the multiply bonded diatomics C2, N2, and BN, PT2SDd and PT2SDo (with double excitations only) outperform AP1roG-LCCSD, fail, however, in predicting accurate dissociation energies. Furthermore, PT2SDd/PT2SDo in the canonical Hartree–Fock basis provides the smallest errors among all investigated PT corrections (similar to MP2) and allows us to cheaply model the thermochemistry of main group elements ($\mathcal{O}(o^{2}v^{2})$). For strongly-correlated systems, however, the molecular orbital basis needs to be optimized before a PT2SDd or PT2SDo correction is applied. Rotating the orbital basis increases the computational cost due to the four-index transformation and some additional prefactor of the orbital optimization.

VI Acknowledgement

K.B. acknowledges financial support from a SONATA BIS grant of the National Science Centre, Poland (no. 2015/18/E/ST4/00584), a Marie-Skłodowska-Curie Individual Fellowship project no. 702635–PCCDX, and the Foundation for Polish Science for a START2016 stipend. P.T. thanks an OPUS grant of the National Science Centre, Poland, no. DEC-2013/11/B/ST4/00771 and an OPUS grant of the National Science Centre, Poland, UMO-2015/19/B/ST4/02707.

Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (http://wcss.pl), grant no. 411 and grant no. 412.