of 15
J. Chem. Phys.
151
, 064112 (2019);
https://doi.org/10.1063/1.5109882
151
, 064112
© 2019 Author(s).
Analytical gradients for projection-based
wavefunction-in-DFT embedding
Cite as: J. Chem. Phys.
151
, 064112 (2019);
https://doi.org/10.1063/1.5109882
Submitted: 13 May 2019 . Accepted: 03 July 2019 . Published Online: 09 August 2019
Sebastian J. R. Lee
, Feizhi Ding
, Frederick R. Manby
, and
Thomas F. Miller
COLLECTIONS
This paper was selected as an Editor’s Pick
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
Analytical gradients for projection-based
wavefunction-in-DFT embedding
Cite as: J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
Submitted: 13 May 2019
Accepted: 3 July 2019
Published Online: 9 August 2019
Sebastian J. R. Lee,
1
Feizhi Ding,
1
Frederick R. Manby,
2
and Thomas F. Miller III
1,a)
AFFILIATIONS
1
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, USA
2
Centre for Computational Chemistry, School of Chemistry, University of Bristol, Bristol BS8 1TS, United Kingdom
a)
tfm@caltech.edu
ABSTRACT
Projection-based embedding provides a simple, robust, and accurate approach for describing a small part of a chemical system at the level of
a correlated wavefunction (WF) method, while the remainder of the system is described at the level of density functional theory (DFT). Here,
we present the derivation, implementation, and numerical demonstration of analytical nuclear gradients for projection-based wavefunction-
in-density functional theory (WF-in-DFT) embedding. The gradients are formulated in the Lagrangian framework to enforce orthogonality,
localization, and Brillouin constraints on the molecular orbitals. An important aspect of the gradient theory is that WF contributions to the
total WF-in-DFT gradient can be simply evaluated using existing WF gradient implementations without modification. Another simplifying
aspect is that Kohn-Sham (KS) DFT contributions to the projection-based embedding gradient do not require knowledge of the WF calcula-
tion beyond the relaxed WF density. Projection-based WF-in-DFT embedding gradients are thus easily generalized to any combination of WF
and KS-DFT methods. We provide a numerical demonstration of the method for several applications, including a calculation of a minimum
energy pathway for a hydride transfer in a cobalt-based molecular catalyst using the nudged-elastic-band method at the coupled-cluster single
double-in-DFT level of theory, which reveals large differences from the transition state geometry predicted using DFT.
Published under license by AIP Publishing.
https://doi.org/10.1063/1.5109882
.,
s
I. INTRODUCTION
The theoretical description of many chemical processes
demands accurate,
ab initio
electronic structure theories. However,
the study of complex reactive processes, including those arising in
inorganic and enzyme catalysis, gives rise to the need for a compro-
mise between accuracy and the ability to complete the computation
in a reasonable amount of time. For systems in which the com-
plicated chemical rearrangements (e.g., bond breaking and form-
ing) occur in a local spatial region, an effective strategy for bal-
ancing accuracy and computational cost is to employ one of the
various multiscale embedding strategies.
1–25
Generally, embedding
methodologies hinge on the condition that a system can be effi-
ciently partitioned into a local subsystem that demands a high-level
treatment and an environment that can be treated with a lower (and
computationally less expensive) level of theory.
The current paper focuses on projection-based embedding,
9,26
a density functional theory (DFT)-based embedding theory in
which subsystem partitioning is performed in terms of localized
Kohn-Sham (KS) molecular orbitals (LMOs). The method describes
subsystem interactions at the level of KS and allows for the par-
titioning of the subsystems across covalent and even conjugated
bonds, and it enables the use of relatively small subsystem sizes for
an embedded wavefunction (WF) description. A recent review of
projection-based wavefunction-in-density functional theory (WF-
in-DFT) embedding is available in Ref. 26.
Projection-based embedding has proven to be a useful tool
in a wide range of chemical contexts including transition-metal
complexes,
19,20,24,27
protein active sites,
21,23
excited states,
28–30
and
condensed phase systems,
16
among others.
31–37
The development
of analytical nuclear gradients for projection-based embedding will
expand its applicability to include geometry optimization, transi-
tion state searches, and potentially
ab initio
molecular dynamics.
Analytical nuclear gradients already exist for a number of other
embedding methodologies, including the incremental molecular
fragmentation method,
38
fragment molecular orbital method,
39–41
quantum
mechanics/molecular
mechanics
(QM/MM),
42–44
ONIOM,
45–48
embedded mean-field theory (EMFT),
17,49
frozen
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-1
Published under license by AIP Publishing
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
density embedding,
50–53
and subsystem DFT.
54–56
However, the
projection-based approach provides a number of advantages for
WF-in-DFT embedding calculations and leads to a distinct analyt-
ical gradient theory, which we derive and numerically demonstrate
in several applications.
In Sec. II A, we outline projection-based WF-in-DFT embed-
ding, and in Sec. II B, we provide the derivation of its analytical
nuclear gradients. Section IV numerically validates the analytical
nuclear gradient theory and its implementation in Molpro
57
via
comparison with finite difference calculations, as well as present-
ing results for optimizing geometries in benchmark systems and the
calculation of a minimum energy profile for an organometallic reac-
tion using the nudged-elastic-band (NEB) method. We addition-
ally provide the analytical nuclear gradient theory for WF-in-DFT
embedding with atomic orbital (AO) truncation
18
in Appendixes C
and D.
II. PROJECTION-BASED EMBEDDING ANALYTICAL
NUCLEAR GRADIENTS
A. Projection-based embedding energy theory
Projection-based WF-in-DFT embedding relies on the parti-
tioning of the LMOs of a system into two subsystems. Subsys-
tem A contains the LMOs that are treated using the WF method,
and subsystem B contains the remaining LMOs that are treated
using KS. This WF-in-DFT procedure is accomplished by first per-
forming a KS calculation on the full system to obtain a set of
KS MOs. The occupied KS MOs are then localized and parti-
tioned into subsystems A and B. Finally, subsystem A is treated
using the WF method in the presence of the embedding poten-
tial created by the frozen LMOs of subsystem B. Note that the
cost of the KS calculation on the full system is typically neg-
ligible in comparison with the subsystem WF calculation. This
results in our working equation for projection-based WF-in-DFT
embedding,
9
E
WF-in-DFT
[
̃
Ψ
A
;
γ
A
,
γ
B
]
=
E
WF
[
̃
Ψ
A
]
+ tr
[(
̃
d
A
γ
A
)
v
emb
[
γ
A
,
γ
B
]
]
+
E
DFT
[
γ
A
+
γ
B
]
E
DFT
[
γ
A
]
+
μ
tr
[
̃
d
A
P
B
]
,
(1)
where
̃
Ψ
A
and
E
WF
[
̃
Ψ
A
]
are the WF and energy of subsystem A,
̃
d
A
is the subsystem A one-particle reduced density matrix that corre-
sponds to
̃
Ψ
A
,
E
DFT
is the KS energy, and
γ
A
and
γ
B
are, respectively,
the KS subsystem A and B one-particle densities that equal the full
system KS density,
γ
, when summed together. Throughout, we shall
use a tilde to indicate quantities that have been calculated using the
WF method. The embedding potential,
v
emb
, is defined as
v
emb
[
γ
A
,
γ
B
]
=
g
[
γ
A
+
γ
B
]
g
[
γ
A
]
,
(2)
where
g
includes all KS two-electron terms,
(
g
[
γ
])
κν
=
λσ
γ
λσ
((
κν
λσ
)
1
2
x
f
(
κλ
νσ
))
+
(
v
xc
[
γ
])
κν
,
(3)
and where
κ
,
ν
,
λ
, and
σ
label atomic orbital basis functions, (
κν
|
λσ
)
are the two-electron repulsion integrals,
x
f
is the fraction of exact
exchange, and
v
xc
is the exchange-correlation (XC) potential matrix.
The level-shift operator,
μ
P
B
, is given by
μ
P
B
=
μ
S
γ
B
S
,
(4)
where
S
is the overlap matrix. In the limit of
μ
, the LMOs that
make up subsystems A and B are enforced to be exactly orthogonal,
eliminating the nonadditive kinetic energy present in other embed-
ding frameworks.
58,59
In practice, finite values of
μ
in the range of 10
4
hartree to 10
7
hartree are found to provide accurate results regard-
less of chemical system.
9
If greater accuracy is needed, a perturbative
correction outlined in Ref. 9 can be added to the WF-in-DFT energy
expression to account for the finiteness of
μ
, but in practice, this cor-
rection is found to contribute negligibly to the total energy and is
thus neglected here.
Projection-based embedding can also be used for DFT-in-DFT
embedding via a simplified version of Eq. (1). The working equation
for projection-based DFT-in-DFT embedding is
9
E
DFT-in-DFT
[
̃
γ
A
;
γ
A
,
γ
B
]
=
E
DFT
[
̃
γ
A
]
+ tr
[(
̃
γ
A
γ
A
)
v
emb
[
γ
A
,
γ
B
]]
+
E
DFT
[
γ
A
+
γ
B
]
E
DFT
[
γ
A
]
+
μ
tr
[
̃
γ
A
P
B
]
.
(5)
The only differences between WF-in-DFT and DFT-in-DFT embed-
ding is that the first term on the right-hand side (RHS) of Eq. (1) is
replaced with the KS energy on subsystem A,
E
DFT
[
̃
γ
A
]
, and in the
second and last terms,
̃
d
A
is reduced to the subsystem A KS density
matrix,
̃
γ
A
.
B. Projection-based WF-in-DFT embedding
gradient theory
Since projection-based embedding is a nonvariational theory,
its analytical gradient is conveniently derived using a Lagrangian
approach. We first construct a Lagrangian based on the projection-
based WF-in-DFT energy. We then minimize the Lagrangian with
respect to the variational parameters in the embedding energy,
which include the subsystem A WF and the LMO coefficients. Then,
we show how to solve for each of the Lagrange multipliers and
provide the working equation for the gradient of the total energy.
For consistency in notation, the MO coefficient matrix
C
refers
to the entire set of KS MOs (occupied and virtual). The submatrix
of
C
that refers to the (occupied) LMOs is denoted as
L
with column
indices
i
,
j
,
k
,
l
. The submatrix of
C
that refers to the canonical virtual
space is denoted as
C
v
with column indices
a
,
b
,
c
,
d
. The indices
m
,
n
,
p
,
q
are used to index generic molecular orbitals.
C. Total energy Lagrangian
We now derive the total energy Lagrangian for projection-
based WF-in-DFT embedding. Where appropriate, we will provide
WF method specific examples (e.g., MP2) of general terms outlined
in the equations. The WF-in-DFT Lagrangian is
L
[
C
,
̃
Ψ
A
,
Λ
,
x
,
z
loc
,
z
]
=
E
WF-in-DFT
[
̃
Ψ
A
;
γ
A
,
γ
B
]
+
s
Λ
WF,A
s
c
s
+
pq
x
pq
(
C
SC
1
)
pq
+
i
>
j
z
loc
ij
r
ij
+
ai
z
ai
(
F
[
γ
A
+
γ
B
])
ai
.
(6)
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-2
Published under license by AIP Publishing
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
The first term on the right-hand side (RHS) of Eq. (6) is
the projection-based WF-in-DFT embedding energy described by
Eq. (1). The second term on the RHS of Eq. (6) contains any con-
straints,
c
s
, and the corresponding Lagrange multipliers,
Λ
WF,A
s
, that
arise from ensuring that the Lagrangian is variational with respect
to parameters in the WF method. The third term on the RHS con-
strains the KS MOs,
C
, to be orthonormal, which accounts for the
basis set being atom centered; this term is commonly referred to as
the Pulay force
60
and arises from the atomic orbital basis set being
atom centered. The localization conditions,
r
ij
= 0, take into account
how the KS MOs are localized before being selected for subsystems
A and B. This is important because the LMOs will have a differ-
ent dependence on nuclear perturbation than canonical MOs. In
this work, we use Pipek-Mezey localization
61
to obtain LMOs. Gen-
eralization to other localization methods (e.g., Boys
62
and intrinsic
bond orbitals
63
) is straightforward. The localization conditions for
Pipek-Mezey are
r
ij
=
C
(
S
C
ii
S
C
jj
)
S
C
ij
=
0
for all
i
>
j
,
(7)
where
C
corresponds to an atom in the molecule. The matrices
S
C
are defined as
S
C
kl
=
α
C
β
(
L
α
k
S
αβ
L
β
l
+
L
α
l
S
αβ
L
β
k
)
,
(8)
where the summation over
α
is restricted to basis functions at atom
C
. The Brillouin conditions,
(
F
[
γ
A
+
γ
B
]
)
ai
=
0, reflect how the KS
MOs are optimized before being used to construct subsystems A and
B. The Brillouin conditions are only needed because subsystem B is
frozen at the KS level of theory. However, due to the nonadditivity
of the XC potential, the Lagrange multipliers,
z
, span the full virtual-
occupied space.
The type and number of constraints applied to the WF method
depend on the chosen method. For example, if the WF method is
MP2, then the constraints are
s
Λ
MP2,A
s
c
s
=
pq
̃
x
pq
(
̃
C
A
S
̃
C
A
1
)
pq
+
ai
̃
z
ai
(
F
A
)
ai
i
A
,
(9)
where the first term on the RHS of Eq. (9) constrains the Hartree-
Fock (HF) MOs,
̃
C
A
, to be orthonormal, the condition
i
A restricts
the sum to occupied MOs in subsystem A, and the second term on
the RHS is the Brillouin conditions using the embedded Fock matrix,
F
A
. The embedded Fock matrix is defined as
9
F
A
=
h
+
g
[
̃
γ
A
]
+
v
emb
[
γ
A
,
γ
B
]
+
μ
P
B
,
(10)
where
h
is the standard one-electron Hamiltonian,
g
includes all of
the usual HF two-electron terms, and
̃
γ
A
is the subsystem A HF one-
particle density. These constraints also arise in the derivation of the
MP2 analytical nuclear gradient.
64
For the projection-based WF-in-DFT energy to be equal the
Lagrangian, the Lagrangian must be minimized with respect to all of
its parameters, including
̃
Ψ
A
,
C
, and all of the Lagrange multipliers.
D. Minimizing the Lagrangian with respect
to the variational parameters of the WF method
Upon minimizing the WF-in-DFT Lagrangian with respect to
̃
Ψ
A
, only terms associated with the first two terms on the RHS of
Eq. (6) survive, all of which are familiar from the WF Lagrangian for
the corresponding WF gradient theories,
L
̃
Ψ
A
=
E
WF
[
̃
Ψ
A
]
̃
Ψ
A
+ tr
̃
d
A
̃
Ψ
A
v
emb
[
γ
A
,
γ
B
]
+
μ
tr
̃
d
A
̃
Ψ
A
P
B
+
̃
Ψ
A
s
Λ
WF,A
s
c
s
=
0.
(11)
Since the embedding potential is independent of
̃
Ψ
A
, the Z-vector
coupled perturbed Hartree-Fock (Z-CPHF) equations of any post-
HF method are only impacted through the eigenvalues of the sub-
system A HF WF. Therefore, the solutions for the WF Lagrange
multipliers [e.g.,
̃
x
and
̃
z
for MP2 in Eq. (9)] are obtained using the
standard implementation of the WF gradient no matter what KS
method is selected to describe subsystem B. However, if an alter-
native embedding potential is used that depends on the subsystem
A WF, such as the Huzinaga constraint (i.e., Ref. 65), then the for-
mulation of the WF gradient is changed; the Z-CPHF equations for a
general WF method would need to be modified to include the contri-
butions from the derivative of the embedding potential with respect
to the subsystem A WF,
̃
Ψ
A
.
E. Minimizing the Lagrangian with respect
to the MO coefficients
The remaining Lagrange multipliers,
z
loc
,
x
, and
z
in Eq. (6), are
determined by minimizing the WF-in-DFT Lagrangian with respect
to the variational parameters of the KS method, namely, the MO
coefficients,
C
. Differentiation of the Lagrangian with respect to
these parameters yields
μ
C
μ
p
L
C
μ
q
=
E
pq
+
(
a
[
z
loc
])
pq
+
(
D
[
z
])
pq
+ 2
x
pq
=
0,
(12)
where
E
pq
=
μ
C
μ
p
E
WF-in-DFT
[
̃
Ψ
A
;
γ
A
,
γ
B
]
C
μ
q
+
C
μ
q
s
Λ
WF,A
s
c
s
,
(13)
(
a
[
z
loc
])
pq
=
μ
C
μ
p
k
>
l
z
loc
kl
r
kl
C
μ
q
=
k
>
l
B
pqkl
z
loc
kl
q
occ
,
(14)
(
D
[
z
])
pq
=
μ
C
μ
p
ak
z
ak
(
F
[
γ
A
+
γ
B
]
)
ak
C
μ
q
=
ak
D
pqak
z
ak
=
(
F
[
γ
A
+
γ
B
]
z
)
pq
q
occ
+
(
F
[
γ
A
+
γ
B
]
z
)
pq
q
vir
+ 2
(
V
[
̄
z
])
pq
q
occ
,
(15)
and
2
x
pq
=
μ
C
μ
p
mn
x
mn
S
mn
C
μ
q
.
(16)
The 4-dimensional tensors,
B
and
D
, are expanded in Appendixes A
and B, respectively,
̄
z
corresponds to
z
+
z
, and
V
[
̄
z
]
includes all
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-3
Published under license by AIP Publishing
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
two-electron terms of the generalized Fock matrix and is shown
explicitly in Appendix B. Since the embedded Fock matrix,
F
A
, con-
tains the embedding potential,
v
emb
, its derivative with respect to the
MO coefficients,
C
, is nonzero resulting in the WF relaxed density
being needed to construct
E
in Eq. (13), which is explicitly shown
in Appendix B. Therefore, the subsystem A WF gradient only affects
the embedding contributions to the gradient through the WF relaxed
density.
We now show that solving for the Lagrange multipliers leads
to familiar coupled perturbed equations. Combining the stationary
conditions described by Eq. (12) with the auxiliary conditions
x
=
x
yields the linear Z-vector equations,
64
(
1
P
pq
)(
E
+
D
[
z
]
+
a
[
z
loc
])
pq
=
0,
(17)
where
P
pq
permutes the indices
p
and
q
, which is used to solve for
z
and
z
lo c
. The matrix
x
is then obtained as
x
pq
=
1
4
(
1 +
P
pq
)(
E
+
D
[
z
]
+
a
[
z
loc
])
pq
.
(18)
The Lagrange multipliers
z
loc
pertain to the occupied-occupied MO
space; considering only the occupied-occupied part of Eq. (17)
yields
(
1
P
ij
)(
E
+
D
[
z
]
+
a
[
z
loc
])
ij
=
0.
(19)
Using the Brillouin conditions and the knowledge that
z
ab
=
z
ij
=
z
ia
= 0, Eq. (19) can be further simplified by showing that
(
1
P
ij
)(
D
[
z
])
ij
=
0.
(20)
The solutions,
z
lo c
, are thus independent of
z
such that Eq. (19)
reduces to
E
ij
E
ji
+
k
>
l
(
B
ijkl
B
jikl
)
z
loc
kl
=
0.
(21)
These are the Z-vector coupled perturbed localization (Z-CPL)
equations, which are used to solve for
z
lo c
. Subsequently,
a
[
z
loc
]
can
be computed according to Eq. (14).
The Lagrange multipliers
z
pertain to the virtual-occupied
MO space; considering only the virtual-occupied part of Eq. (17)
yields
(
1
P
ai
)(
E
+
D
[
z
]
+
a
[
z
loc
])
ai
=
0,
(22)
which further simplifies to
(
E
+
a
[
z
loc
]
+
F
[
γ
A
+
γ
B
]
z
zF
[
γ
A
+
γ
B
]
+ 2
V
[
̄
z
]
)
ai
=
0.
(23)
These are the Z-vector coupled perturbed Kohn-Sham (Z-CPKS)
equations. Having solved the Z-CPL and Z-CPKS equations, the
remaining Lagrangian multipliers associated with the orthogonality
constraints,
x
, can be obtained from Eq. (18).
F. Gradient of the total energy
Once the Lagrangian is minimized with respect to all varia-
tional parameters, the gradient of the total energy takes the form
d
E
WF-in-DFT
d
q
=
d
L
d
q
=
L
q
+
L
̃
Ψ
A
̃
Ψ
A
q
+
L
C
C
q
=
L
q
.
(24)
Since the Lagrangian is minimized with respect to the subsystem
A WF and the KS LMO coefficients, calculation of the WF and
KS LMO responses to nuclear perturbation,
̃
Ψ
A
/
q
and
C
/
q
,
respectively, are avoided. This yields the following expression for the
gradient:
E
(
q
)
WF-in-DFT
=
E
(
q
)
WF
[
̃
Ψ
A
]
+
λν
̃
d
A
λν
(
v
(
q
)
emb
)
λν
+
μ
λν
̃
d
A
λν
(
P
B,
(
q
)
)
λν
+
s
Λ
WF,A
s
c
(
q
)
s
E
(
q
)
DFT
[
γ
A
]
+
E
(
q
)
DFT
[
γ
A
+
γ
B
]
λν
γ
A
λν
(
v
(
q
)
emb
)
λν
+
i
>
j
z
loc
ij
r
(
q
)
ij
+
ai
z
ai
(
F
[
γ
A
+
γ
B
])
(
q
)
ai
+
mn
x
mn
S
(
q
)
mn
,
(25)
where the superscript (
q
) denotes the explicit derivative of the quan-
tity with respect to a nuclear coordinate. Equation (25) can be fur-
ther simplified by folding
s
Λ
WF,A
s
c
(
q
)
s
into the first three terms on
the RHS of Eq. (25), yielding
E
(
q
)
WF-in-DFT
=
E
q
WF
[
̃
Ψ
A
]
+
λν
(
̃
d
A
rel
)
λν
(
v
(
q
)
emb
)
λν
+
μ
λν
(
̃
d
A
rel
)
λν
(
P
B,
(
q
)
)
λν
E
(
q
)
DFT
[
γ
A
]
+
E
(
q
)
DFT
[
γ
A
+
γ
B
]
λν
γ
A
λν
(
v
(
q
)
emb
)
λν
+
i
>
j
z
loc
ij
r
(
q
)
ij
+
ai
z
ai
(
F
[
γ
A
+
γ
B
])
(
q
)
ai
+
mn
x
mn
S
(
q
)
mn
.
(26)
Here,
E
q
WF
[
̃
Ψ
A
]
denotes the total derivative of the subsystem A WF
energy with respect to nuclear coordinate, which can be directly cal-
culated using existing WF gradient implementations, and
̃
d
A
rel
is the
WF-relaxed density for subsystem A. For example, the MP2-relaxed
density is
̃
d
A
rel
=
̃
d
A
+
̃
C
A
̃
z
̃
C
A,
=
̃
γ
A
+
d
(
2
)
+
̃
C
A
̃
z
̃
C
A,
,
(27)
which contains the subsystem A Hartree-Fock density,
̃
γ
A
, the MP2
density matrix,
d
(2)
, and the solutions of the subsystem A Brillouin
conditions,
̃
C
A
̃
z
̃
C
A,
. Equation (26) can be expressed in terms of
the WF gradient on subsystem A and the derivative AO integrals,
yielding our final expression for the projection-based WF-in-DFT
analytical gradient,
E
(
q
)
WF-in-DFT
=
E
q
WF
[
̃
Ψ
A
]
+ tr
[
d
a
h
(
q
)
]
+ tr
[
XS
(
q
)
]
+
1
2
μνλσ
D
μνλσ
(
μν
λσ
)
(
q
)
+
E
(
q
)
xc
[
γ
A
+
γ
B
]
E
(
q
)
xc
[
γ
A
]
+ tr
[(
̃
d
A
rel
γ
A
)(
v
(
q
)
xc
[
γ
A
+
γ
B
]
v
(
q
)
xc
[
γ
A
]
)]
.
(28)
The effective one-particle density
d
a
and effective two-particle den-
sity
D
are defined
d
a
=
γ
B
+
CzC
,
(29)
and
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-4
Published under license by AIP Publishing
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
D
μνλσ
=
(
γ
A
+
γ
B
)
μν
(
d
b
)
λσ
γ
A
μν
(
d
c
)
λσ
1
2
x
f
(
(
γ
A
+
γ
B
)
μλ
(
d
b
)
νσ
γ
A
μλ
(
d
c
)
νσ
)
.
(30)
The effective one-particle densities
d
b
and
d
c
are defined
d
b
=
γ
A
+
γ
B
+ 2
CzC
+ 2
̃
d
A
rel
2
γ
A
,
(31)
and
d
c
=
γ
A
+ 2
̃
d
A
rel
.
(32)
The matrix
X
is defined
X
=
CxC
+
i
>
j
r
ij
S
μν
z
loc
ij
=
X
loc
1
2
L
(
E
+ 2
V
[
̄
z
]
)
L
1
2
(
C
v
(
zF
)
L
+
(
C
v
(
zF
)
L
)
)
+
μ
(
̃
d
A
rel
S
γ
B
+
γ
B
S
̃
d
A
rel
)
,
(33)
where
(
X
loc
)
μν
=
1
2
(
La
[
z
loc
]
L
)
μν
+
i
>
j
r
ij
S
μν
z
loc
ij
.
(34)
The second term on the RHS of Eq. (34) is expanded in Appendix A.
The analytical nuclear gradient expression for projection-based
DFT-in-DFT closely follows that for WF-in-DFT, with regard to
evaluation of both the Lagrange multipliers [Eq. (12)] and the
final gradient [Eq. (28)]. To obtain the corresponding DFT-in-DFT
expressions,
̃
d
A
rel
becomes the subsystem A KS density
̃
d
A
rel
=
̃
γ
A
,
(35)
which affects the evaluation of
E
in Eq. (12) [expanded in
Eq. (B1)] and the evaluation of the final gradient expression
[Eq. (28)]. Additionally, the first term on the RHS of the final gra-
dient expression [Eq. (28)], is replaced with the subsystem A KS
gradient,
E
q
DFT
[
̃
γ
A
]
.
III. COMPUTATIONAL DETAILS
The implementation of projection-based WF-in-DFT embed-
ding gradients is available in the 2019 general release of Mol-
pro.
57
In all embedding calculations reported here, unless other-
wise specified, the Pipek-Mezey localization method
61
is used with
the core and occupied MOs localized together. The subsystem A
region is chosen by including any LMOs with a net Mulliken pop-
ulation larger than 0.4 on the atoms associated with subsystem
A although more sophisticated partitioning algorithms have been
introduced.
27
A level-shift parameter of
μ
= 10
6
hartree is used for
all embedding calculations. The perturbative correction to using a
finite value of
μ
in Eq. (5) is less than 20 microhartrees for the
applications presented here and thus not included (accomplished
by specifying the option
HF
_
COR
= 0). Throughout this work, all
embedding calculations are described using the nomenclature “(WF
method)-in-DFT/basis,” where the WF method describes subsystem
A and the KS method describes subsystem B. For some embed-
ding calculations, a mixed-basis set is used and is denoted by “(WF
method)-in-DFT/large-basis:small-basis,” where the large basis is
used to describe subsystem A and the small basis is used to describe
subsystem B.
All SCF calculations employ a tighter threshold than default for
MO convergence by specifying the option
ORBITAL
= 1
×
10
7
a.u. in
Molpro. All KS calculations used in projection-based embedding are
done without density fitting, employing the local-density approxi-
mation (LDA),
66,67
Perdew-Burke-Ernzerhof (PBE),
68
PBE0,
69
and
LDAX functionals with the def2-TZVPP, def2-SVP, def2-ASVP,
70,71
cc-pVDZ,
72
and 6-31G
73
basis sets. Note that the def2-ASVP basis
set used in Molpro is constructed by adding one set of even tem-
pered diffuse functions to the def2-SVP basis set. The LDAX func-
tional is constructed by including 50% exact exchange and reducing
the weight of the DIRAC functional to 50% in the LDA functional.
For the calculations in Secs. IV A and IV B 2, the XC functional
is evaluated on a fixed-pruned grid with index 7 (Ref. 74). For the
optimized geometries shown in Sec. IV B, and the malondialdehyde
calculations in Sec. IV C, the XC functional is evaluated on an adap-
tively generated quadrature grid that reproduces the energy of the
Slater-Dirac functional to a specified threshold accuracy of 10
10
E
h
.
All WF calculations are performed with the frozen-core approxima-
tion, without density fitting, employing the MP2,
75
coupled-cluster
singles and doubles (CCSD),
76,77
and coupled-cluster singles, dou-
bles, and perturbative triples [CCSD(T)]
78
correlation treatments
with the def2-TZVPP, def2-SVP, cc-pVDZ, and 6-31G basis sets.
Even though the density fitting approximation is not used for the
WF methods in this study, density fitted gradients are available for
the aforementioned WF methods.
79,80
The default values for integral
screening were used in Molpro. For all Z-CPKS calculations, an iter-
ative subspace solver employing the Davidson algorithm
81,82
is used
with a convergence threshold of 1
×
10
6
a.u. For all Z-CPHF calcu-
lations needed for the subsystem A WF gradient, an iterative solver
with a convergence threshold of 1
×
10
7
a.u. is used. Grid weight
derivatives are included for all gradient calculations involving the
XC functional and potential.
For all geometry optimizations, the number of LMOs in sub-
system A is kept unchanged throughout the optimization. A natural
way of enforcing this in future work is to employ even-handed par-
titioning
27
although this was not needed in the examples studied
here; the default procedure based on net Mulliken population suf-
ficed to keep subsystem A unchanged. All geometries are optimized
using the translation-rotation-internal coordinate system devised by
Wang and Song,
83
which is available in the GeomeTRIC package.
84
Convergence parameters for the geometry optimizations follow the
default parameters used by Molpro, namely, that the maximum gra-
dient value becomes less than 3
×
10
4
hartree/bohrs and the energy
change between adjacent steps becomes less than 1
×
10
6
hartree
or the maximum component of the step displacement becomes less
than 3
×
10
4
bohrs. The maximum gradient value is evaluated in the
Cartesian basis. All geometries are provided in the supplementary
material.
Nudged elastic band (NEB)
85
calculations are run using the
implementation of the method in the atomic simulation environ-
ment (ASE) package.
86
All NEB calculations use Molpro forces
which are provided through a Molpro calculator interface within the
ASE package.
The intramolecular proton transfer of malondialdehyde is
modeled with an NEB consisting of 15 images connected by
springs with spring constants of 0.1 eV/Å
2
. The CCSD/def2-aSVP,
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-5
Published under license by AIP Publishing
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
CCSD-in-LDA/def2-aSVP, and LDA/def2-aSVP optimized NEBs
used the image dependent pair potential (IDPP) method
87
as
the initial guess for the band with reactant and product geome-
tries previously optimized at the corresponding level of theory.
All NEB calculations for malondialdehyde are converged with the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) update of the Hessian
and by enforcing that the maximum gradient value is less than
0.01 eV/Å
2
.
The intramolecular proton transfer of the organometallic cobalt
complex is modeled with an NEB consisted of images connected
by springs with spring constants of 9 eV/Å
2
. The PBE0/cc-pVDZ
climbing image NEB,
88
consisting of 26 images, used the IDPP
method as the initial guess for the band with the reactant geome-
try previously optimized. The CCSD-in-PBE0/cc-pVDZ NEB con-
sisting of 23 images used its optimized reactant, and the climbing
image NEB converged at the PBE0/cc-pVDZ level of theory as its
initial guess. Since the product is spatially far away from the reac-
tant, an intermediate geometry between the transition state and
the product is used as the endpoint of the NEB. This intermedi-
ate geometry is determined by initially converging a NEB with an
extra image such that the maximum force dropped below 0.3 eV/Å
2
.
Then, the second to last image is used as the new endpoint and a
new NEB is converged. The PBE0/cc-pVDZ climbing image NEB is
optimized using the FIRE
89
algorithm using a convergence criteria
of 0.05 eV/Å
2
for the maximum gradient value. The projection-
based WF-in-DFT embedding NEB is optimized at the CCSD-in-
PBE0/cc-pVDZ level of theory using the FIRE algorithm with a
convergence criterion of 0.25 eV/Å
2
for the maximum gradient
value.
The CCSD-in-PBE0/cc-pVDZ calculations used for the NEB
optimization are performed by specifying 51 occupied MOs to
be in subsystem A using the N_ORBITALS option and by using
AO truncation with a threshold of 1
×
10
3
a.u. An even-handed
selection of AOs was used along the NEB by creating a union
of the AOs that were selected for each image by the truncation
procedure to ensure that the NEB traversed a smooth potential
energy surface. The final, reported energies of the WF-in-DFT
NEB are performed using the PNO-LCCSD
90
/cc-pVDZ and PNO-
LCCSD-in-PBE0/cc-pVDZ levels of theory. Both the PNO-LCCSD
and PNO-LCCSD-in-PBE0 calculations are performed with den-
sity fitting using the cc-pVTZ/JKFIT
91
(the def2-TZVPP/JKFIT
92
basis set was used for cobalt since the cc-pVTZ/JKFIT basis set
was not available) and the cc-pVTZ/MP2FIT
93
density fitting basis
sets. Tighter domain approximations were employed for all PNO-
LCCSD calculations by specifying the
DOMOPT
=
TIGHT
option.
Additionally, the Boughton-Pulay completeness criterion was used
for the selection of the primary projected atomic orbital domain
by specifying the option
THRBP = 1
, and the Pipek-Mezey local-
ization method was used. For the PNO-LCCSD-in-PBE0 calcu-
lations, AO truncation is not used, the core and valence DFT
molecular orbitals are localized separately using the Pipek-Mezey
localization method, and the subsystem A orbitals are selected
using the default procedure based on the Mulliken population
threshold.
All calculations using AO truncation
18
ensure that at least
one AO is kept per atom (specified by option AO_PER_ATOM)
to make evaluation of the integral derivative contributions from
the one electron Hamiltonian simpler within Molpro. This adds a
negligible amount of AO functions than would have been selected
using only the density threshold parameter
18
for the systems stud-
ied in this paper. In all embedding geometry optimizations that
employ AO truncation, the number of truncated AOs is fixed using
the STOREAO option to ensure smoothness of the potential energy
function. Upon convergence, the truncated AO list is reevaluated
using the same density threshold parameter; if the number of kept
AOs remains a subset of the original list of truncated AOs, then the
optimization is converged.
IV. RESULTS AND DISCUSSION
A. Comparison of analytical and numerical gradients
The implementation of the projection-based WF-in-DFT ana-
lytical gradient is tested by comparison with the gradient eval-
uated by numerical finite difference for a distorted geometry of
ethanol. The finite difference gradients are evaluated using a four-
point central difference formula with a base step size of 0.01 bohr.
The mean absolute error (MAE) between the analytical and finite
difference gradients is reported for a range of embedding calcu-
lations in Table I. These results show that the analytical nuclear
gradient for projection-based WF-in-DFT embedding is essentially
numerically exact with respect to the gradients calculated by finite
difference. A comparison of the results obtained using HF over
the full system vs using LDA over the full system illustrates that
some of the finite difference error comes from the DFT exchange-
correlation grid. A comparison of the HF-in-HF results with full HF
and of the LDA-in-LDA results with full LDA illustrates the mod-
est effect of using a large-but-finite value for the level-shift operator
in projection based embedding. These results confirm the correct
implementation of projection-based WF-in-DFT analytical nuclear
gradients.
TABLE I
. Mean absolute error between the analytically and numerically determined
embedding nuclear gradient for a distorted geometry of ethanol. The basis set 6-
31G is used for all calculations. The distorted geometry of ethanol is provided in the
supplementary material.
Method
MAE (hartree/bohrs)
HF
5.00
×
10
9
HF-in-HF
4.61
×
10
8
LDA
1.48
×
10
8
LDA-in-LDA
7.23
×
10
8
HF-in-LDA
5.24
×
10
8
MP2-in-LDA
5.37
×
10
8
CCSD-in-LDA
5.36
×
10
8
CCSD(T)-in-LDA
5.26
×
10
8
CCSD-in-LDA (AO)
a
3.48
×
10
8
CCSD(T)-in-LDA (AO)
a
3.40
×
10
8
CCSD(T)-in-PBE0
5.26
×
10
8
CCSD(T)-in-PBE0 (AO)
a
1.12
×
10
7
a
Calculations were performed with AO truncation with a density threshold of
1
×
10
1
a.u.
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-6
Published under license by AIP Publishing
The Journal
of Chemical Physics
ARTICLE
scitation.org/journal/jcp
FIG. 1
. Optimized geometry for ethanol using projection-based CCSD-in-LDA/6-
31G. The solid atoms (O and H) are in subsystem A, and the transparent atoms
are in subsystem B.
TABLE II
. Selected bond lengths and angles for ethanol (pictured in Fig. 1) optimized
at different levels of theory. Bond lengths are reported in units of angstroms, and
angles are reported in units of degrees.
Method
r(O–H)
C
1
OH
r(C
1
–C
2
)
r(C
1
–O)
LDA/6-31G
0.988
110.4
1.503
1.439
CCSD-in-LDA/6-31G
0.979
110.7
1.506
1.435
CCSD/6-31G
0.979
110.6
1.532
1.475
B. Optimized geometries
1. Ethanol
As a proof of concept, CCSD-in-LDA/6-31G analytical nuclear
gradients are employed to determine the ground state geometry of
ethanol, which is shown in Fig. 1. For this simple case, the O–H
moiety is treated by CCSD and the remainder of the molecule is
treated by LDA. Table II shows that the O–H bond length within
subsystem A reproduces the CCSD predicted bond length of 0.979 Å
and the remaining bonds within subsystem B reproduce the LDA
predicted bond lengths. This indicates that the potential energy sur-
face produced by projection-based embedding varies smoothly from
CCSD-like interactions for subsystem A and LDA-like interactions
TABLE III
. Selected bond lengths for the organometallic complex pictured in Fig. 2
optimized at different levels of theory and their absolute difference (|
Δ
|). Bond lengths
are reported in units of angstroms.
LDAX
CCSD-in-LDAX
|
Δ
|
Subsystem A
r(Co–N
1
)
1.836
1.846
0.010
r(Co–N
2
)
1.893
1.883
0.010
r(Co–N
3
)
1.932
1.951
0.019
r(Co–N
4
)
1.900
1.926
0.026
r(Co–N
5
)
1.978
2.026
0.048
Boundary
r(N
5
–O
1
)
1.317
1.355
0.038
r(N
5
–O
2
)
1.262
1.301
0.039
r(C
1
–N
5
)
1.131
1.150
0.019
r(N
2
–C
2
)
1.430
1.458
0.028
r(N
3
–C
3
)
1.428
1.458
0.030
Subsystem B
r(O
1
–H)
1.025
1.030
0.005
for subsystem B. Interestingly, the C–O bond located at the bound-
ary between subsystems A and B closely reproduces the LDA bond
length and is not an interpolation between the CCSD and LDA bond
lengths.
2. Cobalt-based organometallic complex
As a demonstration of embedding gradients with AO trun-
cation, the geometry of the cobalt-based organometallic complex,
shown in Fig. 2, is optimized. Figure 2(a) shows the CCSD-in-
LDAX/def2-TZVPP:def2-SVP optimized structure of the cobalt
complex where the solid atoms are included in subsystem A and
the transparent atoms are included in subsystem B. In Fig. 2(b),
the optimized structures evaluated at the CCSD-in-LDAX/def2-
TZVPP:def2-SVP (solid) and the LDAX/def2-TZVPP:def2-SVP
(transparent) levels of theory are overlaid. While only modest dif-
ferences are seen in the overall structure, Table III shows that
the optimized bond lengths do change between the two levels of
theory, both for the region within subsystem A and at the sub-
system boundary. This indicates that the WF method is capa-
ble of relaxing the atoms in subsystem A even when they are
strongly coordinated with subsystem B. It is also seen that the bond
lengths across the boundary of subsystems A and B also differ from
FIG. 2
. (a) The optimized geometry for
the cobalt-based organometallic com-
plex performed with projection-based
CCSD-in-LDAX/def2-TZVPP:def2-SVP
with AO truncation. The solid atoms (Co,
N
1
, N
2
, N
3
, N
4
, N
5
, and C
1
) are included
in subsystem A, and the transparent
atoms are included in subsystem B. (b)
The LDAX/def2-TZVPP:def2-SVP opti-
mized geometry (transparent) and the
projection-based CCSD-in-LDAX/def2-
TZVPP:def2-SVP with AO truncation
optimized geometry (solid).
J. Chem. Phys.
151
, 064112 (2019); doi: 10.1063/1.5109882
151
, 064112-7
Published under license by AIP Publishing