Electron correlation methods based on pair natural orbitals (PNOs) have gained an increasing degree of interest in recent years, as they permit energy calculations to be performed on systems containing up to many hundred atoms, while maintaining chemical accuracy for reaction energies. We present an approach for taking exact analytical first derivatives of the energy contributions in the simplest method of the family of Domain-based Local Pair Natural Orbital (DLPNO) methods, closed-shell DLPNO-MP2. The Lagrangian function contains constraints to account for the relaxation of PNOs. RI-MP2 reference geometries are reproduced accurately, as exemplified for four systems with a substantial degree of nonbonding interactions. By the example of electric field gradients, we demonstrate that omitting PNO-specific constraints can lead to dramatic errors for orbital-relaxed properties.

Local correlation methods based on pair natural orbitals (PNOs) permit accurate calculations at many-body perturbation theory or coupled cluster level to be performed routinely for systems containing hundreds of atoms.1–12 Our approach, the Domain Local Pair Natural Orbital (DLPNO) family of methods,4 combines PNOs with the domain concept for projected atomic orbitals (PAOs), which were introduced by Pulay and Saebø.13–15 The success of PNO-based methods in calculating single-point energies suggests that analytical derivatives could be equally useful for other molecular properties.

Werner, Schütz, and co-workers implemented an analytical gradient and computation of NMR shifts for PAO-based local MP2,16–21 as well as an analytical gradient for the PAO-based local CC2 method.22,23 More recently, Frank and Hättig presented a gradient implementation for PNO-based MP2;24 however, the derivatives of the PNO-MP2 energy are approximate, as the relaxation of orbital-specific virtuals (OSVs) and PNOs was not taken into account.

In this communication, we report a Lagrangian approach25–27 for exact first derivatives of the closed-shell DLPNO-MP2 energy7 and demonstrate its performance for electric field gradients and for geometry optimizations. We investigate the importance of the individual constraints in the Lagrangian in order to find out whether more approximate, less complex schemes might hold promise. A detailed description of the working equations and their implementation, accompanied by more extensive calculations, will be the subject of a forthcoming article.

The DLPNO-MP2 correlation energy EDLPNO is obtained by minimizing the Hylleraas functional in the PNO basis. Furthermore, the total correlation energy EC is corrected with a contribution ΔEPNO for PNO space truncation, and with an estimate ΔEPre for the energy of the screened-out orbital pairs,

(1)

We express the Hylleraas functional for EDLPNO using the virtual unrelaxed pair density Dij and the occupied unrelaxed difference density Dij,

(2)
(3)
(4)

Pair natural orbitals ãij are a set of orthonormal functions representing a truncated virtual space for the pair of localized occupied orbitals ij. We drop the subscript for the sake of simplicity. In Eq. (4), PNOs ã and d̃ span the virtual space of pair ik, whereas b̃ and c̃ belong to kj. The contravariant amplitudes T̃ij are defined as 1+δij14Tij2Tij.

In order to construct pair natural orbitals, we introduce the orbital-unrelaxed semicanonical pair density with the virtual component Dij,

(5)

Semicanonical amplitudes Tij are defined in the basis of non-redundant pseudo-canonical PAOs μ̃,ν̃ belonging to the domain of pair ij,

(6)

We label the full set of eigenvectors of the matrix Dij as ã (they span the same space as the non-redundant PAOs of the pair). PNOs ã with an eigenvalue larger than the threshold TCutPNO are kept for the energy calculation as a compact representation of the virtual space, while the remaining eigenvectors, referred to as the complementary PNOs ã, are discarded. To estimate the error incurred from PNO truncation, the semicanonical amplitudes are projected onto the PNO basis, and the untruncated and truncated pair energies subtracted,

(7)

We derive the gradient of the DLPNO-MP2 energy using the Lagrangian function in Eq. (8). The constraints expressed through s, Dij, and Rij, and associated with the respective multipliers z, vij, and wij, will be explained after introducing the Lagrangian,

(8)

