of 13
Analytical Gradients for Projection-Based Wavefunction-in-DFT Embedding
Sebastian J. R. Lee,
1
Feizhi Ding,
1
Frederick R. Manby,
2
and Thomas F. Miller III
1
1)
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125,
United States
2)
Centre for Computational Chemistry, School of Chemistry, University of Bristol, Bristol BS8 1TS,
United Kingdom
(Dated: 15 March 2019)
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 method while the remainder of the system
is described at the level of density functional theory. 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 calculation 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 numerical demonstration of
the method for benchmark systems, including calculation of a minimum energy pathway for intramolecular
hydrogen-atom transfer reaction in malondialdehyde using the nudged-elastic-band method at the CCSD-in-
DFT level of theory.
I.
INTRODUCTION
The theoretical description of many chemical processes
demands accurate,
ab initio
electronic structure theories.
However, the study of complex reactive processes, includ-
ing those arising in inorganic and enzyme catalysis, gives
rise to the need for a compromise between accuracy and
the ability to complete the computation in a reasonable
amount of time. For systems in which the complicated
chemical rearrangements (e.g. bond breaking and form-
ing) occurs in a local spatial region, an effective strategy
for balancing accuracy and computational cost is to em-
ploy one of various multiscale embedding strategies.
1–24
Generally, embedding methodologies hinge on the con-
dition that a system can be efficiently 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
a DFT-based embedding theory in which
subsystem partitioning is performed in terms of local-
ized Kohn-Sham (KS) molecular orbitals (LMOs). The
method describes subsystem interactions at the level of
KS and allows for the partitioning of the subsystems
across covalent and even conjugated bonds, and it en-
ables the use of relatively small subsystem sizes for an
embedded WF description.
Projection-based embedding has proven to be a use-
ful tool in a wide range of chemical contexts in-
cluding transition-metal complexes,
19,20,24,25
protein ac-
tive sites,
21,23
excited states
26–28
and condensed phase
systems,
16
among others.
29–35
The development of ana-
lytical nuclear gradients for projection-based embedding
will expand its applicability to include geometry opti-
mization, transition state searches, and potentially
ab in-
tio
molecular dynamics. Analytical nuclear gradients al-
ready exist for a number of other embedding methodolo-
gies, including the incremental molecular fragmentation
method,
36
QM/MM,
37–39
ONIOM,
40,41
EMFT,
17,42
and
frozen density embedding.
43–45
However, the projection-
based approach provides a number of advantages for WF-
in-DFT embedding calculations and leads to a distinct
analytical gradient theory, which we derive and numeri-
cally demonstrate in benchmark systems.
In section II A we outline projection-based WF-in-DFT
embedding and in section II B we provide the deriva-
tion of its analytical nuclear gradients. Section IV nu-
merically validates the analytical nuclear gradient the-
ory and its implementation in Molpro
46
via compari-
son with finite difference calculations, as well as present-
ing results for optimizing geometries in benchmark sys-
tems and the calculation of a minimum energy profile for
an organometallic reaction using the nudged-elastic-band
(NEB) method. We additionally provide the analytical
nuclear gradient theory for WF-in-DFT embedding with
atomic orbital (AO) truncation
18
in Appendices 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
partitioning the LMOs of a system into two subsystems.
Subsystem A contains the LMOs that are treated using
the WF method and subsystem B contains the remain-
ing LMOs that are treated using KS. This WF-in-DFT
arXiv:1903.05830v1 [physics.chem-ph] 14 Mar 2019
2
procedure is accomplished by first performing a KS cal-
culation on the full system to obtain a set of KS MOs.
The occupied KS MOs are then localized and partitioned
into subsystems A and B. Finally, subsystem A is treated
using the WF method in the presence of the embedding
potential created by the frozen LMOs of subsystem B.
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 sub-
system A,
̃
d
A
is the subsystem A one-particle reduced
density matrix that corresponds to
̃
Ψ
A
,
E
DFT
is the KS
energy, and
γ
A
and
γ
B
are respectively the KS subsys-
tem A and B one-particle densities that equal the full
system KS density,
γ
, when summed together. Through-
out, 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 func-
tions, (
κν
|
λσ
) 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 opera-
tor,
μ
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 non-additive ki-
netic energy present in other embedding frameworks.
47,48
In practice, finite values of
μ
in the range of 10
4
a.u. to
10
7
a.u. 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 correction 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 em-
bedding 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 embedding is that the first term on the 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 re-
duced to the subsystem A KS density matrix,
̃
γ
A
.
B.
Projection-based WF-in-DFT Embedding Gradient
Theory
Since projection-based embedding is a non-variational
theory, its analytical gradient is conveniently derived us-
ing a Lagrangian approach. We first construct a La-
grangian based on the projection-based WF-in-DFT en-
ergy. We then minimize the Lagrangian with respect
to the variational parameters in the embedding energy,
which include the subsystem A WF and the LMO co-
efficients. 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 ma-
trix
C
refers to the entire set of KS MOs (occupied and
virtual). The submatrix of
C
that refers to the (occu-
pied) 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 or-
bitals.
C.
Total Energy Lagrangian
We now derive the total energy Lagrangian for
projection-based WF-in-DFT embedding. Where appro-
priate 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)
The first term on the right hand side (RHS) of Eq. 6 is
the projection-based WF-in-DFT embedding energy de-
scribed by Eq. 1. The second term on the RHS of Eq. 6
3
contains any constraints,
c
s
, and the corresponding La-
grange 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 ac-
counts for the basis set being atom centered; this term is
commonly referred to as the Pulay force
49
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 sub-
systems A and B. This is important because the LMOs
will have a different dependence on nuclear perturba-
tion than canonical MOs. In this work, we use Pipek-
Mezey localization
50
to obtain LMOs. Generalization
to other localization methods (e.g. Boys
51
and intrin-
sic bond orbitals
52
) 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 ba-
sis functions at atom
C
.
The Brillouin conditions,
(
F
[
γ
A
+
γ
B
])
ai
= 0, reflect how the KS MOs are opti-
mized before being used to construct subsystems A and
B. The Brillouin conditions are only needed because sub-
system B is frozen at the KS level of theory. However, due
to the non-additivity 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 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 are 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
in-
cludes all of the usual HF two-electron terms and
̃
γ
A
is the subsystem A HF one-particle density. These con-
straints also arise in the derivation of the MP2 analytical
nuclear gradient.
53
For the projection-based WF-in-DFT energy to 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 re-
spect to
̃
Ψ
A
, only terms associated with the first two
terms on the RHS of Eq. 6 survive, all of which are fa-
miliar 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 subsystem 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 alternative embedding potential is used
that depends on the subsystem A WF, such as the Huz-
inaga constraint (i.e. Ref. 54), then the formulation of
the WF gradient is changed; the Z-CPHF equations for a
general WF method would need to be modified to include
the contributions 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 La-
grangian with respect to the variational parameters of the
KS method, namely the MO coefficients,
C
. Differentia-
tion of the Lagrangian with respect to these parameters
yields
μ
C
μ,p
L
∂C
μ,q
=
E
pq
+
(
a
[
z
loc
])
pi
+ (
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)
4
(
a
[
z
loc
])
pi
=
μ
C
μ,p
(
k>l
z
loc
kl
∂r
kl
∂C
μ,q
)
=
k>l
B
pi,kl
z
loc
kl
,
(14)
(
D
[
z
])
pq
=
μ
C
μ,p
(
ak
z
ak
(
F
[
γ
A
+
γ
B
])
ak
∂C
μ,q
)
=
ak
D
pq,ak
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
Appendices A and B respectively,
̄
z
corresponds to
z
+
z
,
and
V
[
̄
z
] includes all two-electron terms of the general-
ized Fock matrix and is shown explicitly in Appendix B.
Since the embedded Fock matrix,
F
A
, contains the em-
bedding 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. Combin-
ing the stationary conditions described by Eq. 12 with
the auxiliary conditions
x
=
x
yields the linear Z-vector
equations
53
(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
loc
. 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
loc
, are thus independent of
z
, such that
Eq. 19 reduces to
E
ij
E
ji
+
k>l
(
B
ij,kl
B
ji,kl
)
z
loc
kl
= 0.
(21)
These are the Z-vector coupled perturbed localization (Z-
CPL) equations, which are used to solve for
z
loc
. Subse-
quently,
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
variational 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, calcula-
tion of the WF and KS LMO responses to nuclear per-
turbation,
̃
Ψ
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
]
+ tr
[
̃
d
A
v
(
q
)
emb
]
+
μ
tr
[
̃
d
A
P
B
,
(
q
)
]
+
s
Λ
WF,A
s
c
(
q
)
s
E
(
q
)
DFT
[
γ
A
] +
E
(
q
)
DFT
[
γ
A
+
γ
B
]
tr
[
γ
A
v
(
q
)
emb
]
+
ij
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 deriva-
tive of the quantity with respect to a nuclear coordinate.
Eq. 25 can be further simplified by folding
s
Λ
WF,A
s
c
(
q
)
s
5
into the first three terms on the RHS of Eq. 25, yielding
E
(
q
)
WF-in-DFT
=
E
q
WF
[
̃
Ψ
A
]
+ tr
[
̃
d
A
rel
v
(
q
)
emb
]
+
μ
tr
[
̃
d
A
rel
P
B
,
(
q
)
]
E
(
q
)
DFT
[
γ
A
] +
E
(
q
)
DFT
[
γ
A
+
γ
B
]
tr
[
γ
A
v
(
q
)
emb
]
+
ij
z
loc
ij
r
(
q
)
ij
+
ai
z
ai
(
F
[
γ
A
+
γ
B
])
(
q
)
ai
+
mn
x
mn
S
(
q
)
mn
,
(26)
where
E
q
WF
[
̃
Ψ
A
] denotes the total derivative of the sub-
system A WF energy with respect to nuclear coordinate,
which can be directly calculated using existing WF gradi-
ent 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
,
. Eq. 26
can be expressed in terms of the WF gradient on subsys-
tem A and the derivative AO integrals, yielding our final
expression for the projection-based WF-in-DFT analyti-
cal 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 density
D
are defined
d
a
=
γ
B
+
CzC
,
(29)
and
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 ex-
pressions,
̃
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 expres-
sion, Eq. 28. Additionally, the first term on the RHS of
the final gradient 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
embedding gradients is available in the 2019 general re-
lease of Molpro.
46
In all embedding calculations reported
here, the Pipek-Mezey localization method,
50
with the
core and occupied MOs localized together. The subsys-
tem A region is chosen by including any LMOs with a
net Mulliken population larger than 0.4 on the atoms
associated with subsystem A, although more sophisti-
cated partitioning algorithms have been introduced.
25
A level-shift parameter of
μ
= 10
6
a.u. is used for
all embedding calculations. The perturbative correc-
tion 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 op-
tion
HF
COR
= 0). Throughout this work, all embed-
ding calculations are described using the nomenclature
“(WF method)-in-DFT/basis,” where the WF method
describes subsystem A and the KS method describes sub-
system B. For some embedding 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.
6
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 LDA,
55,56
PBE,
57
PBE0,
58
and LDAX functionals with the def2-TZVPP, def2-SVP,
def2-aSVP,
59,60
cc-pVDZ,
61
and 6-31G
62
basis sets. The
LDAX functional is constructed by including 50% exact
exchange and reducing the weight of the DIRAC func-
tional to 50% in the LDA functional. For the calcula-
tions in section IV A, the XC functional is evaluated on
a fixed-pruned grid with index 7 (Ref. 63). For the opti-
mized geometries shown in section IV B, and the malon-
dialdehyde calculations in section IV C the XC functional
is evaluated on an adaptively 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 approx-
imation, without density fitting, employing the MP2,
64
CCSD,
65,66
and CCSD(T)
67
correlation treatments with
the def2-TZVPP, def2-SVP, cc-pVDZ and 6-31G basis
sets. The default values for integral screening were used
in Molpro. For all Z-CPKS calculations an iterative sub-
space solver employing the Davidson algorithm
68,69
is
used with a convergence threshold of 1
×
10
6
a.u. For
all Z-CPHF calculations 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 subsystem A is kept unchanged throughout the op-
timization. A natural way of enforcing this in future
work is to employ even-handed partitioning,
25
although
this was not needed in the examples studied here. All
geometries are optimized using the translation-rotation-
internal coordinate system devised by Wang and Song,
70
which is available in the GeomeTRIC package.
71
Conver-
gence parameters for the geometry optimizations follow
the default parameters used by Molpro, namely that the
maximum gradient value becomes less than 3
×
10
4
a.u.
and the energy change between adjacent steps becomes
less than 1
×
10
6
a.u. or the maximum component of
the step displacement becomes less than 3
×
10
4
a.u.
The maximum gradient value is evaluated in the Carte-
sian basis. All geometries are provided in the supporting
information.
Nudged elastic band (NEB)
72
calculations are run us-
ing the implementation of the method in the atomic sim-
ulation environment (ASE) package.
73
All NEB calcu-
lations 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 con-
nected by springs with spring constants of 0
.
1 eV/
̊
A
2
.
The CCSD/def2-aSVP, CCSD-in-LDA/def2-aSVP and
LDA/def2-aSVP optimized NEBs used the image depen-
dent pair potential (IDPP) method
74
as the initial guess
for the band with reactant and product geometries pre-
viously optimized at the corresponding level of theory.
All NEB calculations for malondialdehyde are converged
with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) up-
date of the Hessian and by enforcing that the maximum
gradient value is less than 0
.
01 eV/
̊
A
2
.
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 deriva-
tive contributions from the one electron Hamiltonian sim-
pler 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 studied 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 analytical gradient is tested by comparison with
the corresponding gradient evaluated by numerical finite
difference for a distorted geometry of ethanol. The finite
difference gradients are evaluated using a four-point cen-
tral 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 calculations 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 differ-
ence. Comparison of the results obtained using HF over
the full system versus using LDA over the full system
illustrate that some of the finite difference error comes
from the DFT exchange-correlation grid. Comparison of
the HF-in-HF results with full HF and of the LDA-in-
LDA results with full LDA illustrate the modest 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.
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
7
TABLE I. Mean absolute error between the analytically and
numerically determined embedding nuclear gradient for a dis-
torted geometry of ethanol. The basis set 6-31G is used for all
calculations. The distorted geometry of ethanol is provided
in the supporting information.
Method
MAE (a.u.)
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.
C
1
O
H
C
2
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.
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 subsys-
tem A reproduces the CCSD predicted bond length of
0
.
979
̊
A and the remaining bonds within subsystem B re-
produce the LDA predicted bond lengths. This indicates
that the potential energy surface produced by projection-
based embedding varies smoothly from CCSD-like inter-
actions for subsystem A and LDA-like interactions for
subsystem B. Interestingly, the C-O bond located at the
boundary between subsystems A and B closely repro-
duces the LDA bond length and is not an interpolation
between the CCSD and LDA bond lengths.
TABLE II. Selected bond lengths and angles for ethanol (pic-
tured 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
2.
Cobalt-based Organometallic Complex
As a demonstration of embedding gradients with
AO truncation, the geometry of the cobalt-based
organometallic complex, shown in Fig. 2, is optimized.
Panel (a) of Fig. 2 shows the CCSD-in-LDAX/def2-
TZVPP:def2-SVP optimized structure of the cobalt com-
plex where the solid atoms are included in subsystem
A and the transparent atoms are included in subsystem
B. In panel (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. The two geometries significantly dif-
fer in the structure of the ring network around the cobalt
center, with the difference being most pronounced at the
peripheral carbon centers. Additionally, Table III shows
that there is a large difference in bond lengths obtained
with the two levels of theory, both for the region within
subsystem A and at the boundary. This indicates that
the WF method is capable of relaxing the atoms in sub-
system 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 the
LDAX geometry since the bonds in question experience
the effects of both the WF and KS methods. Finally, if
a bond length associated with atoms in subsystem B is
considered, such as the O
1
-H bond, it is found to closely
match the LDAX predicted bond length.
(a)
(b)
O
1
N
1
O
2
N
2
N
3
N
4
N
5
C
1
C
2
C
3
FIG. 2.
Panel (a) shows the optimized geometry
for the cobalt-based organometallic complex 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. Panel
(b) shows the LDAX/def2-TZVPP:def2-SVP optimized ge-
ometry (transparent) and the projection-based CCSD-in-
LDAX/def2-TZVPP:def2-SVP with AO truncation optimized
geometry (solid). The root-mean-square deviation in atomic
position between LDAX/def2-TZVPP:def2-SVP and CCSD-
in-LDAX/def2-TZVPP:def2-SVP is 0
.
314
̊
A.
C.
Malondialdehyde: Minimum Energy Reaction Pathway
The minimum energy reaction pathway for the pro-
ton transfer in malondialdehyde is determined using the