Multi-site Reaction Dynamics Through Multi-fragment Density Matrix Embedding

The practical description of disordered chemical reactions, where reactions involve multiple species at multiple sites, is presently challenge using correlated electronic structure methods. Here we describe the gradient theory of multi-fragment density matrix embedding theory, which potentially provides a minimal computational framework to model such processes at the correlated electron level. We present the derivation and implementation of the gradient theory, its validation on model systems and chemical reactions using density matrix embedding, and its application to a molecular dynamics simulation of proton transport in a small water cluster, a simple example of multi-site reaction dynamics.


I. INTRODUCTION
Electron correlation plays a crucial role in the energetics of chemical reactions, but is expensive to model in large and complex systems. Quantum embedding methods 1-4 present a divide and conquer approach to the electron correlation problem. They proceed by defining small embedded regions that are described accurately, surrounded by environments that are described approximately. Typically, the environment is either represented by some type of potential acting on the embedded system, or is modeled by explicit, but reduced, quantum degrees of freedom (i.e. a bath) that help capture the fluctuations between the embedded region and its exterior. The density matrix embedding theory (DMET) [5][6][7][8][9][10][11] is a quantum embedding method in the second category. DMET provides an exact description of embedded fragments at a mean-field level and an approximate description at the correlated level, with the capability to describe systems broken into multiple fragments of arbitrary chemical nature and size, while ensuring a continuous potential energy surface (PES).
In the current work, we derive the analytical nuclear gradient theory of DMET. Although such a gradient theory, similar to the gradient theory of other embedding methods 12,13 , has multiple applications, our motivation is to study reactions in settings where the "reaction site" is not necessarily easily defined, e.g. in a disordered process where the atoms move freely and chemical reactions can occur at multiple locations, such as in nanoreactors 14 , combustion processes 15,16 , and heterogeneous catalysis 17,18 . Here, the ability of DMET to partition problems into arbitrary fragments while retaining a continuos PES, potentially positions it as a minimal theoretical framework to describe such phenomena. Below we first describe the working equations and implementation of the DMET gradient, then illustrate its applications to some simple chemical reactions, including a reactive molecular dynamics (MD) simulation of proton transfer in a small water cluster, a prototypical example of multi-site reaction dynamics.

II. THEORY
In the following, we adopt the Einstein summation convention, i.e., repeated indices imply summation unless otherwise noted.