where EHF is the closed-shell Hartree-Fock energy and EC is the correlation energy defined in Eq. (1). In this preliminary work, however, we omit the dipole energy contribution of the screened-out pairs ΔEPre. Nonetheless, we emphasize that the derivatives of all the other energy components are taken exactly, while the energy of the screened-out pairs is considerably smaller than the errors introduced by PNO truncation and domain selection.

Perturbed molecular orbital coefficients c are related to the occupied and virtual unperturbed orbitals c(0) through an exponential parameterization with the antisymmetric matrix κ=κ (as described by Helgaker and Almlöf28),

(9)

where S is the overlap matrix of the unperturbed orbitals c(0) in the basis of perturbed atomic basis functions (AOs). Separation of the occupied and virtual molecular orbitals (i and a, respectively) is ensured through the Brillouin condition in the Lagrangian

(10)

Because of the local approximations, the correlation energy is not invariant with respect to unitary transformations among occupied orbitals, necessitating a constraint specific to the localization method. El Azhary and co-workers showed for Pipek-Mezey orbitals29 that the localization function needs to be stationary with respect to unitary transformations among the occupied orbitals.18 We apply the same argument to the Foster-Boys orbitals30,31 used in this work, which thus need to satisfy the following equation:

(11)

After the initial publication of the DLPNO-MP2 method, it was found out that a tighter PNO threshold is needed for the core region.32 In the current implementation, TCutPNO is scaled with a factor of 10−2 for all pairs with at least one core orbital. Results from the first article on DLPNO-MP2 remain unaffected, as the frozen-core approximation was used throughout. The core and valence orbitals are localized independently to avoid mixing; therefore, Eq. (11) applies separately to both groups of orbitals in the Lagrangian. In addition, the core-valence blocks of the Fock matrix are constrained to zero.

Derivatives of PNOs with respect to the orbital parameters κpq cannot be taken straightforwardly—unlike for redundant PAOs, which can be expressed as a simple function of molecular orbitals. This necessitates an additional set of parameters that we introduce in analogy to the MOs,

(12)

The coefficients dμ̃ãij represent the full set of PNOs ã in the basis of redundant PAOs μ̃ of the pair domain {ij}. Therefore, the columns of dij0 span the linearly independent subspace of the PAOs. A finite perturbation of the system can lead to a domain reassignment or change the rank of the PAO coefficient subset. However, an infinitesimal perturbation cannot induce a finite change, provided that the one-particle density matrix of the reference determinant is a continuous function of the external perturbation. For this reason, dij0 is a non-orthogonal representation of the linearly independent subspace spanned by the infinitesimally perturbed PAOs. This is followed by symmetric orthogonalization of the full PNO set and a unitary transformation, where θij is an antisymmetric, pair-specific matrix.

Pair natural orbitals are selected by their eigenvalues of the semicanonical pair density matrix. However, a constraint requiring the perturbed PNOs to be eigenvectors of the semicanonical density matrix is not appropriate. Since there is no unique choice of eigenvectors for degenerate eigenvalues, derivatives of repeated eigenvalues are not generally well defined. A related problem affects the derivatives of canonical electron correlation methods: enforcing canonicity of the perturbed orbitals leads to singularities in the coupled-perturbed self-consistent field (CP-SCF) equations.33 Helgaker and co-workers showed that the Brillouin condition [Eq. (10)] should be used instead of a stronger constraint Fpq=δpqεp for methods that are invariant to occupied-occupied or virtual-virtual orbital rotations.26 

The DLPNO-MP2 energy is invariant to rotations among the truncated PNOs so that it is sufficient to prevent mixing between the spaces of retained PNOs ã and of the complementary PNOs b̃ that were discarded. We choose the third constraint in the Lagrangian analogously to the Brillouin condition

(13)

The eigenvectors of this block-diagonal matrix are linear combinations of either only the PNOs ã with eigenvalues above TCutPNO or only the complementary PNOs b̃ with eigenvalues below TCutPNO.

Equation (13) introduces the semicanonical amplitudes Tij as an additional set of variables, distinct from the PNO basis amplitudes Tij in the Hylleraas functional [Eq. (2)]. Therefore, the fourth constraint in the Lagrangian is the generalization of Eq. (6) to a non-canonical virtual basis,

(14)

We still label the amplitudes Tij as semicanonical because off-diagonal Fock matrix elements are neglected only for the occupied part. Taking Tij as parameters of the Lagrangian leads to redundancies in the Z-vector equations, which is resolved by a reparameterization,

(15)

To summarize, the independent parameters in the Lagrangian are the MO rotation matrices κ, the PNO rotation matrices θij, the PNO basis amplitudes Tij, and the amplitudes Tij defined in Eq. (15). As the semicanonical amplitudes are fully constrained, the PNO truncation correction ΔEPNO, as specified in Eq. (7), is taken as part of the energy functional in the Lagrangian. The multipliers for the PNO and semicanonical conditions are determined from stationarity of the Lagrangian with the parameters θij and amplitudes Tij, respectively, resulting in independent sets of equations for each orbital pair. Making the Lagrangian stationary with respect to orbital rotations κ leads to Z-vector equations for the Brillouin and localization constraints.

We have implemented the DLPNO-MP2 gradient in a development version of the ORCA34 package. Besides the missing energy contribution for the screened-out pairs, the implementation does not yet contain AO-based coefficient truncation for the additional RI integrals and their derivatives. All calculations in this work were performed without the frozen core approximation. Geometric derivatives of dipole integrals, needed for the gradient of the Foster-Boys localization condition, are evaluated numerically on a grid. Except for these minor limitations, we implement the complete derivatives of the Hylleraas functional EDLPNO and of the PNO correction ΔEPNO. The localized orbitals need to satisfy Eq. (11) to a sufficient degree of precision, but we found that tight convergence is often difficult to achieve with the default procedure of 2 × 2 rotations. To circumvent this, we use a localization implementation based on rational function optimization,35 in combination with Davidson diagonalization36,37 for large systems.

The complexity introduced by the four conditions in the Lagrangian makes it desirable to establish if all of them are indeed necessary to obtain accurate results. We performed calculations using four types of setups: first, the respective exact derivatives of the energy containing the PNO energy correction ΔEPNO (1) and of the uncorrected DLPNO-MP2 energy (2). Second, we consider approximate derivatives of the energy without PNO correction: with the Brillouin and localization conditions, but without the PNO and semicanonical constraints in the Lagrangian (3), and with the Brillouin condition only (4). For the two types of approximate derivatives (3) and (4), the PNO parameterization in Eq. (12) cannot be used, as orthogonalizing the untruncated PNO vectors leads to uncontrollable errors without a constraint on the rotations θij. In those cases, we span the PNO space by the truncated set of the frozen coefficients dμ̃ãij0 instead (an additional rotation is unnecessary because of the energy invariance)

(16)

In preliminary calculations, we found that the influence of the PNO-specific conditions on electrostatic properties is often relatively small, but sometimes their omission can lead to drastic errors, as we will demonstrate for electric field gradients (EFGs) at nuclear positions. Figure 1 shows the MP2 contribution to the EFG principal components Vzz, determined as the difference between Vzz values computed at MP2 and Hartree-Fock levels, for chlorothiazide and a taurine zwitterion. Only the PNO selection threshold TCutPNO is varied while adjusting the core-level threshold to 10−2 × TCutPNO, whereas all other thresholds are set to zero to eliminate their associated approximations.

FIG. 1.

MP2 contribution to the largest principal component of the electric field gradient as a function of TCutPNO. All other thresholds are set to zero. (a) Cl nucleus of chlorothiazide, (b) S nucleus of taurine. Basis: cc-pwCVTZ.38 

FIG. 1.

MP2 contribution to the largest principal component of the electric field gradient as a function of TCutPNO. All other thresholds are set to zero. (a) Cl nucleus of chlorothiazide, (b) S nucleus of taurine. Basis: cc-pwCVTZ.38 

Close modal