A. The DMET energy expression
The theory of DMET is documented elsewhere 7,19 ; we only summarize some general ideas here to establish the energy expression and solution constraints in order to take the gradient.
The basic idea in DMET is to decompose a large system into a set of embedded fragments or impurities (we use the terms interchangeably) coupled to baths. Each bath represents the environment surrounding each fragment. As a practical approximation, DMET constructs the baths using a mean-field (i.e. Slater determinant) representation of the global ground state. The information in the global Slater determinant is contained in its mean-field one-particle reduced density matrix (1-RDM), . Here , , , index the computational atomic orbital (AO) basis; the mean-field canonical orbitals will be referred to as MOs. To facilitate the definition of impurities, we transform all quantities into an orthogonalized AO basis, which we term the localized orbitals (LOs), labeled , , , . As an example, the mean-field 1-RDM elements rotated into the LO basis are labelled .
An impurity is defined as a subset of the full LOs. We partition the full LO basis into multiple non-overlapping impurities. For impurity , the set of LOs not in defines its environment, which in DMET is later replaced by a set of bath and core orbitals. The bath of is obtained from the singular value decomposition (SVD) of the environment-impurity mean-field density matrix where ̃ transforms the environment LOs into the bath orbitals. Similarly, the core orbitals for are defined from the environment-environment block of , as the eigenvectors with double occupancy. Together, the set of impurity, bath, and core orbitals define the embedded impurity problem for ; the impurity and bath orbitals are usually referred to as the embedding orbitals (EOs) which we label with tilde ("∼") labels̃ ,̃ ,̃ ,̃ ,̃ ; the corresponding core orbitals are labelled̃ ,̃ . The embedding and core orbitals are, of course, different for each , but where possible we suppress the additional label to reduce notation. The transformation matrix from the LOs to the EOs of fragment can be summarized as The embedding Hamiltonian for impurity is defined by projecting the original Hamiltonian into the space of determinants spanned by {̃ }⊗|core⟩ (where |core⟩ is the doubly occupied core Slater determinant constructed from the core orbitals associated with ). Since there are different EOs for each impurity, there is a different for each impurity. In second quantization, we write where we have included a chemical potential (see below). The one-electron part contains the Coulomb and exchange contributions from the core orbitals, where ̃ ̃ ̃ ̃ is the anti-symmetrized electron repulsion integral (ERI) ((̃ ̃ |̃ ̃ ) − 1 2 (̃ ̃ |̃ ̃ ) in chemists' notation) and ̃ ̃ , ̃ ̃ correspond to the bare one-electron and Fock operator rotated into the EO basis respectively. The chemical potential in Eq. 3 appears only on the impurity orbitals (i.e. not on the bath orbitals) and is the same for all fragments. In the absence of the full self-consistency of DMET, it is needed to ensure the correct global number of electrons.
Given , we solve for the ground-state of each embedding Hamiltonian using a quantum chemistry solver. The ground-state solutions of each embedded impurity problem give embedding 1-RDMs Γ ̃ ̃ and two-body reduced density matrices (2-RDMs) Γ ̃ ̃ ̃ ̃ . We then evaluate the energy in a one-shot manner using the following recipe. We first assemble a correlated global density matrix as a weighted sum over the fragment density matrices Γ ̃ ̃ : The weight factor ̃ ̃ is 1 for̃ ,̃ ∈ imp, 0 for̃ ,̃ ∈ bath, and 1 2 otherwise; these factors are needed to avoid double counting of energy contributions from different impurities, using the democratic partitioning convention. From the global density matrix we can solve for the chemical potential such that tr Γ ′ = elec , where elec is the desired global number of electrons. We then update the one-electron part of the embedding Hamiltonian, and compute the total DMET energy as an average over all fragments, where the four-index weight ̃ ̃ ̃ ̃ is defined as 1 4 times the number of impurity labels in the set {̃ ,̃ ,̃ ,̃ }.

B. The DMET energy gradient
In a variational theory, the energy gradient depends only on the AO integral derivatives. However, from the above, we see that the DMET energy is not a variational expression, and thus there are additional contributions. Because the energy is not variational in the solver wavefunctions, the DMET gradient contains terms analogous to those appearing in the gradient of perturbation theories in quantum chemistry, i.e. the response of the impurity wavefunctions contributes. In addition, the DMET energy depends on a set of basis transformations. The AO→LO and LO→EO transformations are both atomic-position-dependent, and thus their response also contributes to the energy gradient. In particular, the LO→EO transformation is defined from the mean-field 1-RDM and as a result, the gradient requires the global MO response. Depending on the orthogonalization method, the AO→LO transformation can be a function of AO overlap (such as for Löwdin orthogonalization 20 ) and/or the mean-field solution or additional equations (such as for intrinsic atomic orbitals 21,22 ). In both cases, the AO→LO response depends on the AO overlap nuclear gradient.
Thus, in total, the DMET energy is a function of the set of impurity ground-state solutions = { } (these could denote, e.g. the configuration interaction (CI) vectors for a full CI (FCI) solver, or the amplitudes plus embedding MO coefficients for a coupled cluster (CC) solver), the global mean-field solution (MO coefficients) , and the AO integrals (the AO 1-electron and 2-electron integrals  and AO overlap integrals ). The gradient (with respect to nuclear positions or any perturbation) can be expressed in terms of partial derivatives with respect to these variables: Note that DMET does not depend on the chemical potential explicitly, but affects the energy indirectly through its effect on . The d ∕d and d ∕d contributions are related to d ∕d since and are implicit functions of as defined by the solution conditions (the mean-field and impurity ground-state equations) they satisfy. Written in this form, computing the DMET gradient naively requires one to solve for d ∕d and d ∕d for the Cartesian coordinates of each atom, i.e. 3 atom times. This complexity can be greatly reduced, however, by employing the so-called -vector technique, also known as the Lagrangian method 23,24 .
The DMET ground-state solutions involve three sets of working equations: the solution conditions for the impurity groundstates where each is the corresponding impurity solution condition (such as the CI, CC ground-state equations, etc.); the mean-field solution condition; and the chemical potential solution condition; We define a Lagrangian associated with the above constraints, and consider the stationary points w.r.t. all the variables. Vanishing derivatives of  w.r.t , and give the working equations of DMET, while vanishing derivatives of  w.r.t , and give equations for the minimizers of the multipliers, 0 , 0 , 0 respectively, If no chemical potential is used, Eq. 13 reduces to † 0 = DMET −1 (15) which can be solved for each impurity independently due to the block-diagonal nature of ∕ , i.e., . At the stationary points of ,  = DMET and d DMET ∕d = d∕d . However, the nuclear gradient of  no longer contains the contributions from the partial derivatives w.r.t or as in Eq. 8, as their implicit dependence on is accounted for by the multipliers. Thus, We note the repeated occurrence of DMET − † 0 in Eqs. 14 and 16, thus it is useful to define  emb = DMET − † 0 , then its derivative w.r.t. ∈ { , } is which takes the form of a contraction between density matrices and Hamiltonian derivatives, and where we have introduced a set of auxiliary RDMs, defined as ̃ ̃ , which account for all solverspecific details. In the Appendix, we give the corresponding expressions for the Hartree-Fock (HF) and FCI solvers as examples. Choosing to be the AO integrals  in Eq. 17 gives the "skeleton" derivative of the DMET gradient, while choosing to be the AO overlap leads to more complications, due to the orbital dependence on in localization, as discussed in Appendix VI D. These two components sum to the first two terms of Eq. 16. The remaining task is to solve for from  emb ∕ using the coupled-perturbed Hartree-Fock (CP-HF) equations 24,25 , and to consider the contributions from ∕ . Further details can be found in Appendix VI C.

III. RESULTS AND DISCUSSION
We first verify our derivation and implementation using some simple examples, then proceed to some prototypical reactions.
In all examples, we use Löwdin orthogonalized orbitals as our localized orbitals. We consider DMET calculations both with a chemical potential (DMET) and without a chemical potential (DMET (no )).

A. Numerical validation
Since DMET provides an exact embedding at the level of mean-field theory, the correctness of the DMET gradient using HF impurity solvers can be verified by comparing the DMET gradient to the gradient of a conventional HF calculation on the whole problem. The DMET gradients for correlated solvers can be validated by comparing to finite differences of the DMET energies, or by numerical integration of the analytic gradients to compare to the energies directly.

Water trimer
We take a water trimer as our first example (configuration shown in Fig. 1A). We used the 6-31G** basis set 26 and a HF impurity solver, with the DMET impurities defined as the individual atoms. The difference between the DMET energy and global HF energy is less than 10 −11 a.u., and the mean absolute error (MAE) in the gradient is less than 10 −8 a.u., demonstrating the correctness of the DMET energy implementation (recovering the exact mean-field embedding) as well as that of the gradient. We next test the FCI solver using a STO-3G basis 27 on the same system and similar atomic impurities. The deviation of the analytic gradient from the finite difference gradient is shown in Fig. 1A. We see that the finite difference and analytical gradients clearly agree, with the deviation between the two decreasing with the finite different displacement ( ).

Hydrogen ring
We next consider the symmetric dissociation of a hydrogen ring 6 . We treat the H 10 ring using a STO-3G basis, a FCI solver, and atomic impurities. We numerically integrate the DMET gradient along the dissociation profile to recover the PES. The correctness of the gradient is verified by the perfect match between the integrated PES and the directly calculated DMET energies (Fig. 1B).
Since the electron number in DMET (without ) changes only slightly along the symmetric stretching process (MAE = ∼0.02 with a standard deviation of ∼0.02) fitting here only slightly affects the energy.

B. Diels-Alder reaction
We next use the prototypical Diels-Alder reaction 28 to illustrate using the DMET gradient to optimize structures and locate reaction transition states. Although the reaction in this case takes place at a single "site" (making it easy to define specific atoms and orbitals involved in the reaction) we still treat it computationally as a multi-site problem, breaking the problem into multiple fragments in the calculation. The specific Diels-Alder reaction we choose is the addition of methyl vinyl ketone (MVK) to cyclopentadiene (CP) (see Fig. 2), which has previously been used to test a QM/MM transition-state searching protocol 29 . We used the 6-31G basis set 27 for all atoms and each hydrogen was defined as an impurity while for heavy atoms, the atomic orbitals were split into two sets constituting two impurities: the first containing the 1s 2s 2p orbitals, and the second, the 3s 3p orbitals. All impurity correlation problems were treated using FCI. Since DMET does not capture long-range fluctuations such as those giving rise to dispersion, we also investigated adding an empirical D3 30 dispersion correction (parameterized for HF) to the energy. The reactant and product states were then optimized using the DMET-D3 gradient (including the fitted chemical potential ) and the gradient was further used in a nudged elastic band (NEB) [31][32][33] calculation to find the reaction transition state. The obtained transition state structure was verified via vibrational frequency analysis using the numerical Hessian computed with the DMET-D3 gradient. All the optimization was performed using ASE 34 . Fig. 2 shows the DMET-D3 optimized structures and associated energies (see also Tab. I). In this case, DMET without fitting gives a large error in the total electron number (0.734 for the reactants, 0.704 at the transition state, and 0.706 for the products), thus we report only DMET results which fit the chemical potential. We use CC singles and doubles (CCSD) with renormalized triples, namely R-CCSD(T) 35 , as the reference energy value; we note that adding the triples correction results in a relatively large shift in the reaction barrier, highlighting the open-shell character of the transition state. After adding the empirical D3 correction, the DMET reaction barrier is similar to the R-CCSD(T) value, while DMET without D3 yields a result comparable to the CCSD barrier. Regardless of whether the dispersion correction is used, DMET yields a better reaction barrier compared to its underlying mean-field approach (HF), and compared to a commonly used density functional theory (DFT) approach (M06-2X 36 ) (although the latter comparison is complicated by the different basis set dependence). In terms of the relative stability of the reactants and products, the DMET results show worse agreement with the reference R-CCSD(T) energies. This shows that using very small fragments is unfortunately not always sufficient for quantitative energetics.
In addition to the energies, a frequency analysis at the CCSD level at the DMET-D3 transition-state shows that the DMET-D3 transition state is a good approximation to the transition state of CCSD, as there is only a single large CCSD imaginary frequency (larger than 1 cm −1 ) at the geometry. The imaginary frequency obtained with the DMET-D3 numerical Hessian (562 cm −1 ) is in good agreement with the imaginary frequency obtained with CCSD (581 cm −1 ).

C. Proton transport in a water cluster
We finally consider proton transport (PT) in water as an example of a disordered chemical reaction where the reactions can involve multiple components and multiple locations, making it difficult to isolate a single reaction center to be treated by highlevel quantum methods. In PT, the multi-site character arises due to the Grotthuss proton hopping mechanism 37 where an excess proton is transported by forming and breaking its surrounding bonds with the solvent water molecules, with only minimal physical movement of the atoms. As a result, the identity of the excess proton changes dynamically and the associated reaction center cannot be tied to any specific proton or water molecule. Fig. 3 shows an example of such a process where one excess proton is solvated and transported in a water tetramer. The trajectory was generated by MD simulations at 50K using the DMET gradient (without fitting the chemical potential, to reduce cost) using a STO-3G basis for oxygen and the 6-31G basis for hydrogen, with each atom defined as an impurity to be treated by the FCI solver (the small oxygen basis was chosen to allow for a tractable FCI cost). The reference energy curve was generated at the CCSD(T) level using MD sampled configurations. Similar to the hydrogen ring system, the electron number fluctuated in a small range that is close to the correct number (MAE = ∼0.06 with standard deviation ∼0.002) thus chemical potential fitting is seen to change the energy only slightly. We used parallel-biased metadynamics 40 to enhance the sampling of PT and the morphology of the cluster via biasing the 10 pairwise distances between the four water oxygens and a virtual atom that denotes the excess proton position 41,42 . Along these biased trajectories, the excess proton constantly changes its identity between the protons shared by the two waters on the right and those shared by the two waters on top, similar to processes seen in aqueous solution and in bio-molecules 31,43-48 making the definition of a fixed reaction center almost impossible. In such scenarios, if one wished to use a single reaction center then one might seek to dynamically define it to follow where the reaction takes place. For example, we might change the reaction center from the right H 5 O 2 + complex (i.e., whole system is partitioned vertically) to the upper H 5 O 2 + complex (i.e., partitioned horizontally) when the proton hops. As shown in Fig. 3, when two different partitioning schemes are used corresponding to the above two choices, one commonly used quantum embedding method, the projection-based wavefunction-in-DFT embedding method 2 provides an accurate description in the corresponding reaction scenarios respectively. However, relatively large errors are observed when the active excess proton goes beyond the defined fragment. Thus although each DFT embedding curve is continuous, one would have to accept a discontinuous energy jump between two different PES when changing the fragment definition to follow the excess proton identity. Another segment of the MD trajectory (Fig. 4) shows more clearly the difficulty of defining an appropriate reaction center. The excess proton remains as a shared proton between the upper two waters in the period of 30.7 ps to 32.3 ps, suggesting a horizontal partition is appropriate, while the hydrogen bond between the two lower waters is broken. We see a notable discrepancy between the DFT embedding curve and the whole-system CCSD(T) potential energy curve in the region around 30.7 ps. This is likely due to the poor description of the hydrogen bonds in the underlying mean-field (PBE) treatment in the DFT embedding. To improve this, we would then need to describe the whole problem using a high-level correlated quantum treatment, thus removing any cost savings from performing the "quantum embedding" in the DFT embedding approach.
In contrast, DMET produces a continuous PES that closely follows the reference CCSD(T) curve for both segments shown in Figs. 3, 4, and shows a significant improvement upon the underlying mean-field (HF). This illustrates the power of the DMET bath construction which finds the most important orbitals to correlate regardless of the position of the reaction center, and demonstrates the potential of DMET as a minimal framework for multi-site reactive dynamics. An important caveat is that the quantitative accuracy of the DMET energy throughout the reaction trajectory still relies on the correlations being reasonably local. However, this appears to be the case in this simple system.

IV. CONCLUSION
In this work we derived and implemented the analytic nuclear gradient of the density matrix embedding theory (DMET) and validated our implementation via several numerical examples. In particular, we demonstrated the use of the gradient to perform geometry optimization, transition state searches and reactive molecular dynamics simulations. In the various scenarios, we used  Fig. 3. The DMET, HF and DFT embedding/vertical trajectories are aligned to CCSD(T) using the averaged energy in the first 2.5 ps segment. The DFT embedding/horizontal trajectory is aligned using the last 2.5 ps segment average. DMET with small embedded fragments (in some cases, just a group of atomic orbitals) and observed that we could qualitatively reproduce high-level correlated calculations carried out on the full system. Importantly, the DMET potential energy surface remains continuous and thus allows for qualitatively correct reaction dynamics even for problems without a well-defined reaction center, as we demonstrated in a small water-based proton transport simulation.
Although the DMET calculations are qualitatively reasonable across the reaction trajectory, we did not always obtain results of quantitative accuracy, due to the use of small fragments, as seen when modeling a Diels-Alder reaction. For molecular dynamics, however, small embedded fragments are desired to limit the cost of a time-step. Thus the development of empirical corrections to account for the missing correlations in DMET (perhaps along the lines of dispersion corrections) are needed. For other applications, such as single-point determinations of reaction barrier heights, more expensive impurity calculations are potentially feasible. The generalization of the gradient theory described in this work to such extensions of density matrix embedding is then an interesting future direction.

V. ACKNOWLEDGEMENTS
This work, including the development of embedding gradient theory, was primarily supported by the US Department of Energy, Office of Science, via Award No. DOE-SC0023318. Additional work on the PySCF backend was supported by the US Department of Energy, Office of Science, via Award No. DOE-SC0019390. GKC is a Simons Investigator. The authors thank Dr. Zhihao Cui for helpful discussions.

VI. APPENDIX
In the following derivations, we will intentionally use the component form, rather than the matrix form, of the equations where appropriate to avoid any confusion about basis and representation. The basis transformation matrices will be identified by their indexes, e.g. ̃ is interpreted as the AO→EO transformation such that |̃ ⟩ = | ⟩ ̃ . We will also assume the orbital coefficients are real, although the equations can be straightforwardly generalized to complex-valued orbitals.
A. DMET gradient terms arising from the impurity solution The DMET energy gradient w.r.t the solution of the impurity problem is needed to obtain in Eq. 13. Consider , the impurity solution of fragment . Then differentiating Eq. 7 gives where we have a sum over in the last term because ℎ ′ depends on the solutions for all fragments (see Eq. 6), , via its dependence on the updated global 1-RDM Γ ′ . Explicitly, (no summation over ) where we have summed over first in Γ ̃ ̃ ̃ ̃ ̃ ̃ , used Eq. 5, and performed some dummy variable changes, and defined the effective potential ′ = Γ ′ . Inserting back into Eq. 18, (no summation over ) where the terms in the square brackets define ℎ ′′ .
In the case of no chemical potential fitting, is directly obtained after a contraction with ( ∕ ) −1 (where the contraction is performed by solving the response equation of the underlying solver method, such as coupled perturbed CI (CP-CI) 24 for a CI solver, or the Λ-equations for CC solvers 49 ). If chemical potential fitting is performed, the solver responses of different impurities are connected via a common chemical potential. Thus their 's are solved for at the same time as , the multiplier associated with the electron number constraint (Eq. 13).

Hartree-Fock solver response
The Hartree-Fock impurity solver response enters into the DMET gradient theory in two ways: (1) it is the starting point for the response theory of more sophisticated post-Hartree-Fock solvers such as CC, truncated CI and complete active space (CAS)-CI, (2) it enters directly if Hartree-Fock is used to treat some of the impurities while others are treated by correlated methods. The DMET energy gradient w.r.t the HF solution of the impurity problem (Eq. 20) is (we consider ( DMET ∕ ) instead of DMET ∕ to connect to the well-known CP-HF equation) In the above, we usẽ and̃ to label the Hartree-Fock impurity orbitals, i.e. the MOs in the embedding space (̃ ,̃ will be used for occupied and̃ ,̃ for virtual orbitals) whilẽ ,̃ ,̃ and̃ , consistent with the main text, denote embedding orbitals, i.e., the computational basis in the embedding space. The embedding MO coefficient matrix is thus written as ̃ ̃ or ̃ ̃ , and the embedding "Fock" matrix is defined as The CP-HF equation is usually written as where  and in our language are  = − and = d d with ̃ ̃ = ̃ ̃ = 0 the Hartree-Fock equation of the embedding problem. The solution is related to MO response via −1 (d ∕d ). Comparing to Eq. 15, if no -fit is performed, for the HF solver is obtained from the CP-HF equation by replacing with ′ , i.e., The auxiliary RDMs for the HF solver (defined after Eq. 17) are then obtained by contracting with ℎ̃ ̃ and (̃ ̃ |̃ ̃ ) If a -fitting is performed, the above auxiliary RDM expressions still hold but should be obtained from Eq. 13 using all the impurity response equations simultaneously. The other information needed in Eq. 13 is and

FCI solver response
The FCI vector  is the solution of the equations, where we use  and  to label an electronic configuration (e.g. a determinant). The RDMs are constructed from  via where Note that since the FCI solution only depends on the space of , we can perform the CI calculation directly in the EO basis without the need to perform an embedded HF calculation (to determine MOs). Then Eq. 20 becomes which has the form of a CI Hamiltonian contracted with the CI vector 2  to give a new CI vector. Note that if there are no weights and if ℎ instead of ℎ ′ is used to compute the energy, Eq. 32 simply becomes 2   = 2 CI  , which vanishes when dotted with d  ∕d due to the orthogonality between and d ∕d . It is obvious that such a Hellmann-Feynman condition is not fulfilled, and this is why the CI response is needed for the DMET gradient even though CI is a variational method. Differentiating the CI equation (Eq. 29) w.r.t  and after some linear algebra, we have In the case of no -fitting, we obtain i.e., is the solution of a CP-CI equation with  . Given that the auxiliary RDMs defined in Eq. 17 become To solve for with -fitting from Eq. 13, we also need where ̃ ∈ = 1 if and only if̃ is occupied in configuration , and

B. Orbital transformation gradient
The DMET gradient includes contributions from the response of the orbital transformations. A position change induces both AO→LO and LO→EO changes that contribute to a change in the energy. This response is encoded in the derivative of  emb w.r.t the AO→LO and LO→EO coefficients, while both of them, via the chain rule, are related to the derivative  emb ∕ ̃ . Following our conventions for basis indexing, ̃ = ̃ is the transformation matrix from AO to EO, and represents the AO→LO and ̃ represents the LO→EO transformation as defined in Eq. 2.
We write the three components of  emb ∕ ̃ appearing in Eq. 17 as C. DMET gradient contribution from the global mean-field The DMET energy response w.r.t the mean-field solution (  emb ∕ ) is required to solve for (Eq. 14). We label the global occupied MOs by , and the virtuals by , . satisfies = 0 (44) Similar to the treatment of the HF solver in section VI A 1, it is more convenient to consider  emb ∕ def ==  than  emb ∕ . Assuming we use Löwdin orthogonalization to obtain the LOs, then there is no mean-field contribution to the AO→LO transformation gradient. As such, ℎ ′̃ ̃ and (̃ ̃ |̃ ̃ ) depend on only through the LO→EO transformation ̃ , which is obtained from the SVD of . On the other hand, ℎ̃ ̃ not only depends on ̃ due to the basis transformation but also has the extra "skeleton" dependency on , due to its explicit appearance in the global Fock matrix used to define ℎ̃ ̃ (Eq. 4).
We first consider the "non-skeleton" part of  , arising from the LO→EO transformation, expressed as a chain of partial derivatives, The and indices in Eq. 45 are restricted because the bath is constructed from the env-imp block of . The and̃ indices are restricted to be in the environment and bath since the other elements of the LO→EO transformation matrix are 0 or 1 (Eq. 2) and do not contribute to the derivative. The partial derivative ̃ ∕ corresponds to the standard expression for the derivative of the SVD vectors w.r.t the SVD input matrix, and the prefactor  emb ∕ ̃ in Eq. 45 is given in Appendix VI B. As discussed earlier, the "skeleton" part of  has the contribution from ℎ only: Adding Eq. 45 and Eq. 46 completes the DMET "Fock" matrix  , and when combined with the CP-HF equation, we obtain the third term of the DMET gradient (Eq. 16) as where [ ] is the "skeleton" nuclear gradient of the MO overlap matrix, i.e., , and the expressions for  and can be found in the CP-HF literature, such as Ref. 24 .

D. DMET gradient arising from the AO overlap and AO integrals
We consider the first two terms of the DMET gradient (Eq. 16) to be composed of two contributions: The first term comes from computing the embedding Lagrangian using AO integral derivatives ( and ( | ) ) in place of the original integrals: Regarding the second term in Eq. 48, the Lagrangian  emb depends on in two ways: (1) from the explicit appearance of in ℎ̃ ̃ when transforming the AO 1-RDM into the EO basis using ̃ = ̃ , (2) and the implicit dependency on of the LOs. The first contribution, which can be regarded as a "skeleton" derivative, is simply