Electric field gradients obtained from the fully relaxed densities show good convergence with the PNO truncation threshold towards the RI-MP2 result. The errors are reduced upon including the energy correction for PNO truncation ΔEPNO. On the other hand, omitting the constraints for the PNOs and for the semicanonical amplitudes introduces unacceptably large errors and erratic, non-monotonous convergence. With the threshold 10−8, used in the default settings, the errors in the correlation contribution to the EFGs amount to 20% in the approximate derivatives schemes, but only to ca. 1% with the full set of constraints. We demonstrate in the supplementary material that this behavior is not an artifact of core-level PNO truncation.

To examine the errors in a realistic context, Fig. 2 compares DLPNO-MP2 results for nuclear quadrupole coupling constants calculated with existing accuracy settings to the canonical reference values and to experimental results deduced from solid-state NMR spectra of the two the compounds.39,40 Calculations for the 35Cl nucleus of chlorothiazide were performed on an isolated molecule in the gas phase using the cc-pwCVQZ basis set.38 At the 33S nucleus in taurine, which exists entirely in its zwitterionic form in the solid state, there is a substantial intermolecular contribution to the EFG. Therefore, we created a model of four molecules embedded in point charges. Deviations of the DLPNO-MP2 results from the RI-MP2 reference are small: 0.05 MHz with NormalPNO and 0.03 MHz with TightPNO for taurine, and even less for chlorothiazide. On the other hand, omitting PNO-specific constraints leads to excessive errors relative to RI-MP2 with any of the chosen accuracy settings.

FIG. 2.

Nuclear quadrupole coupling constants for (a) chlorothiazide and (b) taurine. The horizontal dotted lines mark results for Hartree-Fock and RI-MP2, and experimental values from solid-state NMR experiments.39,40 Grouping of DLPNO-MP2 results: (1) full derivative with ΔEPNO, (2) full derivative without ΔEPNO, (3) no PNO-specific constraints, (4) no PNO and localization constraints.

FIG. 2.

Nuclear quadrupole coupling constants for (a) chlorothiazide and (b) taurine. The horizontal dotted lines mark results for Hartree-Fock and RI-MP2, and experimental values from solid-state NMR experiments.39,40 Grouping of DLPNO-MP2 results: (1) full derivative with ΔEPNO, (2) full derivative without ΔEPNO, (3) no PNO-specific constraints, (4) no PNO and localization constraints.

Close modal

In conclusion, orbital-relaxed properties are only reliable if the relaxation of PNOs is taken into account. An alternative route to first-order properties exists within the DLPNO framework through the orbital-unrelaxed density matrix for coupled cluster with single and double excitations (DLPNO-CCSD).41 A finding, shared with several other authors, is that wave function truncation makes unrelaxed molecular properties prone to substantial errors.41–44 Our results do not lead us to conclude that PNO relaxation alone would have a significant effect on orbital-unrelaxed coupled cluster densities: to the contrary, numerically computed unrelaxed CCSD dipole moments showed that a deeper understanding of contributions beyond PNO relaxation might be necessary to improve accuracy.41 

We have performed geometry optimizations for four systems with a substantial influence of nonbonding interactions (Fig. 3): 1,5-dibromo-4,8-dichloronaphthalene and RESVAN were taken from the LB12 set of Grimme and co-workers,46 whereas the structures of 1,1′-bi-2-naphthol (BINOL) and the adenine-thymine dimer were pre-optimized with the HF-3c method.47 Most optimizations converged with the TightOpt settings of ORCA; only for the adenine-thymine dimer with LoosePNO thresholds, the less rigorous default settings were used. The RI-MP2 geometries are reproduced accurately at NormalPNO and TightPNO levels, while errors in covalent bond lengths are usually significantly below 0.1 pm. Geometry optimizations performed with approximate derivatives are discussed in the supplementary material.

FIG. 3.

Geometries optimized with DLPNO-MP2 and RI-MP2 in comparison. The calculations were performed with the def2-TZVP basis.45 

FIG. 3.

Geometries optimized with DLPNO-MP2 and RI-MP2 in comparison. The calculations were performed with the def2-TZVP basis.45 

Close modal

We show for linear alkane chains in Fig. 4 that the crossover with the canonical RI-MP2 implementation occurs around C15H32. A DLPNO-MP2 gradient calculation remains feasible throughout all system sizes for which an SCF calculation can be performed, even if it is accelerated with RIJCOSX.48 Timing results are provided in the supplementary material.

FIG. 4.

Scaling of the total cost of single-point SCF and MP2 gradient calculations for linear alkane chains with NormalPNO settings and the RIJCOSX approximation. Basis set: def2-TZVP.45 The sum of the times needed for the SCF and CP-SCF steps only is shown for comparison.

FIG. 4.

Scaling of the total cost of single-point SCF and MP2 gradient calculations for linear alkane chains with NormalPNO settings and the RIJCOSX approximation. Basis set: def2-TZVP.45 The sum of the times needed for the SCF and CP-SCF steps only is shown for comparison.

Close modal

See supplementary material for detailed description of the computational procedures with tables of numerical results, further analysis with additional calculations, molecular geometries, and geometry optimization trajectories.

The authors gratefully acknowledge financial support of this work by the Max-Planck-Society and Fonds der Chemischen Industrie. We thank Róbert Izsák and Yang Guo for stimulating discussions.

1.
W.
Meyer
,
Int. J. Quantum Chem.
5
,
341
(
1971
).
2.
W.
Meyer
,
J. Chem. Phys.
58
,
1017
(
1973
).
3.
F.
Neese
,
F.
Wennmohs
, and
A.
Hansen
,
J. Chem. Phys.
130
,
114108
(
2009
).
4.
C.
Riplinger
and
F.
Neese
,
J. Chem. Phys.
138
,
034106
(
2013
).
5.
C.
Riplinger
,
B.
Sandhoefer
,
A.
Hansen
, and
F.
Neese
,
J. Chem. Phys.
139
,
134101
(
2013
).
6.
C.
Riplinger
,
P.
Pinski
,
U.
Becker
,
E. F.
Valeev
, and
F.
Neese
,
J. Chem. Phys.
144
,
024109
(
2016
).
7.
P.
Pinski
,
C.
Riplinger
,
E. F.
Valeev
, and
F.
Neese
,
J. Chem. Phys.
143
,
034108
(
2015
).
8.
H.-J.
Werner
,
G.
Knizia
,
C.
Krause
,
M.
Schwilk
, and
M.
Dornbach
,
J. Chem. Theory Comput.
11
,
484
(
2015
).
9.
M.
Schwilk
,
Q.
Ma
,
C.
Köppl
, and
H.-J.
Werner
,
J. Chem. Theory Comput.
13
,
3650
(
2017
).
10.
G.
Schmitz
,
B.
Helmich
, and
C.
Hättig
,
Mol. Phys.
111
,
2463
(
2013
).
11.
G.
Schmitz
,
C.
Hattig
, and
D. P.
Tew
,
Phys. Chem. Chem. Phys.
16
,
22167
(
2014
).
12.
G.
Schmitz
and
C.
Hättig
,
J. Chem. Phys.
145
,
234107
(
2016
).
13.
14.
S.
Saebø
and
P.
Pulay
,
Chem. Phys. Lett.
113
,
13
(
1985
).
15.
P.
Pulay
and
S.
Saebø
,
Theor. Chim. Acta
69
,
357
(
1986
).
16.
M.
Schütz
,
G.
Hetzer
, and
H. J.
Werner
,
J. Chem. Phys.
111
,
5691
(
1999
).
17.
H. J.
Werner
,
F. R.
Manby
, and
P. J.
Knowles
,
J. Chem. Phys.
118
,
8149
(
2003
).
18.
A.
El Azhary
,
G.
Rauhut
,
P.
Pulay
, and
H.-J.
Werner
,
J. Chem. Phys.
108
,
5185
(
1998
).
19.
M.
Schütz
,
H.-J.
Werner
,
R.
Lindh
, and
F. R.
Manby
,
J. Chem. Phys.
121
,
737
(
2004
).
20.
J.
Gauss
and
H.-J.
Werner
,
Phys. Chem. Chem. Phys.
2
,
2083
(
2000
).
21.
S.
Loibl
and
M.
Schütz
,
J. Chem. Phys.
137
,
084107
(
2012
).
22.
D.
Kats
and
M.
Schütz
,
J. Chem. Phys.
131
,
124117
(
2009
).
23.
K.
Ledermüller
and
M.
Schütz
,
J. Chem. Phys.
140
,
164113
(
2014
).
24.
M. S.
Frank
,
G.
Schmitz
, and
C.
Hättig
,
Mol. Phys.
115
,
343
(
2017
).
25.
P.
Jørgensen
and
T.
Helgaker
,
J. Chem. Phys.
89
,
1560
(
1988
).
26.
T.
Helgaker
,
P.
Jørgensen
, and
N. C.
Handy
,
Theor. Chim. Acta
76
,
227
(
1989
).
27.
J.
Gauss
, in
Modern Methods and Algorithms of Quantum Chemistry
, edited by
J.
Grotendorst
(
John von Neumann Institute for Computing, Forschungszentrum Jülich
,
Germany
,
2000
), p.
541
.
28.
T. U.
Helgaker
and
J.
Almlöf
,
Int. J. Quantum Chem.
26
,
275
(
1984
).
29.
J.
Pipek
and
P. G.
Mezey
,
J. Chem. Phys.
90
,
4916
(
1989
).
30.
J. M.
Foster
and
S. F.
Boys
,
Rev. Mod. Phys.
32
,
300
(
1960
).
31.
S. F.
Boys
, in
Quantum Theory of Atoms, Molecules, and the Solid State
, edited by
P.-O.
Löwdin
(
Academic Press
,
New York
,
1966
), p.
253
.
32.
G.
Bistoni
,
C.
Riplinger
,
Y.
Minenkov
,
L.
Cavallo
,
A. A.
Auer
, and
F.
Neese
,
J. Chem. Theory Comput.
13
,
3220
(
2017
).
33.
N. C.
Handy
,
R. D.
Amos
,
J. F.
Gaw
,
J. E.
Rice
, and
E. D.
Simandiras
,
Chem. Phys. Lett.
120
,
151
(
1985
).
34.
F.
Neese
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1327
(
2018
).
35.
A.
Banerjee
,
N.
Adams
,
J.
Simons
, and
R.
Shepard
,
J. Phys. Chem.
89
,
52
(
1985
).
36.
E. R.
Davidson
,
J. Comput. Phys.
17
,
87
(
1975
).
37.
M.
Crouzeix
,
B.
Philippe
, and
M.
Sadkane
,
SIAM J. Sci. Comput.
15
,
62
(
1994
).
38.
K. A.
Peterson
and
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
117
,
10548
(
2002
).
39.
F. A.
Perras
and
D. L.
Bryce
,
Angew. Chem., Int. Ed.
51
,
4227
(
2012
).
40.
L. A.
O’Dell
and
C. I.
Ratcliffe
,
J. Phys. Chem. A
115
,
747
(
2011
).
41.
D.
Datta
,
S.
Kossmann
, and
F.
Neese
,
J. Chem. Phys.
145
,
114101
(
2016
).
42.
T.
Korona
,
K.
Pflüger
, and
H.-J.
Werner
,
Phys. Chem. Chem. Phys.
6
,
2059
(
2004
).
43.
H. R.
McAlexander
and
T. D.
Crawford
,
J. Chem. Theory Comput.
12
,
209
(
2016
).
44.
A.
Kumar
and
T. D.
Crawford
,
J. Phys. Chem. A
121
,
708
(
2017
).
45.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
46.
S.
Grimme
,
J. G.
Brandenburg
,
C.
Bannwarth
, and
A.
Hansen
,
J. Chem. Phys.
143
,
054107
(
2015
).
47.
R.
Sure
and
S.
Grimme
,
J. Comput. Chem.
34
,
1672
(
2013
).
48.
F.
Neese
,
F.
Wennmohs
,
A.
Hansen
, and
U.
Becker
,
Chem. Phys.
356
,
98
(
2009
).

Supplementary Material