# Vision Transformer Compression with Structured Pruning and Low Rank
Approximation

Ankur Kumar
Department of Computer Science
University of California, Los Angeles
<EMAIL_ADDRESS>

###### Abstract

Transformer architecture has gained popularity due to its ability to scale
with large dataset. Consequently, there is a need to reduce the model size and
latency, especially for on-device deployment. We focus on vision transformer
proposed for image recognition task (Dosovitskiy et al.,, 2021), and explore
the application of different compression techniques such as low rank
approximation and pruning for this purpose. Specifically, we investigate a
structured pruning method proposed recently in Zhu et al., (2021) and find
that mostly feedforward blocks are pruned with this approach, that too, with
severe degradation in accuracy. We propose a hybrid compression approach to
mitigate this where we compress the attention blocks using low rank
approximation and use the previously mentioned pruning with a lower rate for
feedforward blocks in each transformer layer. Our technique results in 50%
compression with 14% relative increase in classification error whereas we
obtain 44% compression with 20% relative increase in error when only pruning
is applied. We propose further enhancements to bridge the accuracy gap but
leave it as a future work.

##  1 Introduction

BERT (Devlin et al.,, 2019) and GPT (Radford and Narasimhan,, 2018) inspired a
trend of pre-training huge neural networks on large scale dataset for very
accurate results on downstream tasks. Large pre-training corpus and a model
that scales well on such dataset are critical in this process. Kaplan et al.,
(2020) and Henighan et al., (2020) show that transformer architecture (Vaswani
et al.,, 2017) strikes a balance between computation efficiency and
scalability. Consequently, training transformer models with several billion
parameters has become a trend to achieve high accuracy (Ramesh et al.,, 2021;
Radford et al.,, 2021). Vision transformer (Dosovitskiy et al.,, 2021) show
that large scale pre-training of transformers can benefit computer vision
tasks also. The problem, however, is the resource constrained environment for
deployment, specially for on device settings.

Model compression is a well researched topic. Several techniques such as
knowledge distillation, pruning, low rank approximation and quantization have
been studied extensively to reduce model size without significantly degrading
accuracy. However, few works have applied these compression techniques to
vision transformer (Zhu et al.,, 2021; Yu and Wu,, 2021; Yang et al.,, 2021;
Chen et al.,, 2021; Hou and Kung,, 2022). Most of these works propose
variations of structured pruning, which does not require specialize hardware
to run the pruned model as opposed to unstructured pruning.

In this work, we follow Zhu et al., (2021) because of their simplistic and
structured pruning approach. However, we question their importance score
estimation for pruning groups of neurons, and analyze the pruned network for
this purpose. We conduct experiments on CIFAR-10 dataset to provide empirical
evidence that the pruning strategy may be indeed suboptimal. In this process,
we come up with a hybrid compression technique where we apply low rank
approximation in combination with the above pruning techniqe to mitigate some
of the accuracy degradation. We further discuss that advanced techniques for
pruning (Molchanov et al.,, 2019) and low rank approximation (Chen,, 2021) can
be used in our _modular_ framework to further improve model performance.

Next, we present a brief overview of prior works on model compression,
including compressing transformer models, in section 2. We present our method
in section 3. Section 4 and 5 contain details of our experiments. Finally, we
discuss shortcomings of the pruning technique mentioned earlier with possible
future directions in section 6, and conclude in section 7.

##  2 Related Work

###  2.1 Knowledge Distillation

Hinton et al., (2015) proposed to train a smaller neural network by first
training a much bigger model and then training the smaller one using soft
probabilities of the previous one. Since then, this technique has been
successfully applied in multiple applications (Gou et al.,, 2021). However, we
leave out this type of compression technique since it involves training the
_teacher_ model first and then running it again to train a _student_ model,
which may require significant computational resources.

###  2.2 Quantization

A simple way to compress a model is to use low precision for weights and
activations. With 8-bit quantization, model size can be reduced by a factor of
4 with negligible accuracy loss (Krishnamoorthi,, 2018; Nagel et al.,, 2021).
Multiple efforts try to further reduce the precision, train a model with low
precision etc (Gupta et al.,, 2015). However, it is difficult to lower the
precision beyond 16-bit without significant accuracy degradation and
specialized hardware (Wang et al.,, 2018). It has become a trend to use half
precision for mixed precision training to avoid the accuracy drop and make use
of the available hardware (Micikevicius et al.,, 2018).

###  2.3 Pruning

Over-parameterization of neural networks is considered to be one of the
reasons behind its superior performance. Many works have also established that
a lot of _unimportant_ weights can be removed to compress the models without
sacrificing accuracy (Han et al.,, 2015; Blalock et al.,, 2020). There are
several criteria to decide the importance score (LeCun et al.,, 1989; Hassibi
and Stork,, 1992; Molchanov et al.,, 2019). A common and convenient way is to
zero out the weights with small magnitude. This is a type of _unstructured_
pruning. However, simply zeroing out does not improve the model size and speed
unless there is specialized hardware. Another method, called _structured_
pruning (Li et al.,, 2017; Wen et al.,, 2016), avoids this issue by dropping a
group of weights such as entire convolution filter or attention head.
Therefore, we consider structured pruning in this work and combine it with
other compression techniques (Han et al.,, 2016).

###  2.4 Low Rank Approximation

Decomposition of weight matrices into low rank matrices is yet another
technique to reduce the model size without affecting its structure. In case of
convolution, factorizing channel and spatial convolution weights is used in
place of traditional convolutional operation due to lesser number of
parameters. Feedforward weight matrices can be factorized into two lower rank
matrices to reduce the number of parameters (Sainath et al.,, 2013). However,
it is not necessary that weight matrices are low rank, which we confirm for
vision transformer from our experiments. The low rank structure can still be
used if input feature is believed to lie in smaller dimensions (Chen,, 2021).

###  2.5 Transformer Compression

With the success of BERT, several works have studied how to make transformer
models smaller. DistilBERT (Sanh et al.,, 2019) is a widely used compressed
version of BERT models. Another line of work tries to make transformer models
more efficient since computational complexity of transformer model grows
quadratically in the number of input tokens. This is challenging for natural
language processing, specially speech processing. Wang et al., (2020) show
that attention matrix is typically low rank, and therefore, can benefit from
low rank approximation.

Most of the works in computer visoin domain focus on removing uninformative
patch tokens to speed up the training and inference steps (Wang et al.,,
2021). Vision transformer pruning (VTP) (Zhu et al.,, 2021) removes
unimportant dimensions (columns or rows) of matrices in a transformer block.
Inspired by several works combining multiple compression strategies (Han et
al.,, 2016; Mao et al.,, 2020), we build upon VTP and first show the
shortcomings of this technique. Then we propose to use low rank approximation
in combination with VTP to mitigate some of the issues. Our work is orthogonal
to token pruning, quantization and knowledge distillation and can be combined
with these approaches.

##  3 Method

###  3.1 Vision Transformer

An input image $x\in\mathbb{R}^{H\times W\times C}$ is divided and reshaped to
have $P\times P$ patches, each having $P^{2}C$ pixels. The reshaped input
$x_{P}$ has shape $\mathbb{R}^{N\times P^{2}C}$ where $N=HW/P^{2}$. Each patch
is linearly projected into a $d$ dimensional space using a learnable embedding
matrix. Following BERT, a learnable embedding corresponding to special token
$[class]$ is prepended to $x_{P}$ making total sequence length $N+1$. This 1-d
sequence of $d$ dimensional token can be fed to standard transformer model.
Dosovitskiy et al., (2021) use learnable 1-d positional encoding which is
added to the token embeddings. Structure of a transformer block is as shown in
figure 1(a). We do not show GELU non-linearity which is applied after each of
the top two linear layers in the block.

Output embedding corresponding to the $[class]$ token after the final
transformer block serves as the input image representation. A classification
head is attached on top of this to train the model. The head may contain
additional hidden layer.

(a) Transformer block

(b) VTP (norm layer not shown)

Figure 1: Transformer block with and without pruning layers. GELU non-
linearlity is not shown in the feedforward block (consisting of top two linear
layers) but applied on top of the linear layer outputs. We do not show the
normalization layer in pruning case for simplicity.

###  3.2 Vision Transformer Pruning

Each transformer block has almost all of the parameters in weight matrices of
linear layer and query, key, value of self-attention layers. Vision
transformer pruning Zhu et al., (2021) aims to drop unimportant columns or
rows from a matrix. To achieve this, dimension pruning layers are inserted
between two layers containing weight matrices as shown in figure 1(b). The
pruning layer is a binary mask with dimension equal to the input dimension of
next layer (or output dimension of previous layer) and is multiplied with the
output of previous layer. Therefore, if $i^{th}$ entry in the mask is 0, then
it is equivalent to dropping $i^{th}$ column from the next layer and $i^{th}$
row from the previous layer. If the previous layer happens to be a self-
attention layer, then $i^{th}$ row is dropped from each of the key, query and
value matrices.

Since binary mask is non-differentiable, the above method relaxes mask to be a
continuous value during training. More precisely, all mask values in a
dimension pruning layer is initialized to 1 indicating that initially no
column or row is dropped at the beginning. The following loss is optimized to
learn the weights and mask

$\displaystyle L_{CE}+\lambda\sum_{l=1}^{L}|m_{l}|_{L1}$

where $L_{CE}$ is the regular cross entropy loss for classification task and
$|m_{l}|_{L1}$ is the L1 norm of the mask or dimenison pruning layer in
$l^{th}$ transformer block. Intuitively, the L1 loss promotes sparsity making
some of the mask values close to 0, which results in lower contribution of
corresponding dimensions in the overall activation. Therefore, we can drop the
entire column to obtain less parameters. During inference, we again discretize
the mask as follows

$\displaystyle m_{l}=m_{l}>\tau$

where $\tau$ is the threshold determined from desired pruning rate.

###  3.3 Hybrid Compression Framework

Figure 2(a) shows the effect of applying VTP to vision transformer. We observe
that mainly matrices from feedforward modules (top two linear layers in a
transformer block) are pruned whereas attention weight matrices and projection
layers are majorly unpruned. We discuss the reasons in more detail in section
6. However, this motivates us to compress self-attention block (query, key,
value matrices in scaled dot-product layer and the following linear layer) in
combination with above pruning approach for extra compression. We apply low
rank approximation (LRA) for this purpose.

(a)

(b)

Figure 2: (a) Effect of pruning with VTP in vision transformer layers.
Feedforward layers are pruned and self-attention layers remain unaffected (b)
Hybrid compression framework. Layer norm is not shown but applied before
projection matrix and first linear layer in feedforward block

One issue with low rank approximation is that new rank $k$ must be much
smaller than original rank $d$ since reducing the rank from $d$ to $k=d/2$
does not lead to lesser parameters if the original matrix has shape $d\times
d$. We partly avoid this problem by applying LRA only to the attention block,
where we share one of the low rank matrices among key, query and value
matrices. It can be implemented by projecting $d-$dimensional input of a
transformer block to a lower dimensional space of size $k$ by inserting a
linear layer or projection matrix as shown by blue layer in figure 2(b)(b).
This changes input dimension of query, key and value matrices from $d$ to $k$.
We can now compute attention weights in either $d-$dimensional space or a
lower $k-$dimensional space. We refer to the former as LRA($d\times k$) and
the later as LRA($k\times k$). LRA($d\times k$) does not affect the size of
linear layer after self-attention layer. However, LRA($d\times k$) reduces the
input dimension of this layer also to $k$.

We apply LRA only to the attention block and not the feedforward block. This
is because by design, transformer feedforward block is a factorized matrix
where hidden dimension in between them is kept higher. We experiment with only
LRA applied to vision transformer. We also apply VTP discussed earlier in the
feedforward block (as shown in dashed green layer in figure 2(b)(b)) with LRA
in attention block. This gives our hybrid compression framework.

##  4 Experiments

We conduct experiments on CIFAR-10 dataset which has 10 output classes, 50,000
training images and 10,000 test images. We do not use any data augmentation.
All the models are trained on a single gpu with a batch size of 64 and Adam
optimizer with $10^{-4}$ learning rate. We also apply dropout of 0.1 to input
feature embedding (after adding positional embedding) and all the weight
matrices in the network except attention key, query and value matrices. Unless
mentioned otherwise, we train for 15 epochs due to resource constraint.
However, we also train a model for 50 epochs to validate our hypothesis which
we discuss in section 6. Further, a pruned network is normally finetuned for
few epochs to recover any accuracy drop due to pruning. However, we do not do
finetuning because we believe that our model is not fully trained. As a
result, it may be difficult to credit improvement in performance to either
finetuning process or the fact that models are trained for extra epochs while
comparison with the unpruned model.

##  5 Results

###  5.1 Vision Transformer Training

We conduct hyper-parameter tuning for vision transformer for which results are
shown in table 1. We experiment with patch size $P$ of 2, 4 and 8. Results are
shown in the first 3 rows. Patch size of 8 performs the worst and patch size
of 2 performs the best. However, we choose patch size of 4 because of better
training efficiency since transformer computational complexity grows
quadratically with number of input tokens.

Then we varied the attention dimension and hidden dimension of feedforward
blocks. Attention dimension of 512 and feedforward hidden dimension of 1024
worked the best as given in $6^{th}$ row. Finally, we changed the number of
transformer blocks. We did not find significant accuracy difference between 6
and 9 layers. Therefore, we chose 6 transformer blocks in our network for
better computation time. All the models compute self-attention with 8 heads
(512/8 = 64 dimension for each head). Our final choice of hyper-parameter is
shown in the second last row of table 1. It achieves 72.61% accuracy with 13M
parameters on CIFAR-10 test set.

Patch Size | Dim | MLP Dim | Layers | Accuracy | Parameters
---|---|---|---|---|---
2 (16x16) | 512 | 512 | 6 | 72.81% | 10M
4 (8x8) | 512 | 512 | 6 | 70.68% | 10M
8 (4x4) | 512 | 512 | 6 | 63.93% | 10M
4 (8x8) | 256 | 512 | 6 | 66.91% | 3.5M
4 (8x8) | 512 | 512 | 6 | 70.68% | 10M
4 (8x8) | 512 | 1024 | 6 | 72.61% | 13M
4 (8x8) | 512 | 1536 | 6 | 72.06% | 16.6M
4 (8x8) | 512 | 1024 | 3 | 68.64% | 7M
4 (8x8) | 512 | 1024 | 6 | 72.61% | 13M
4 (8x8) | 512 | 1024 | 9 | 73.03% | 19.5M
Table 1: Hyper-parameter tuning for vision transformer trained on CIFAR-10
dataset

###  5.2 Vision Transformer Pruning

For pruning using VTP, we again train vision transformer with the final choice
of hyperparameters. We insert the dimension pruning layers in each transformer
block as shown in figure 1(a)(b). We have only one hyper-parameter $\lambda$,
which measures the strength of the sparsity or L1 loss. We experimented with
values of $10^{-2},10^{-3}\text{ and }10^{-4}$. We found that higher values of
$\lambda$ significantly degraded the model performance but $10^{-4}$ worked
well with a very small drop in model accuracy. This value is also the same as
that in the original work although for ImageNet experiments. We show the
results for VTP with different pruning rate in table 2. Prune rate of $\alpha$
means we want to zero out $\alpha$ fraction of all the mask values combined
across all the pruning layers inserted in our model. Note that we have to re-
train only one model with pruning layers as pruning rate is specified only at
the inference time.

Model | Test Accuracy | Increase in Error | Model Size | Compression
---|---|---|---|---
Base | 72.26% | - | 13.2M (51MB) | -
VTP (prune rate = 0.1) | 68.4% | 14% | 12.3M (47MB) | 7%
VTP (prune rate = 0.3) | 69.4% | 10% | 10.4M (40MB) | 21%
VTP (prune rate = 0.5) | 69.06% | 11.5% | 8.5M (33MB) | 36%
VTP (prune rate = 0.6) | 66.57% | 20.5% | 7.4M (29MB) | 44%
Table 2: Vision Transformer Pruning

We can see from the second last row that we can compress the model by 36% with
around 12% increase in error. However, the error goes up significantly if we
compress the model further. We make two observations in this table. First,
prune rate of 0.1 has higher error than prune rate of 0.3 which is counter-
intuitive. We discuss this in detail in section 6. Second, pruning rate is not
exactly the same as fractions of weights pruned because pruning rate is
specified for fraction of dimensions to be pruned and not weights. We also
observed that prune rate greater than 0.6 gave error perhaps because all
dimensions of some weight matrix was getting dropped. This is indeed possible
with VTP whereas low rank approximation in our method systematically
compresses all attention blocks.

###  5.3 Hybrid Compression

We chose $k=d/4$ for low rank compression. The results for LRA are shown in
table 3. We note that we train a new model for all the LRA experiments (with
and without pruning). We observe that only applying LRA is better than VTP for
similar compression. For example, LRA $(d\times k)$ compresses the model by
24% and has 5.8% increase in error whereas pruning has 10% increase in error
for similar compression. Similarly, LRA $(k\times k)$ has 39% compression with
8% increase in error whereas pruning has 11.5% increase in error. When we use
pruning in feedforward block and LRA in attention block (hybrid compression),
we achieve 50% compression with 14% increase in error whereas pruning by
itself has 20.5% increase in error for only 44% compression, which is about 2%
absolute accuracy gap.

Model | Test Accuracy | Increase in Error | Model Size | Compression
---|---|---|---|---
Base | 72.26% | - | 13.2M (51MB) | -
VTP | 69.06% | 11.5% | 8.5M (33MB) | 36%
VTP | 66.57% | 20.5% | 7.4M (29MB) | 44%
LRA ($d$ x $k$) | 70.65% | 5.8% | 10M (39MB) | 24%
LRA ($k$ x $k$) | 70.03% | 8% | 8M (31MB) | 39%
VTP (FF only) + LRA ($d$ x $k$) | 68.46% | 14% | 6.6M (25MB) | 50%
Table 3: Hybrid Compression for Vision Transformer

##  6 Discussion

In this section, we discuss multiple shortcomings of VTP approach. The first
issue, mentioned in section 3.3, is the type of layers getting pruned.
Specifically, mostly dimensions for feedforward weight matrices are dropped
and attention blocks are relativley unpruned as shown in figure 2(b). This
means that mask values for attention blocks would be much higher as compared
to mask values for feedforward block since smaller mask values are pruned. We
plot the histogram of mask values for attention blocks and feedforward blocks
separately as show in figure 3. It shows that, indeed, most attention mask
values are greater than 0.8 whereas feedforward mask values roughly follow a
normal distribution around 0.5.

(a)

(b)

Figure 3: Histogram of all mask values in attention and feedforward blocks
after training vision transformer with pruning

The second issue as mentioned in section 5.2 is VTP having lower accuracy with
prune rate of 0.1 as compared to prune rate of 0.3. One reason for this could
be that mask values are not close to 0 after training for 15 epochs. However,
we hypothesize that mask values after training may not actually reflect the
importance of corresponding dimensions. So, while pruning, we may be removing
an important feature dimension whose absolute mask value is the lowest,
resulting in significant drop in accuracy (72.26% to 68.4%). In addition, if
the dimensions being pruned were indeed of lowest importance, the accuracy
should not have improved to 69.4% on increasing the pruning rate to 0.3.

We verify our hypothesis that importance of a dimension may not be correlated
with importance of corresponding mask values by training a model for 50
epochs. We expect many mask values to be close to 0 this time due to L1 loss
and longer training. The distribution of mask values are shown in figure 4.

(a)

(b)

Figure 4: Histogram of all mask values in attention and feedforward blocks
after training vision transformer with pruning for 50 epochs

We find that mask values in attention blocks are not affected much. But many
mask values for feedforward blocks indeed become 0 or close to 0. Therefore,
we would expect that VTP should perform well for this case. It turns out that
the accuracy drop here is even more severe with at least 48% increase in error
as shown in table 4. We still observe the irregular pattern that accuracy
increases and then decreases upon increasing the pruning rate from 0.1 to 0.6.
This is because of hard thresholding to get binary mask values during
inference as compared to soft mask values during training. For a small prune
rate such as 0.1 or less, we may be converting many mask values, which are
close to 0, to 1 resulting in a severe degradation in accuracy.

Model | Test Accuracy | Increase in Error | Model Size | Compression
---|---|---|---|---
Base | 80.71% | - | 13.2M (51MB) | -
VTP (prune rate = 0.1) | 56.11% | 127.5% | 12.3M (47MB) | 7%
VTP (prune rate = 0.3) | 62.48% | 94.5% | 10.4M (40MB) | 21%
VTP (prune rate = 0.5) | 71.45% | 48% | 8.5M (33MB) | 36%
VTP (prune rate = 0.6) | 70.26% | 54.2% | 7.4M (29MB) | 44%
Table 4: Vision Transformer Pruning (50 epochs)

This brings to one of the future directions. We would like to augment or
replace VTP with pruning tecniques that actually consider the importance of a
group of weights (dimension in VTP is one type of grouping) for final
accuracy. For example, Taylor pruning (Molchanov et al.,, 2019) estimates the
importance using first order Taylor approximation of error when a group of
weights is dropped. It also removes the need to threshold mask values during
inference due to which we observed severe degradation in accuracy for VTP.

Second direction that we would like to explore is improving the low rank
approximation. We observe more than 2% absolute drop in accuracy when we use
LRA (row 4 and 5 in table 3). This is a significant drop for compression of
40% and a model which is far from the best accuracy that can be achieved.
Chen, (2021) propose a data-driven approach for low rank approximation,
although in natural language processing, and show that the method is provably
optimal. We can benefit from the modularity of our compression method and
apply both the enhancements for better model performance. We get the extra
benefit that Taylor pruning and data-driven LRA can be readily applied to any
pre-trained models, thereby eliminating the need to re-train a model from
scratch as in current scenario.

##  7 Conclusion

We presented a hybrid compression framework for vision transformer using
structured pruning and low rank approximation. The proposed LRA has better
compression and accuracy tradeoff as compared to VTP. The combined or hybrid
approach had even better compression with higher accuracy than VTP. We also
discussed a major shortcoming of VTP approach and validated it on CIFAR-10
dataset. Finally, we proposed various ways in which the hybrid compression
technique can be further improved using advanced pruning and low rank
approximation techniques.

## References

  * Blalock et al., (2020) Blalock, D. W., Ortiz, J. J. G., Frankle, J., and Guttag, J. V. (2020).  What is the state of neural network pruning?  ArXiv, abs/2003.03033.
  * Chen, (2021) Chen, P. H. (2021).  Drone: Data-aware low-rank compression for large nlp models.
  * Chen et al., (2021) Chen, T., Cheng, Y., Gan, Z., Yuan, L., Zhang, L., and Wang, Z. (2021).  Chasing sparsity in vision transformers: An end-to-end exploration.  ArXiv, abs/2106.04533.
  * Devlin et al., (2019) Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K. (2019).  Bert: Pre-training of deep bidirectional transformers for language understanding.  ArXiv, abs/1810.04805.
  * Dosovitskiy et al., (2021) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., and Houlsby, N. (2021).  An image is worth 16x16 words: Transformers for image recognition at scale.
  * Gou et al., (2021) Gou, J., Yu, B., Maybank, S. J., and Tao, D. (2021).  Knowledge distillation: A survey.  ArXiv, abs/2006.05525.
  * Gupta et al., (2015) Gupta, S., Agrawal, A., Gopalakrishnan, K., and Narayanan, P. (2015).  Deep learning with limited numerical precision.  In ICML.
  * Han et al., (2016) Han, S., Mao, H., and Dally, W. J. (2016).  Deep compression: Compressing deep neural network with pruning, trained quantization and huffman coding.  arXiv: Computer Vision and Pattern Recognition.
  * Han et al., (2015) Han, S., Pool, J., Tran, J., and Dally, W. J. (2015).  Learning both weights and connections for efficient neural network.  ArXiv, abs/1506.02626.
  * Hassibi and Stork, (1992) Hassibi, B. and Stork, D. G. (1992).  Second order derivatives for network pruning: Optimal brain surgeon.  In NIPS.
  * Henighan et al., (2020) Henighan, T. J., Kaplan, J., Katz, M., Chen, M., Hesse, C., Jackson, J., Jun, H., Brown, T. B., Dhariwal, P., Gray, S., Hallacy, C., Mann, B., Radford, A., Ramesh, A., Ryder, N., Ziegler, D. M., Schulman, J., Amodei, D., and McCandlish, S. (2020).  Scaling laws for autoregressive generative modeling.  ArXiv, abs/2010.14701.
  * Hinton et al., (2015) Hinton, G. E., Vinyals, O., and Dean, J. (2015).  Distilling the knowledge in a neural network.  ArXiv, abs/1503.02531.
  * Hou and Kung, (2022) Hou, Z. and Kung, S. Y. (2022).  Multi-dimensional model compression of vision transformer.  ArXiv, abs/2201.00043.
  * Kaplan et al., (2020) Kaplan, J., McCandlish, S., Henighan, T. J., Brown, T. B., Chess, B., Child, R., Gray, S., Radford, A., Wu, J., and Amodei, D. (2020).  Scaling laws for neural language models.  ArXiv, abs/2001.08361.
  * Krishnamoorthi, (2018) Krishnamoorthi, R. (2018).  Quantizing deep convolutional networks for efficient inference: A whitepaper.  ArXiv, abs/1806.08342.
  * LeCun et al., (1989) LeCun, Y., Denker, J. S., and Solla, S. A. (1989).  Optimal brain damage.  In NIPS.
  * Li et al., (2017) Li, H., Kadav, A., Durdanovic, I., Samet, H., and Graf, H. P. (2017).  Pruning filters for efficient convnets.  ArXiv, abs/1608.08710.
  * Mao et al., (2020) Mao, Y., Wang, Y., Wu, C., Zhang, C., Wang, Y., Yang, Y., Zhang, Q., Tong, Y., and Bai, J. (2020).  Ladabert: Lightweight adaptation of bert through hybrid model compression.  In COLING.
  * Micikevicius et al., (2018) Micikevicius, P., Narang, S., Alben, J., Diamos, G. F., Elsen, E., García, D., Ginsburg, B., Houston, M., Kuchaiev, O., Venkatesh, G., and Wu, H. (2018).  Mixed precision training.  ArXiv, abs/1710.03740.
  * Molchanov et al., (2019) Molchanov, P., Mallya, A., Tyree, S., Frosio, I., and Kautz, J. (2019).  Importance estimation for neural network pruning.  2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 11256–11264.
  * Nagel et al., (2021) Nagel, M., Fournarakis, M., Amjad, R. A., Bondarenko, Y., van Baalen, M., and Blankevoort, T. (2021).  A white paper on neural network quantization.  ArXiv, abs/2106.08295.
  * Radford et al., (2021) Radford, A., Kim, J. W., Hallacy, C., Ramesh, A., Goh, G., Agarwal, S., Sastry, G., Askell, A., Mishkin, P., Clark, J., Krueger, G., and Sutskever, I. (2021).  Learning transferable visual models from natural language supervision.  In ICML.
  * Radford and Narasimhan, (2018) Radford, A. and Narasimhan, K. (2018).  Improving language understanding by generative pre-training.
  * Ramesh et al., (2021) Ramesh, A., Pavlov, M., Goh, G., Gray, S., Voss, C., Radford, A., Chen, M., and Sutskever, I. (2021).  Zero-shot text-to-image generation.  ArXiv, abs/2102.12092.
  * Sainath et al., (2013) Sainath, T. N., Kingsbury, B., Sindhwani, V., Arisoy, E., and Ramabhadran, B. (2013).  Low-rank matrix factorization for deep neural network training with high-dimensional output targets.  2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 6655–6659.
  * Sanh et al., (2019) Sanh, V., Debut, L., Chaumond, J., and Wolf, T. (2019).  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.  ArXiv, abs/1910.01108.
  * Vaswani et al., (2017) Vaswani, A., Shazeer, N. M., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. (2017).  Attention is all you need.  ArXiv, abs/1706.03762.
  * Wang et al., (2018) Wang, N., Choi, J., Brand, D., Chen, C.-Y., and Gopalakrishnan, K. (2018).  Training deep neural networks with 8-bit floating point numbers.  In NeurIPS.
  * Wang et al., (2020) Wang, S., Li, B. Z., Khabsa, M., Fang, H., and Ma, H. (2020).  Linformer: Self-attention with linear complexity.  ArXiv, abs/2006.04768.
  * Wang et al., (2021) Wang, Y., Huang, R., Song, S., Huang, Z., and Huang, G. (2021).  Not all images are worth 16x16 words: Dynamic vision transformers with adaptive sequence length.  ArXiv, abs/2105.15075.
  * Wen et al., (2016) Wen, W., Wu, C., Wang, Y., Chen, Y., and Li, H. H. (2016).  Learning structured sparsity in deep neural networks.  In NIPS.
  * Yang et al., (2021) Yang, H., Yin, H., Molchanov, P., Li, H., and Kautz, J. (2021).  Nvit: Vision transformer compression and parameter redistribution.  ArXiv, abs/2110.04869.
  * Yu and Wu, (2021) Yu, H. and Wu, J. (2021).  A unified pruning framework for vision transformers.  ArXiv, abs/2111.15127.
  * Zhu et al., (2021) Zhu, M., Tang, Y., and Han, K. (2021).  Vision transformer pruning.


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

Chenghan Li  Junjie Yang  Xing Zhang  Garnet Kin-Lic Chan  Division of
Chemistry and Chemical Engineering, California Institute of Technology,
Pasadena CA 91101

###### Abstract

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 methodsSun and Chan (2016); Lee  _et al._ (2019a); Jones  _et al._
(2020); Ma  _et al._ (2021) 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)Knizia and Chan (2012, 2013); Cui, Zhu, and Chan
(2019); Pham, Hermes, and Gagliardi (2019); Hermes and Gagliardi (2019);
Tsuchimochi, Welborn, and Van Voorhis (2015); Ye  _et al._ (2018) 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 methodsLee  _et al._ (2019b); Plekhanov, Bonini, and Weber (2021),
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 nanoreactorsWang  _et al._ (2014),
combustion processesZeng  _et al._ (2020); Ashraf and Van Duin (2017), and
heterogeneous catalysisSchlögl (2015); Senftle  _et al._ (2016). 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.

###  II.1 The DMET energy expression

The theory of DMET is documented elsewhere Wouters  _et al._ (2016); Cui, Zhu,
and Chan (2019); 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), $\gamma_{\mu\nu}$. Here $\mu,\nu,\lambda,\sigma$ 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 $p,q,r,s$. As an example, the mean-field
1-RDM elements rotated into the LO basis are labelled $\gamma_{pq}$.

An impurity $I$ is defined as a subset of the full LOs. We partition the full
LO basis into multiple non-overlapping impurities. For impurity $I$, the set
of LOs not in $I$ defines its environment, which in DMET is later replaced by
a set of bath and core orbitals. The bath of $I$ is obtained from the singular
value decomposition (SVD) of the environment-impurity mean-field density
matrix

$\displaystyle\gamma_{p\centernot\in I,q\in
I}=U_{p\tilde{p}}\sigma_{\tilde{p}}V_{\tilde{p}q}$ (1)

where $U_{p\tilde{p}}$ transforms the environment LOs into the bath orbitals.
Similarly, the core orbitals for $I$ are defined from the environment-
environment block of $\gamma_{pq}$, as the eigenvectors with double occupancy.
Together, the set of impurity, bath, and core orbitals define the embedded
impurity problem for $I$; the impurity and bath orbitals are usually referred
to as the embedding orbitals (EOs) which we label with tilde (“$\sim$”) labels
$\tilde{p},\tilde{q},\tilde{r},\tilde{s},\tilde{t}$; the corresponding core
orbitals are labelled $\tilde{c},\tilde{d}$. The embedding and core orbitals
are, of course, different for each $I$, but where possible we suppress the
additional $I$ label to reduce notation. The transformation matrix from the
LOs to the EOs of fragment $I$ can be summarized as

$\displaystyle
C_{p\tilde{p}}=\begin{cases}\delta_{p\tilde{p}},~{}~{}~{}p\in\text{imp},\tilde{p}\in\text{imp}\\\
U_{p\tilde{p}},~{}~{}p\in\text{env},\tilde{p}\in\text{bath}\\\
0,~{}~{}~{}~{}\text{otherwise}\end{cases}$ (2)

The embedding Hamiltonian $H^{I}$ for impurity $I$ is defined by projecting
the original Hamiltonian into the space of determinants spanned by
$\\{\tilde{p}\\}\otimes|\mathrm{core}\rangle$ (where $|\mathrm{core}\rangle$
is the doubly occupied core Slater determinant constructed from the core
orbitals associated with $I$). Since there are different EOs for each
impurity, there is a different $H^{I}$ for each impurity. In second
quantization, we write

$\displaystyle
H^{I}=h_{\tilde{p}\tilde{q}}a_{\tilde{p}}^{\dagger}a_{\tilde{q}}-\mu
a_{\tilde{p}}^{\dagger}a_{\tilde{p}}~{}(\tilde{p}\in\text{imp})+{\frac{1}{2}}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}a_{\tilde{p}}^{\dagger}a_{\tilde{r}}^{\dagger}a_{\tilde{s}}a_{\tilde{q}}$
(3)

where we have included a chemical potential (see below). The one-electron part
contains the Coulomb and exchange contributions from the core orbitals,

$\displaystyle h_{\tilde{p}\tilde{q}}$
$\displaystyle=t_{\tilde{p}\tilde{q}}+V_{\tilde{p}\tilde{q}\tilde{c}\tilde{d}}\gamma_{\tilde{c}\tilde{d}}$
$\displaystyle=F_{\tilde{p}\tilde{q}}-V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}{\gamma}_{\tilde{r}\tilde{s}}$
(4)

where $V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}$ is the anti-symmetrized
electron repulsion integral (ERI)
($(\tilde{p}\tilde{q}|\tilde{s}\tilde{r})-\frac{1}{2}(\tilde{p}\tilde{r}|\tilde{s}\tilde{q})$
in chemists’ notation) and $t_{\tilde{p}\tilde{q}}$, $F_{\tilde{p}\tilde{q}}$
correspond to the bare one-electron and Fock operator rotated into the EO
basis respectively. The chemical potential $\mu$ 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 $H^{I}$, 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 $\Gamma_{\tilde{p}\tilde{q}}^{I}$ and
two-body reduced density matrices (2-RDMs)
$\Gamma^{I}_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}$. 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 $\Gamma^{I}_{\tilde{p}\tilde{q}}$:

$\displaystyle\Gamma_{pq}^{\prime}=C_{p\tilde{p}}^{I}\Gamma_{\tilde{p}\tilde{q}}^{I}w_{\tilde{p}\tilde{q}}^{I}C_{q\tilde{q}}^{I}$
(5)

The weight factor $w_{\tilde{p}\tilde{q}}$ is $1$ for
$\tilde{p},\tilde{q}\in\mathrm{imp}$, $0$ for
$\tilde{p},\tilde{q}\in\mathrm{bath}$, and $\frac{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 $\mu$ such that
$\mathrm{tr}\ \Gamma^{\prime}=N_{\text{elec}}$, where $N_{\text{elec}}$ is the
desired global number of electrons. We then update the one-electron part of
the embedding Hamiltonian,

$\displaystyle h_{\tilde{p}\tilde{q}}^{\prime}$
$\displaystyle=C_{p\tilde{p}}(t_{pq}+{\frac{1}{2}}V_{pqrs}\Gamma_{rs}^{\prime})C_{q\tilde{q}}-{\frac{1}{2}}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}\Gamma_{\tilde{r}\tilde{s}}$
(6)

and compute the total DMET energy as an average over all fragments,

$\displaystyle{E^{\mathrm{DMET}}}=h_{\tilde{p}\tilde{q}}^{\prime
I}w_{\tilde{p}\tilde{q}}^{I}\Gamma_{\tilde{q}\tilde{p}}^{I}+{\frac{1}{2}}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}^{I}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}$
(7)

where the four-index weight $w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}$ is
defined as $\frac{1}{4}$ times the number of impurity labels in the set
$\\{\tilde{p},\tilde{q},\tilde{r},\tilde{s}\\}$.

###  II.2 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$\rightarrow$LO and LO$\rightarrow$EO transformations
are both atomic-position-dependent, and thus their response also contributes
to the energy gradient. In particular, the LO$\rightarrow$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$\rightarrow$LO transformation can be a function of AO overlap (such as for
Löwdin orthogonalizationLöwdin (1950)) and/or the mean-field solution or
additional equations (such as for intrinsic atomic orbitalsKnizia (2013); West
_et al._ (2013)). In both cases, the AO$\rightarrow$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 $\boldsymbol{c}=\\{\boldsymbol{c}_{I}\\}$ (these could denote,
e.g. the configuration interaction (CI) vectors for a full CI (FCI) solver, or
the $\boldsymbol{t}$ amplitudes plus embedding MO coefficients for a coupled
cluster (CC) solver), the global mean-field solution (MO coefficients)
$\boldsymbol{C}$, and the AO integrals $\boldsymbol{A}$ (the AO 1-electron and
2-electron integrals $\boldsymbol{\mathcal{I}}$ and AO overlap integrals
$\boldsymbol{S}$). The gradient (with respect to nuclear positions or any
perturbation) can be expressed in terms of partial derivatives with respect to
these variables:

$\displaystyle\frac{{\mathrm{d}}{E^{\mathrm{DMET}}}}{{\mathrm{d}}{\boldsymbol{R}}}={\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{c}}}{\frac{{\mathrm{d}}\boldsymbol{c}}{{\mathrm{d}}{\boldsymbol{R}}}}+{\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{C}}}{\frac{{\mathrm{d}}\boldsymbol{C}}{{\mathrm{d}}{\boldsymbol{R}}}}+{\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}$
(8)

Note that ${E^{\mathrm{DMET}}}$ does not depend on the chemical potential
$\mu$ explicitly, but $\mu$ affects the energy indirectly through its effect
on $\boldsymbol{c}$. The
${\mathrm{d}}\boldsymbol{c}/{\mathrm{d}}{\boldsymbol{R}}$ and
${\mathrm{d}}\boldsymbol{C}/{\mathrm{d}}{\boldsymbol{R}}$ contributions are
related to ${\mathrm{d}}\boldsymbol{A}/{\mathrm{d}}{\boldsymbol{R}}$ since
$\boldsymbol{c}$ and $\boldsymbol{C}$ are implicit functions of
$\boldsymbol{A}$ 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
${\mathrm{d}}\boldsymbol{c}/{\mathrm{d}}{\boldsymbol{R}}$ and
${\mathrm{d}}\boldsymbol{C}/{\mathrm{d}}{\boldsymbol{R}}$ for the Cartesian
coordinates of each atom, i.e. $3N_{\text{atom}}$ times. This complexity can
be greatly reduced, however, by employing the so-called $Z$-vector technique,
also known as the Lagrangian methodHandy and Schaefer III (1984); Yamaguchi
_et al._ (1994).

The DMET ground-state solutions involve three sets of working equations: the
solution conditions for the impurity ground-states

$\displaystyle\boldsymbol{X}(\mu,\boldsymbol{c},\boldsymbol{C},\boldsymbol{A})=\\{\boldsymbol{X}^{I}(\mu,\boldsymbol{c}^{I},\boldsymbol{C},\boldsymbol{A})\\}=\boldsymbol{0}$
(9)

where each $\boldsymbol{X}^{I}$ is the corresponding impurity solution
condition (such as the CI, CC ground-state equations, etc.); the mean-field
solution condition;

$\displaystyle{\boldsymbol{M}}(\boldsymbol{C},\boldsymbol{A})=\boldsymbol{0}$
(10)

and the chemical potential solution condition;

$\displaystyle N(\boldsymbol{c})=\mathrm{tr}\ \Gamma^{\prime}=N_{\text{elec}}$
(11)

We define a Lagrangian associated with the above constraints,

$\displaystyle\mathcal{L}={E^{\mathrm{DMET}}}-z(N(\boldsymbol{c})-N_{\mathrm{elec}})-\boldsymbol{\lambda}^{\dagger}\boldsymbol{X}(\mu,\boldsymbol{c},\boldsymbol{C},\boldsymbol{A})-\boldsymbol{\alpha}^{\dagger}{\boldsymbol{M}}(\boldsymbol{C},\boldsymbol{A})$
(12)

and consider the stationary points w.r.t. all the variables. Vanishing
derivatives of $\mathcal{L}$ w.r.t $z$, $\boldsymbol{\lambda}$ and
$\boldsymbol{\alpha}$ give the working equations of DMET, while vanishing
derivatives of $\mathcal{L}$ w.r.t $\mu$, $\boldsymbol{c}$ and
$\boldsymbol{C}$ give equations for the minimizers of the multipliers,
$z_{0}$, $\boldsymbol{\lambda}_{0}$, $\boldsymbol{\alpha}_{0}$ respectively,

$\displaystyle(\boldsymbol{\lambda}_{0}~{}z_{0})\begin{pmatrix}{\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{c}}}&{\frac{\partial\boldsymbol{X}}{\partial\mu}}\\\
{\frac{\partial
N}{\partial\boldsymbol{c}}}&0\end{pmatrix}=\Big{(}{\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{c}}}~{}~{}0\Big{)}$
(13)

and

$\displaystyle\boldsymbol{\alpha}_{0}^{\dagger}=\Big{(}{\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{\boldsymbol{C}}}}-\boldsymbol{\lambda}_{0}^{\dagger}{\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{\boldsymbol{C}}}}\Big{)}\Big{(}{\frac{\partial{\boldsymbol{M}}}{\partial\boldsymbol{\boldsymbol{C}}}}\Big{)}^{-1}$
(14)

If no chemical potential is used, Eq. 13 reduces to

$\displaystyle\boldsymbol{\lambda}_{0}^{\dagger}={\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{c}}}\Big{(}{\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{c}}}\Big{)}^{-1}$
(15)

which can be solved for each impurity independently due to the block-diagonal
nature of $\partial{\boldsymbol{X}}/\partial\boldsymbol{c}$, i.e.,
$\partial{\boldsymbol{X}}/\partial{\boldsymbol{c}}=\text{diag}(\\{\partial{\boldsymbol{X}^{I}}/\partial{\boldsymbol{c}^{I}}\\})$.

At the stationary points of $\mathcal{L}$, $\mathcal{L}={E^{\mathrm{DMET}}}$
and
${\mathrm{d}}{E^{\mathrm{DMET}}}/{\mathrm{d}}{\boldsymbol{R}}={\mathrm{d}}\mathcal{L}/{\mathrm{d}}{\boldsymbol{R}}$.
However, the nuclear gradient of $\mathcal{L}$ no longer contains the
contributions from the partial derivatives w.r.t $\boldsymbol{c}$ or
$\boldsymbol{C}$ as in Eq. 8, as their implicit dependence on
${\boldsymbol{R}}$ is accounted for by the multipliers. Thus,

$\displaystyle{\frac{{\mathrm{d}}{E^{\mathrm{DMET}}}}{{\mathrm{d}}{\boldsymbol{R}}}}={\frac{\partial\mathcal{L}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}$
$\displaystyle={\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}-\boldsymbol{\lambda}_{0}^{\dagger}{\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}-\boldsymbol{\alpha}_{0}^{\dagger}{\frac{\partial{\boldsymbol{M}}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}$
(16)

We note the repeated occurrence of
${E^{\mathrm{DMET}}}-\boldsymbol{\lambda}_{0}^{\dagger}\boldsymbol{X}$ in Eqs.
14 and 16, thus it is useful to define
$\mathcal{L}_{\mathrm{emb}}={E^{\mathrm{DMET}}}-\boldsymbol{\lambda}_{0}^{\dagger}\boldsymbol{X}$,
then its derivative w.r.t.
$\boldsymbol{Y}\in\\{\boldsymbol{C},\boldsymbol{A}$} is

$\displaystyle{\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial\boldsymbol{Y}}}={\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{Y}}}-\boldsymbol{\lambda}_{0}^{\dagger}{\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{Y}}}=\Gamma_{\tilde{q}\tilde{p}}^{I}w_{\tilde{p}\tilde{q}}^{I}{\frac{\partial
h_{\tilde{p}\tilde{q}}^{\prime
I}}{\partial\boldsymbol{Y}}}+{\frac{1}{2}}\Big{(}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}+D_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}\Big{)}{\frac{\partial{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}^{I}}{\partial\boldsymbol{Y}}}+D_{\tilde{q}\tilde{p}}^{I}{\frac{\partial
h_{\tilde{p}\tilde{q}}^{I}}{\partial\boldsymbol{Y}}}$ (17)

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
$D_{\tilde{p}\tilde{q}}\overset{\mathrm{def}}{=\joinrel=}-\boldsymbol{\lambda}_{0}^{\dagger}\partial\boldsymbol{X}/\partial
h_{\tilde{p}\tilde{q}}$ and
$D_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}\overset{\mathrm{def}}{=\joinrel=}-2\boldsymbol{\lambda}_{0}^{\dagger}\partial\boldsymbol{X}/\partial{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}$,
which account for all solver-specific details. In the Appendix, we give the
corresponding expressions for the Hartree-Fock (HF) and FCI solvers as
examples. Choosing $\boldsymbol{Y}$ to be the AO integrals
$\boldsymbol{\mathcal{I}}$ in Eq. 17 gives the “skeleton” derivative of the
DMET gradient, while choosing $\boldsymbol{Y}$ to be the AO overlap leads to
more complications, due to the orbital dependence on $\boldsymbol{S}$ in
localization, as discussed in Appendix VI.4. These two components sum to the
first two terms of Eq. 16. The remaining task is to solve for
$\boldsymbol{\alpha}$ from
$\partial\mathcal{L}_{\mathrm{emb}}/\partial\boldsymbol{C}$ using the coupled-
perturbed Hartree-Fock (CP-HF) equations Peng (1941); Yamaguchi  _et al._
(1994), and to consider the contributions from
$\partial\boldsymbol{M}/\partial\boldsymbol{A}$. Further details can be found
in Appendix VI.3.

##  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 $\mu$)).

###  III.1 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.

####  III.1.1 Water trimer

We take a water trimer as our first example (configuration shown in Fig. 1A).
We used the 6-31G** basis setHariharan and Pople (1973) 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 basisHehre, Stewart,
and Pople (1969) 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 ($\delta x$).

####  III.1.2 Hydrogen ring

We next consider the symmetric dissociation of a hydrogen ring Knizia and Chan
(2013). We treat the $\mathrm{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).

Figure 1: Numerical validation of the DMET gradient. (A) Analytic DMET nuclear
gradient compared with finite difference results using different atomic
position displacements ($\delta x$) in a water trimer (geometry shown). The
gradient with and without chemical potential ($\mu$) fitting is shown in red
and blue respectively. The solid lines indicate the difference in the norm of
the gradient and the dashed lines indicate the MAE. The MAE tends to saturate
with small $\delta x$ due to the accuracy limit of solving the response
equations. (B) PES of hydrogen ring symmetric dissociation. Solid curves
represent direct calculations of the energies, while crosses represent the PES
obtained by numerically integrating the DMET analytic gradient.

Since the electron number in DMET (without $\mu$) changes only slightly along
the symmetric stretching process (MAE = $\sim$0.02 with a standard deviation
of $\sim$0.02) fitting $\mu$ here only slightly affects the energy.

###  III.2 Diels-Alder reaction

We next use the prototypical Diels-Alder reactionFringuelli and Taticchi
(2002) 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 protocolYang, Doubleday, and
Houk (2015). We used the 6-31G basis setHehre, Stewart, and Pople (1969) 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 D3Grimme  _et al._ (2010) 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 $\mu$) and the gradient was further used in a nudged
elastic band (NEB)Berne, Ciccotti, and Coker (1998); Henkelman and Jónsson
(2000); Henkelman, Uberuaga, and Jónsson (2000) 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
ASELarsen  _et al._ (2017).

Figure 2: Optimized structures and energies of a Diels-Alder reaction. R, TS,
and P represent reactants, transition state and products, respectively. Energy
levels for each method are aligned at R. Note that the black and green lines
overlap with each other at TS. The gray and black, and red and purple overlap
at P.

Fig. 2 shows the DMET-D3 optimized structures and associated energies (see
also Tab. 1). In this case, DMET without fitting $\mu$ 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)Kowalski and Piecuch (2000), 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-2XZhao and Truhlar (2006))
(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).

| $E(\text{TS})-E(\text{R})$ | $E(\text{P})-E(\text{R})$
---|---|---
R-CCSD(T) | 20.97 | -16.41
CCSD | 23.78 | -16.26
HF | 33.97 | -13.91
DMET-D3 | 21.60 | -24.05
DMET | 23.91 | -23.12
M06-2X | 13.93 | -23.98
Table 1: Barrier and reaction energies of a simple Diels-Alder reaction in
kcal/mol. R, TS, P denote reactants, transition state, and products.

###  III.3 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 high-level quantum methods. In PT, the multi-
site character arises due to the Grotthuss proton hopping mechanismDe
Grotthuss (2006) 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.

Figure 3: PES of proton transport in a water tetramer along one MD trajectory
segment. DFT embedding vertical means the system is partitioned into H5O2+ on
the right and a water dimer on the left, while horizontal means it is
partitioned into an upper H5O2+ plus a lower water dimer. The DFT embedding is
performed with the projection-based wavefunction-in-DFT method using a CCSD(T)
solver embedded in PBE DFTPerdew, Burke, and Ernzerhof (1996) implemented in
MolproWerner  _et al._ (2020). The DMET, HF, and DFT embedding/vertical
trajectories are aligned to CCSD(T) using the averaged energy in the 0.8-0.9
ps segment. The DFT embedding/horizontal trajectory is aligned using the
0.95-1.05 ps segment energy average. Lower figure shows the impact of chemical
potential fitting on the DMET trajectory energy.

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 =
$\sim$0.06 with standard deviation $\sim$0.002) thus chemical potential
fitting is seen to change the energy only slightly. We used parallel-biased
metadynamicsPfaendtner and Bonomi (2015) 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
positionLi and Swanson (2020); Li and Voth (2021). 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-
moleculesVoth (2006); Berne, Ciccotti, and Coker (1998); Watkins  _et al._
(2019); Lee  _et al._ (2018); Maag  _et al._ (2021); Liu  _et al._ (2021); Li
_et al._ (2022) 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
H5O2+ complex (i.e., whole system is partitioned vertically) to the upper
H5O2+ 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 Lee  _et
al._ (2019a) 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.

Figure 4: PES of proton transport in a water tetramer along one MD trajectory
segment. Figure labels are the same as those in 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.

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
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. $C_{\mu\tilde{p}}$ is interpreted as
the AO$\rightarrow$EO transformation such that $|\tilde{p}\rangle=|\mu\rangle
C_{\mu\tilde{p}}$. We will also assume the orbital coefficients are real,
although the equations can be straightforwardly generalized to complex-valued
orbitals.

###  VI.1 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 $\boldsymbol{\lambda}$ in Eq. 13. Consider $\boldsymbol{c}^{J}$, the
impurity solution of fragment $J$. Then differentiating Eq. 7 gives

$\displaystyle{\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{c}^{J}}}$
$\displaystyle=h_{\tilde{p}\tilde{q}}^{\prime
J}w_{\tilde{p}\tilde{q}}^{J}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}}^{J}}{\partial\boldsymbol{c}^{J}}}+{\frac{1}{2}}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}^{J}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}}{\partial\boldsymbol{c}^{J}}}+\Gamma_{\tilde{p}\tilde{q}}^{I}w_{\tilde{p}\tilde{q}}^{I}{\frac{\partial
h_{\tilde{p}\tilde{q}}^{\prime I}}{\partial\boldsymbol{c}^{J}}}$ (no summation
over $J$) (18)

where we have a sum over $I$ in the last term because $h^{\prime I}$ depends
on the solutions for all fragments (see Eq. 6), $\boldsymbol{c}$, via its
dependence on the updated global 1-RDM $\Gamma^{\prime}$. Explicitly,

$\displaystyle\Gamma_{\tilde{p}\tilde{q}}^{I}w_{\tilde{p}\tilde{q}}^{I}{\frac{\partial
h_{\tilde{p}\tilde{q}}^{\prime I}}{\partial\boldsymbol{c}^{J}}}$
$\displaystyle=\Big{(}{\frac{1}{2}}C_{p\tilde{p}}^{I}C_{q\tilde{q}}^{I}V_{pqrs}{\frac{\partial\Gamma_{rs}^{\prime}}{\partial\boldsymbol{c}^{J}}}-{\frac{1}{2}}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{r}\tilde{s}}^{I}}{\partial\boldsymbol{c}^{J}}}\delta_{IJ}\Big{)}\Gamma_{\tilde{p}\tilde{q}}^{I}w_{\tilde{p}\tilde{q}}^{I}$
$\displaystyle={\frac{1}{2}}\Gamma_{\tilde{p}\tilde{q}}^{I}w_{\tilde{p}\tilde{q}}^{I}C_{p\tilde{p}}^{I}C_{q\tilde{q}}^{I}V_{pqrs}w_{\tilde{r}\tilde{s}}^{J}C_{r\tilde{r}}^{J}C_{s\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{r}\tilde{s}}^{J}}{\partial\boldsymbol{c}^{J}}}-{\frac{1}{2}}\Gamma_{\tilde{p}\tilde{q}}^{J}w_{\tilde{p}\tilde{q}}^{J}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{r}\tilde{s}}^{J}}{\partial\boldsymbol{c}^{J}}}$
$\displaystyle={\frac{1}{2}}v_{pq}^{\prime}w_{\tilde{p}\tilde{q}}^{J}C_{p\tilde{p}}^{J}C_{q\tilde{q}}^{J}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}}^{J}}{\partial\boldsymbol{c}^{J}}}-{\frac{1}{2}}\Gamma_{\tilde{r}\tilde{s}}^{J}w_{\tilde{r}\tilde{s}}^{J}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}}^{J}}{\partial\boldsymbol{c}^{J}}}$
(19) (no summation over $J$)

where we have summed over $I$ first in
$\Gamma^{I}_{\tilde{p}\tilde{q}}w_{\tilde{p}\tilde{q}}^{I}C_{p\tilde{p}}^{I}C_{q\tilde{q}}^{I}$,
used Eq. 5, and performed some dummy variable changes, and defined the
effective potential $v_{rs}^{\prime}=\Gamma_{pq}^{\prime}V_{pqrs}$. Inserting
back into Eq. 18,

$\displaystyle{\frac{\partial{E^{\mathrm{DMET}}}}{\partial\boldsymbol{c}^{J}}}$
$\displaystyle=\Big{[}h_{\tilde{p}\tilde{q}}^{\prime
J}w_{\tilde{p}\tilde{q}}^{J}-{\frac{1}{2}}\Gamma_{\tilde{r}\tilde{s}}^{J}w_{\tilde{r}\tilde{s}}^{J}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}+{\frac{1}{2}}v_{pq}^{\prime}w_{\tilde{p}\tilde{q}}^{J}C_{p\tilde{p}}^{J}C_{q\tilde{q}}^{J}\Big{]}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}}^{J}}{\partial\boldsymbol{c}^{J}}}+{\frac{1}{2}}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}^{J}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}}{\partial\boldsymbol{c}^{J}}}$
$\displaystyle\overset{\mathrm{def}}{=\joinrel=}h^{\prime\prime
J}_{\tilde{p}\tilde{q}}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}}^{J}}{\partial\boldsymbol{c}^{J}}}+{\frac{1}{2}}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}^{J}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{J}}{\partial\boldsymbol{c}^{J}}}$
(20) (no summation over $J$)

where the terms in the square brackets define $h_{ij}^{\prime\prime J}$.

In the case of no chemical potential fitting, $\boldsymbol{\lambda}$ is
directly obtained after a contraction with
$(\partial{\boldsymbol{M}}/\partial{\boldsymbol{c}})^{-1}$ (where the
contraction is performed by solving the response equation of the underlying
solver method, such as coupled perturbed CI (CP-CI)Yamaguchi  _et al._ (1994)
for a CI solver, or the $\Lambda$-equations for CC solversShavitt and Bartlett
(2009)). If chemical potential fitting is performed, the solver responses of
different impurities are connected via a common chemical potential. Thus their
$\boldsymbol{\lambda}$’s are solved for at the same time as $z$, the
multiplier associated with the electron number constraint (Eq. 13).

####  VI.1.1 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. VI.1)
is (we consider
$(\partial{E^{\mathrm{DMET}}}/\partial\boldsymbol{C})\boldsymbol{C}$ instead
of $\partial{E^{\mathrm{DMET}}}/\partial\boldsymbol{C}$ to connect to the
well-known CP-HF equation)

$\displaystyle{\frac{\partial{E^{\mathrm{DMET}}}}{\partial
C_{\tilde{p}\tilde{i}}}}C_{\tilde{p}\tilde{a}}={\frac{\partial{E^{\mathrm{DMET}}}}{\partial\Gamma_{\tilde{p}\tilde{q}}}}{\frac{\partial\Gamma_{\tilde{p}\tilde{q}}}{\partial
C_{\tilde{r}\tilde{i}}}}C_{\tilde{r}\tilde{a}}=(f^{\prime}_{\tilde{a}\tilde{i}}+f^{\prime}_{\tilde{i}\tilde{a}})$
(21)

In the above, we use $\tilde{i}$ and $\tilde{a}$ to label the Hartree-Fock
impurity orbitals, i.e. the MOs in the embedding space ($\tilde{i}$,
$\tilde{j}$ will be used for occupied and $\tilde{a}$, $\tilde{b}$ for virtual
orbitals) while $\tilde{p}$, $\tilde{q}$, $\tilde{r}$ and $\tilde{s}$,
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 $C_{\tilde{p}\tilde{i}}$ or
$C_{\tilde{p}\tilde{a}}$, and the embedding “Fock” matrix is defined as

$\displaystyle
f^{\prime}_{\tilde{p}\tilde{q}}={\frac{\partial{E^{\mathrm{DMET}}}}{\partial\Gamma_{\tilde{p}\tilde{q}}}}=h^{\prime\prime}_{\tilde{p}\tilde{q}}+V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}\Gamma_{\tilde{r}\tilde{s}}$
(22)

The CP-HF equation is usually written as

$\displaystyle\mathcal{A}_{{\tilde{b}}{\tilde{j}},{\tilde{a}}{\tilde{i}}}U_{{\tilde{a}}{\tilde{i}}}^{\boldsymbol{R}}=B^{\boldsymbol{R}}_{{\tilde{b}}{\tilde{j}}}$
(23)

where $\mathcal{A}$ and $B^{\boldsymbol{R}}$ in our language are
$\mathcal{A}=-{\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{C}}}\boldsymbol{C}$
and
$B^{\boldsymbol{R}}={\frac{\partial\boldsymbol{X}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}$
with $X_{\tilde{a}\tilde{i}}=f_{\tilde{a}\tilde{i}}=0$ the Hartree-Fock
equation of the embedding problem. The solution
$\boldsymbol{U}^{\boldsymbol{R}}$ is related to MO response via
$\boldsymbol{C}^{-1}({\mathrm{d}}{\boldsymbol{C}}/{\mathrm{d}}{{\boldsymbol{R}}})$.
Comparing to Eq. 15, if no $\mu$-fit is performed, $\boldsymbol{\lambda}$ for
the HF solver is obtained from the CP-HF equation by replacing
$B^{\boldsymbol{R}}$ with $f^{\prime}$, i.e.,

$\displaystyle-\lambda_{\tilde{a}\tilde{i}}=2f^{\prime}_{\tilde{b}\tilde{j}}\mathcal{A}^{-1}_{\tilde{b}\tilde{j},\tilde{a}\tilde{i}}$
(24)

The auxiliary RDMs for the HF solver (defined after Eq. 17) are then obtained
by contracting with ${\frac{\partial\boldsymbol{X}}{\partial
h_{\tilde{p}\tilde{q}}}}$ and
${\frac{\partial\boldsymbol{X}}{\partial{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}}}$

$\displaystyle D_{\tilde{p}\tilde{q}}$
$\displaystyle=\lambda_{\tilde{a}\tilde{i}}C_{\tilde{p}\tilde{a}}C_{\tilde{q}\tilde{i}}$
(25) $\displaystyle D_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}$
$\displaystyle=D_{\tilde{p}\tilde{q}}\Gamma_{\tilde{r}\tilde{s}}-{\frac{1}{2}}D_{\tilde{p}\tilde{s}}\Gamma_{\tilde{q}\tilde{r}}+\Gamma_{\tilde{p}\tilde{q}}D_{\tilde{r}\tilde{s}}-{\frac{1}{2}}\Gamma_{\tilde{p}\tilde{s}}D_{\tilde{q}\tilde{r}}$
(26)

If a $\mu$-fitting is performed, the above auxiliary RDM expressions still
hold but $\boldsymbol{\lambda}$ should be obtained from Eq. 13 using all the
impurity response equations simultaneously. The other information needed in
Eq. 13 is

$\displaystyle{\frac{\partial N}{\partial
C_{\tilde{p}\tilde{i}}}}=\begin{cases}2C_{\tilde{p}\tilde{i}},~{}~{}~{}~{}\tilde{p}\in\text{imp}\\\
0,~{}~{}~{}~{}~{}~{}~{}~{}\text{otherwise}\end{cases}$ (27)

and

$\displaystyle{\frac{\partial
X_{\tilde{a}\tilde{i}}}{\partial\mu}}=-\sum_{\tilde{p}\in\text{imp}}C_{\tilde{p}\tilde{a}}C_{\tilde{p}\tilde{i}}$
(28)

####  VI.1.2 FCI solver response

The FCI vector $c_{\mathcal{P}}$ is the solution of the equations,

$\displaystyle H_{\mathcal{P}\mathcal{Q}}c_{\mathcal{Q}}$
$\displaystyle=E^{\mathrm{CI}}c_{\mathcal{P}}$
$\displaystyle\sum_{\mathcal{P}}c_{\mathcal{P}}^{2}$ $\displaystyle=1$ (29)

where we use $\mathcal{P}$ and $\mathcal{Q}$ to label an electronic
configuration (e.g. a determinant). The RDMs are constructed from
$c_{\mathcal{P}}$ via

$\displaystyle\Gamma_{\tilde{q}\tilde{p}}=c_{\mathcal{P}}c_{\mathcal{Q}}\Gamma^{\mathcal{PQ}}_{\tilde{q}\tilde{p}}$
(30)
$\displaystyle\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}=c_{\mathcal{P}}c_{\mathcal{Q}}\Gamma^{\mathcal{PQ}}_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}$
(31)

where
$\Gamma^{\mathcal{PQ}}_{\tilde{q}\tilde{p}}=\langle\Phi_{\mathcal{P}}|a^{\dagger}_{\tilde{p}}a_{\tilde{q}}|\Phi_{\mathcal{Q}}\rangle$
and
$\Gamma^{\mathcal{PQ}}_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}=\langle\Phi_{\mathcal{P}}|a^{\dagger}_{\tilde{p}}a^{\dagger}_{\tilde{r}}a_{\tilde{s}}a_{\tilde{q}}|\Phi_{\mathcal{Q}}\rangle$.
Note that since the FCI solution only depends on the space of $\mathcal{P}$,
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. VI.1 becomes

$\displaystyle{\frac{\partial{E^{\mathrm{DMET}}}}{\partial c_{\mathcal{P}}}}$
$\displaystyle=2c_{\mathcal{Q}}\Big{(}h^{\prime\prime}_{\tilde{p}\tilde{q}}\Gamma_{\tilde{p}\tilde{q}}^{\mathcal{Q}\mathcal{P}}+{\frac{1}{2}}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{\mathcal{Q}\mathcal{P}}\Big{)}\overset{\mathrm{def}}{=\joinrel=}b_{\mathcal{P}}$
(32)

which has the form of a CI Hamiltonian contracted with the CI vector
$2c_{\mathcal{Q}}$ to give a new CI vector. Note that if there are no weights
$w$ and if $h$ instead of $h^{\prime}$ is used to compute the energy, Eq. 32
simply becomes
$2c_{\mathcal{Q}}H_{\mathcal{P}\mathcal{Q}}=2E^{\mathrm{CI}}c_{\mathcal{P}}$,
which vanishes when dotted with
$\mathrm{d}c_{\mathcal{P}}/\mathrm{d}{\boldsymbol{R}}$ due to the
orthogonality between $\boldsymbol{c}$ and
$\mathrm{d}\boldsymbol{c}/\mathrm{d}{\boldsymbol{R}}$. 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 $c_{\mathcal{P}}$ and after
some linear algebra, we have

$\displaystyle{\frac{\partial X_{\mathcal{Q}}}{\partial
c_{\mathcal{P}}}}=H_{\mathcal{Q}\mathcal{P}}-E^{\mathrm{CI}}\delta_{\mathcal{P}\mathcal{Q}}+2c_{\mathcal{P}}c_{\mathcal{Q}}\overset{\mathrm{def}}{=\joinrel=}\mathcal{X}_{\mathcal{P}\mathcal{Q}}$
(33)

In the case of no $\mu$-fitting, we obtain

$\displaystyle\lambda_{\mathcal{Q}}=b_{\mathcal{P}}\mathcal{X}^{-1}_{\mathcal{P}\mathcal{Q}}$
(34)

i.e., $\boldsymbol{\lambda}$ is the solution of a CP-CI equation with
$b_{\mathcal{P}}$. Given that

$\displaystyle{\frac{\partial X_{\mathcal{P}}}{\partial
h_{\tilde{p}\tilde{q}}}}$
$\displaystyle=\Gamma_{\tilde{q}\tilde{p}}^{\mathcal{P}\mathcal{Q}}c_{\mathcal{Q}}-\Gamma_{\tilde{q}\tilde{p}}c_{\mathcal{P}}$
(35) $\displaystyle{\frac{\partial
X_{\mathcal{P}}}{\partial{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}}}$
$\displaystyle={\frac{1}{2}}(\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{\mathcal{P}\mathcal{Q}}c_{\mathcal{Q}}-\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}c_{\mathcal{P}})$
(36)

the auxiliary RDMs defined in Eq. 17 become

$\displaystyle
D_{\tilde{q}\tilde{p}}=\lambda_{\mathcal{P}}c_{\mathcal{P}}\gamma_{\tilde{q}\tilde{p}}-\lambda_{\mathcal{P}}c_{\mathcal{Q}}\gamma_{\tilde{q}\tilde{p}}^{\mathcal{P}\mathcal{Q}}$
(37) $\displaystyle
D_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}=\lambda_{\mathcal{P}}c_{\mathcal{P}}\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}-\lambda_{\mathcal{P}}c_{\mathcal{Q}}\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{\mathcal{P}\mathcal{Q}}$
(38)

To solve for $\boldsymbol{\lambda}$ with $\mu$-fitting from Eq. 13, we also
need

$\displaystyle{\frac{\partial N}{\partial
c_{\mathcal{P}}}}=\sum_{\tilde{p}\in\text{imp}}2c_{\mathcal{P}}\delta_{\tilde{p}\in\mathcal{P}}$
(39)

where $\delta_{\tilde{p}\in\mathcal{P}}=1$ if and only if $\tilde{p}$ is
occupied in configuration $\mathcal{P}$, and

$\displaystyle{\frac{\partial
X_{\mathcal{P}}}{\partial\mu}}=\sum_{\tilde{p}\in\text{imp}}\Gamma_{\tilde{p}\tilde{p}}c_{\mathcal{P}}-\Gamma_{\tilde{p}\tilde{p}}^{\mathcal{P}\mathcal{Q}}c_{\mathcal{Q}}$
(40)

###  VI.2 Orbital transformation gradient

The DMET gradient includes contributions from the response of the orbital
transformations. A position change induces both AO$\rightarrow$LO and
LO$\rightarrow$EO changes that contribute to a change in the energy. This
response is encoded in the derivative of $\mathcal{L}_{\text{emb}}$ w.r.t the
AO$\rightarrow$LO and LO$\rightarrow$EO coefficients, while both of them, via
the chain rule, are related to the derivative
$\partial{\mathcal{L}_{\mathrm{emb}}}/\partial{C_{\mu\tilde{p}}}$. Following
our conventions for basis indexing, $C_{\mu\tilde{p}}=C_{\mu p}C_{p\tilde{p}}$
is the transformation matrix from AO to EO, and $C_{\mu p}$ represents the
AO$\rightarrow$LO and $C_{p\tilde{p}}$ represents the LO$\rightarrow$EO
transformation as defined in Eq. 2.

We write the three components of
$\partial{\mathcal{L}_{\mathrm{emb}}}/\partial{C_{\mu\tilde{p}}}$ appearing in
Eq. 17 as

$\displaystyle D_{\tilde{q}\tilde{p}}{\frac{\partial
h_{\tilde{p}\tilde{q}}}{\partial C_{\mu\tilde{t}}}}$
$\displaystyle=2D_{\tilde{q}\tilde{t}}F_{\mu
q}C_{q\tilde{q}}-2D_{\tilde{q}\tilde{p}}V_{\tilde{p}\tilde{q}\tilde{t}\tilde{s}}C_{q\tilde{s}}\gamma_{\nu
q}S_{\mu\nu}$
$\displaystyle-\Big{(}D_{\tilde{q}\tilde{t}}V_{\mu\tilde{q}\tilde{r}\tilde{s}}\gamma_{\tilde{r}\tilde{s}}+D_{\tilde{q}\tilde{t}}V_{\tilde{q}\mu\tilde{r}\tilde{s}}\gamma_{\tilde{r}\tilde{s}}+D_{\tilde{q}\tilde{r}}V_{\tilde{r}\tilde{q}\mu\tilde{s}}\gamma_{\tilde{t}\tilde{s}}+D_{\tilde{q}\tilde{s}}V_{\tilde{s}\tilde{q}\tilde{r}\mu}\gamma_{\tilde{r}\tilde{t}}\Big{)}$
$\displaystyle=2\Big{(}D_{\tilde{q}\tilde{t}}F_{\mu\tilde{q}}-D_{\tilde{q}\tilde{p}}V_{\tilde{p}\tilde{q}\tilde{t}\tilde{s}}\gamma_{\nu\tilde{s}}S_{\mu\nu}-D_{\tilde{q}\tilde{t}}V_{\mu\tilde{q}\tilde{r}\tilde{s}}\gamma_{\tilde{r}\tilde{s}}-\gamma_{\tilde{q}\tilde{t}}V_{\mu\tilde{q}\tilde{r}\tilde{s}}D_{\tilde{r}\tilde{s}}\Big{)}$
(41)
$\displaystyle{\frac{1}{2}}(D_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}+\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}){\frac{\partial{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}}{\partial
C_{\mu\tilde{t}}}}$
$\displaystyle=2(D_{\tilde{t}\tilde{q}\tilde{r}\tilde{s}}+\Gamma_{\tilde{t}\tilde{q}\tilde{r}\tilde{s}}w_{\tilde{t}\tilde{q}\tilde{r}\tilde{s}})(\mu\tilde{q}|\tilde{r}\tilde{s})$
(42)
$\displaystyle\Gamma_{\tilde{q}\tilde{p}}w_{\tilde{p}\tilde{q}}{\frac{\partial
h^{\prime}_{\tilde{p}\tilde{q}}}{\partial C_{\mu\tilde{t}}}}$
$\displaystyle=2\Gamma_{\tilde{q}\tilde{t}}w_{\tilde{t}\tilde{q}}(t_{\mu
q}+v^{\prime}_{\mu q})C_{q\tilde{q}}$
$\displaystyle-{\frac{1}{2}}\Big{(}\Gamma_{\tilde{q}\tilde{t}}w_{\tilde{t}\tilde{q}}V_{\mu\tilde{q}\tilde{r}\tilde{s}}\Gamma_{\tilde{r}\tilde{s}}+\Gamma_{\tilde{q}\tilde{t}}w_{\tilde{t}\tilde{q}}V_{\tilde{q}\mu\tilde{r}\tilde{s}}\Gamma_{\tilde{r}\tilde{s}}+\Gamma_{\tilde{q}\tilde{r}}w_{\tilde{q}\tilde{r}}V_{\tilde{r}\tilde{q}\mu\tilde{s}}\Gamma_{\tilde{t}\tilde{s}}+\Gamma_{\tilde{q}\tilde{s}}w_{\tilde{q}\tilde{s}}V_{\tilde{s}\tilde{q}\tilde{r}\mu}\Gamma_{\tilde{r}\tilde{t}}\Big{)}$
$\displaystyle=2\Gamma_{\tilde{q}\tilde{t}}w_{\tilde{t}\tilde{q}}(t_{\mu
q}+V_{\mu
qrs}\Gamma^{\prime}_{rs})C_{q\tilde{q}}-\Gamma_{\tilde{q}\tilde{t}}w_{\tilde{q}\tilde{t}}V_{\mu\tilde{q}\tilde{r}\tilde{s}}\Gamma_{\tilde{r}\tilde{s}}-\Gamma_{\tilde{q}\tilde{t}}V_{\mu\tilde{q}\tilde{r}\tilde{s}}\Gamma_{\tilde{r}\tilde{s}}w_{\tilde{r}\tilde{s}}$
(43)

###  VI.3 DMET gradient contribution from the global mean-field

The DMET energy response w.r.t the mean-field solution
($\partial\mathcal{L}_{\mathrm{emb}}/\partial\boldsymbol{C}$) is required to
solve for $\boldsymbol{\alpha}$ (Eq. 14). We label the global occupied MOs by
$i$, $j$ and the virtuals by $a$, $b$. $\boldsymbol{C}$ satisfies

$\displaystyle C_{\mu i}F_{\mu\nu}C_{\nu a}=0$ (44)

Similar to the treatment of the HF solver in section VI.1.1, it is more
convenient to consider
$\partial\mathcal{L}_{\mathrm{emb}}/\partial\gamma_{\mu\nu}\overset{\mathrm{def}}{=\joinrel=}\mathcal{F}_{\mu\nu}$
than $\partial\mathcal{L}_{\mathrm{emb}}/\partial\boldsymbol{C}$. Assuming we
use Löwdin orthogonalization to obtain the LOs, then there is no mean-field
contribution to the AO$\rightarrow$LO transformation gradient. As such,
$h^{\prime}_{\tilde{p}\tilde{q}}$ and
${(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}$ depend on $\gamma_{\mu\nu}$ only
through the LO$\rightarrow$EO transformation $C_{p\tilde{p}}$, which is
obtained from the SVD of $\gamma_{pq}$. On the other hand,
$h_{\tilde{p}\tilde{q}}$ not only depends on $C_{p\tilde{p}}$ due to the basis
transformation but also has the extra “skeleton” dependency on $\gamma$, due
to its explicit appearance in the global Fock matrix used to define
$h_{\tilde{p}\tilde{q}}$ (Eq. II.1).

We first consider the “non-skeleton” part of $\mathcal{F}_{\mu\nu}$, arising
from the LO$\rightarrow$EO transformation, expressed as a chain of partial
derivatives,

$\displaystyle{\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial
C_{\mu\tilde{p}}^{I}}}{\frac{\partial
C_{\mu\tilde{p}}^{I}}{\partial\gamma_{rs}}}{\frac{\partial\gamma_{rs}}{\partial\gamma_{\lambda\sigma}}}={\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial
C_{\mu\tilde{p}}^{I}}}C_{\mu p}{\frac{\partial
C_{p\tilde{p}}^{I}}{\partial\gamma_{rs}}}C_{r\lambda}C_{s\sigma}$ ($r\in I$,
$s,p\centernot\in I$ and $\tilde{p}\in$ bath) (45)

The $r$ and $s$ indices in Eq. VI.3 are restricted because the bath is
constructed from the env-imp block of $\gamma$. The $p$ and $\tilde{p}$
indices are restricted to be in the environment and bath since the other
elements of the LO$\rightarrow$EO transformation matrix are $0$ or $1$ (Eq. 2)
and do not contribute to the derivative. The partial derivative $\partial
C^{I}_{p\tilde{p}}/\partial\gamma_{rs}$ corresponds to the standard expression
for the derivative of the SVD vectors w.r.t the SVD input matrix, and the
prefactor $\partial\mathcal{L}_{\text{emb}}/\partial C_{\mu\tilde{p}}$ in Eq.
VI.3 is given in Appendix VI.2.

As discussed earlier, the “skeleton” part of $\mathcal{F}_{\mu\nu}$ has the
contribution from $h$ only:

$\displaystyle
D_{\tilde{q}\tilde{p}}^{I}h_{\tilde{p}\tilde{q}}^{I[\gamma_{\mu\nu}]}=D_{\tilde{q}\tilde{p}}^{I}V_{\tilde{p}\tilde{q}\mu\nu}^{I}-D_{\tilde{q}\tilde{p}}^{I}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}C_{\tilde{r}\mu}^{I}C_{\tilde{s}\nu}^{I}$
(46)

Adding Eq. VI.3 and Eq. 46 completes the DMET “Fock” matrix
$\mathcal{F}_{\mu\nu}$, and when combined with the CP-HF equation, we obtain
the third term of the DMET gradient (Eq. 16) as

$\displaystyle-\boldsymbol{\alpha}_{0}^{\dagger}{\frac{\partial\boldsymbol{M}}{\partial\boldsymbol{A}}}{\frac{{\mathrm{d}}\boldsymbol{A}}{{\mathrm{d}}{\boldsymbol{R}}}}$
$\displaystyle=2\mathcal{F}_{ai}\mathcal{A}^{-1}_{ai,bj}B^{\boldsymbol{R}}_{bj}-\mathcal{F}_{ij}S_{ij}^{[{\boldsymbol{R}}]}$
(47)

where $S^{[{\boldsymbol{R}}]}$ is the “skeleton” nuclear gradient of the MO
overlap matrix, i.e., $S_{\mu\nu}^{\boldsymbol{R}}C_{\mu i}C_{\nu j}$, and the
expressions for $\mathcal{A}$ and $B^{\boldsymbol{R}}$ can be found in the CP-
HF literature, such as Ref. Yamaguchi  _et al._ (1994).

###  VI.4 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:

$\displaystyle{\frac{\partial\mathcal{L}_{\text{emb}}}{\partial\boldsymbol{A}}}\boldsymbol{A}^{\boldsymbol{R}}={\frac{\partial\mathcal{L}_{\text{emb}}}{\partial\boldsymbol{\mathcal{I}}}}\boldsymbol{\mathcal{I}}^{\boldsymbol{R}}+{\frac{\partial\mathcal{L}_{\text{emb}}}{\partial
S_{\mu\nu}}}S_{\mu\nu}^{\boldsymbol{R}}$ (48)

The first term comes from computing the embedding Lagrangian using AO integral
derivatives ($t_{\mu\nu}^{\boldsymbol{R}}$ and
${(\mu\nu|\lambda\sigma)}^{\boldsymbol{R}}$) in place of the original
integrals:

$\displaystyle{\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial\boldsymbol{\mathcal{I}}}}\boldsymbol{\mathcal{I}}^{\boldsymbol{R}}=\Gamma_{\tilde{q}\tilde{p}}^{I}w_{\tilde{p}\tilde{q}}^{I}h_{\tilde{p}\tilde{q}}^{\prime
I[{\boldsymbol{R}}]}+{\frac{1}{2}}\Big{(}w_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}\Gamma_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}+D_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}\Big{)}\big{[}{(\tilde{p}\tilde{q}|\tilde{r}\tilde{s})}^{I}\big{]}^{[{\boldsymbol{R}}]}+D_{\tilde{q}\tilde{p}}^{I}h_{\tilde{p}\tilde{q}}^{I[{\boldsymbol{R}}]}$
(49)

Regarding the second term in Eq. 48, the Lagrangian $\mathcal{L}_{\text{emb}}$
depends on $S_{\mu\nu}$ in two ways: (1) from the explicit appearance of
$S_{\mu\nu}$ in $h_{\tilde{p}\tilde{q}}$ when transforming the AO 1-RDM into
the EO basis using $C_{\tilde{p}\mu}=S_{\mu\nu}C_{\nu\tilde{p}}$, (2) and the
implicit dependency on $S_{\mu\nu}$ of the LOs. The first contribution, which
can be regarded as a “skeleton” derivative, is simply

$\displaystyle
D_{\tilde{q}\tilde{p}}^{I}h_{\tilde{p}\tilde{q}}^{I[S_{\mu\nu}]}S_{\mu\nu}^{\boldsymbol{R}}=-2D_{\tilde{q}\tilde{p}}^{I}V_{\tilde{p}\tilde{q}\tilde{r}\tilde{s}}^{I}\gamma_{\mu\lambda}C_{\tilde{s}\lambda}^{I}C_{\nu\tilde{r}}^{I}S_{\mu\nu}^{\boldsymbol{R}}$
(50)

The “non-skeleton” contribution can be decomposed further into two components
related to $C_{\mu p}$ and $C_{p\tilde{p}}$ respectively:

$\displaystyle{\frac{\partial\mathcal{L}_{\text{emb}}}{\partial
C_{p\tilde{p}}^{I}}}{\frac{\partial C_{p\tilde{p}}^{I}}{\partial
S_{\mu\nu}}}S_{\mu\nu}^{\boldsymbol{R}}+{\frac{\partial\mathcal{L}_{\text{emb}}}{\partial
C_{\mu p}}}{\frac{\partial C_{\mu p}}{\partial
S_{\lambda\sigma}}}S_{\lambda\sigma}^{\boldsymbol{R}}$ (51)

of which the $C_{p\tilde{p}}$ part of gradient can be computed from the chain
rule,

$\displaystyle{\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial
C_{p\tilde{p}}^{I}}}{\frac{\partial
C_{p\tilde{p}}^{I}}{\partial\gamma_{rs}}}{\frac{\partial\gamma_{rs}}{\partial
S_{\mu\nu}}}S_{\mu\nu}^{\boldsymbol{R}}={\frac{\partial\mathcal{L}_{\text{emb}}}{\partial
C_{\lambda\tilde{p}}^{I}}}C_{\lambda p}{\frac{\partial
C_{\lambda\tilde{p}}^{I}}{\partial\gamma_{rs}}}\Big{(}{\frac{{\mathrm{d}}C_{r\mu}}{{\mathrm{d}}{\boldsymbol{R}}}}C_{s\nu}+C_{s\mu}{\frac{{\mathrm{d}}C_{r\nu}}{{\mathrm{d}}{\boldsymbol{R}}}}\Big{)}\gamma_{\mu\nu}$
(52)

where the prefactor before the parentheses is the same as the prefactor in Eq.
VI.3 and the same orbital restrictions are also applied here. The $C_{p\mu}$
gradient is related to the AO$\rightarrow$LO gradient via

$\displaystyle{\frac{{\mathrm{d}}C_{p\nu}}{{\mathrm{d}}{\boldsymbol{R}}}}={\frac{{\mathrm{d}}C_{\mu
p}}{{\mathrm{d}}{\boldsymbol{R}}}}S_{\nu\mu}+C_{\mu
p}S_{\nu\mu}^{\boldsymbol{R}}={\frac{\partial C_{\mu p}}{\partial
S_{\lambda\sigma}}}S_{\lambda\sigma}^{\boldsymbol{R}}S_{\nu\mu}+C_{\mu
p}S_{\nu\mu}^{\boldsymbol{R}}$ (53)

Since we assume we are using Löwdin orthogonalization, $C_{\mu p}$ is a pure
function of the $S_{\mu\nu}$ eigenvalues and eigenvectors, and thus $\partial
C_{\mu p}/\partial S_{\lambda\sigma}$ can be rewritten using well-known
expressions for matrix eigenvalue and eigenvector derivatives. Again using the
chain rule, the $C_{\mu p}$-related overlap response (second term in Eq. 51)
can be written as,

$\displaystyle{\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial C_{\mu
p}}}{\frac{\partial C_{\mu p}}{\partial
S_{\lambda\sigma}}}S_{\lambda\sigma}^{\boldsymbol{R}}={\frac{\partial\mathcal{L}_{\mathrm{emb}}}{\partial
C_{\mu\tilde{p}}^{I}}}C_{p\tilde{p}}^{I}{\frac{\partial C_{\mu p}}{\partial
S_{\lambda\sigma}}}S_{\lambda\sigma}^{\boldsymbol{R}}$ (54)

where the $\mathcal{L}_{\text{emb}}$ derivative w.r.t. $C_{\mu\tilde{p}}$ was
already discussed in Appendix VI.2. Summing up Eqs. 50, 52 and 54 gives the
full ${\frac{\partial\mathcal{L}_{\text{emb}}}{\partial
S_{\mu\nu}}}S_{\mu\nu}^{\boldsymbol{R}}$, and adding
${\frac{\partial\mathcal{L}_{\text{emb}}}{\partial\boldsymbol{\mathcal{I}}}}\boldsymbol{\mathcal{I}}^{\boldsymbol{R}}$
(Eq. 49) completes the
${\frac{\partial\mathcal{L}_{\text{emb}}}{\partial\boldsymbol{A}}}\boldsymbol{A}^{\boldsymbol{R}}$
part of the DMET gradient in Eq. 16. The last piece, related to the mean-field
response, is given in Eq. 47.

## References

  * Sun and Chan (2016) Q. Sun and G. K.-L. Chan, Accounts of chemical research 49, 2705 (2016).
  * Lee  _et al._ (2019a) S. J. Lee, M. Welborn, F. R. Manby, and T. F. Miller III, Accounts of chemical research 52, 1359 (2019a).
  * Jones  _et al._ (2020) L. O. Jones, M. A. Mosquera, G. C. Schatz, and M. A. Ratner, Journal of the American Chemical Society 142, 3281 (2020).
  * Ma  _et al._ (2021) H. Ma, N. Sheng, M. Govoni, and G. Galli, Journal of Chemical Theory and Computation 17, 2116 (2021).
  * Knizia and Chan (2012) G. Knizia and G. K.-L. Chan, Physical review letters 109, 186404 (2012).
  * Knizia and Chan (2013) G. Knizia and G. K.-L. Chan, Journal of chemical theory and computation 9, 1428 (2013).
  * Cui, Zhu, and Chan (2019) Z.-H. Cui, T. Zhu, and G. K.-L. Chan, Journal of Chemical Theory and Computation 16, 119 (2019).
  * Pham, Hermes, and Gagliardi (2019) H. Q. Pham, M. R. Hermes, and L. Gagliardi, Journal of Chemical Theory and Computation 16, 130 (2019).
  * Hermes and Gagliardi (2019) M. R. Hermes and L. Gagliardi, Journal of chemical theory and computation 15, 972 (2019).
  * Tsuchimochi, Welborn, and Van Voorhis (2015) T. Tsuchimochi, M. Welborn, and T. Van Voorhis, The Journal of chemical physics 143, 024107 (2015).
  * Ye  _et al._ (2018) H.-Z. Ye, M. Welborn, N. D. Ricke, and T. Van Voorhis, The Journal of Chemical Physics 149, 194108 (2018).
  * Lee  _et al._ (2019b) S. J. Lee, F. Ding, F. R. Manby, and T. F. Miller III, The Journal of Chemical Physics 151, 064112 (2019b).
  * Plekhanov, Bonini, and Weber (2021) E. Plekhanov, N. Bonini, and C. Weber, Physical Review B 104, 235131 (2021).
  * Wang  _et al._ (2014) L.-P. Wang, A. Titov, R. McGibbon, F. Liu, V. S. Pande, and T. J. Martínez, Nature chemistry 6, 1044 (2014).
  * Zeng  _et al._ (2020) J. Zeng, L. Cao, M. Xu, T. Zhu, and J. Z. Zhang, Nature communications 11, 1 (2020).
  * Ashraf and Van Duin (2017) C. Ashraf and A. C. Van Duin, The Journal of Physical Chemistry A 121, 1051 (2017).
  * Schlögl (2015) R. Schlögl, Angewandte Chemie International Edition 54, 3465 (2015).
  * Senftle  _et al._ (2016) T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng, Y. K. Shin, C. Junkermeier, R. Engel-Herbert, M. J. Janik, H. M. Aktulga, _et al._ , npj Computational Materials 2, 1 (2016).
  * Wouters  _et al._ (2016) S. Wouters, C. A. Jiménez-Hoyos, Q. Sun, and G. K.-L. Chan, Journal of chemical theory and computation 12, 2706 (2016).
  * Löwdin (1950) P.-O. Löwdin, The Journal of Chemical Physics 18, 365 (1950).
  * Knizia (2013) G. Knizia, Journal of chemical theory and computation 9, 4834 (2013).
  * West  _et al._ (2013) A. C. West, M. W. Schmidt, M. S. Gordon, and K. Ruedenberg, The Journal of chemical physics 139, 234107 (2013).
  * Handy and Schaefer III (1984) N. C. Handy and H. F. Schaefer III, The Journal of chemical physics 81, 5031 (1984).
  * Yamaguchi  _et al._ (1994) Y. Yamaguchi, H. F. Schaefer, Y. Osamura, J. Goddard, _et al._ ,  _A new dimension to quantum chemistry: analytic derivative methods in ab initio molecular electronic structure theory_ (Oxford University Press, USA, 1994).
  * Peng (1941) H. Peng, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 178, 499 (1941).
  * Hariharan and Pople (1973) P. C. Hariharan and J. A. Pople, Theoretica chimica acta 28, 213 (1973).
  * Hehre, Stewart, and Pople (1969) W. J. Hehre, R. F. Stewart, and J. A. Pople, The Journal of Chemical Physics 51, 2657 (1969).
  * Fringuelli and Taticchi (2002) F. Fringuelli and A. Taticchi,  _The Diels-Alder reaction: selected practical methods_ (John Wiley & Sons, 2002).
  * Yang, Doubleday, and Houk (2015) Z. Yang, C. Doubleday, and K. Houk, Journal of Chemical Theory and Computation 11, 5606 (2015).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of chemical physics 132, 154104 (2010).
  * Berne, Ciccotti, and Coker (1998) B. J. Berne, G. Ciccotti, and D. F. Coker,  _Classical and quantum dynamics in condensed phase simulations: Proceedings of the International School of Physics_ (World Scientific, 1998).
  * Henkelman and Jónsson (2000) G. Henkelman and H. Jónsson, The Journal of chemical physics 113, 9978 (2000).
  * Henkelman, Uberuaga, and Jónsson (2000) G. Henkelman, B. P. Uberuaga, and H. Jónsson, The Journal of chemical physics 113, 9901 (2000).
  * Larsen  _et al._ (2017) A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, _et al._ , Journal of Physics: Condensed Matter 29, 273002 (2017).
  * Kowalski and Piecuch (2000) K. Kowalski and P. Piecuch, The Journal of Chemical Physics 113, 18 (2000).
  * Zhao and Truhlar (2006) Y. Zhao and D. G. Truhlar, The Journal of chemical physics 125, 194101 (2006).
  * De Grotthuss (2006) C. De Grotthuss, Biochimica et Biophysica Acta (BBA)-Bioenergetics 1757, 871 (2006).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Werner  _et al._ (2020) H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, _et al._ , The Journal of Chemical Physics 152, 144107 (2020).
  * Pfaendtner and Bonomi (2015) J. Pfaendtner and M. Bonomi, Journal of chemical theory and computation 11, 5062 (2015).
  * Li and Swanson (2020) C. Li and J. M. Swanson, The Journal of Physical Chemistry B 124, 5696 (2020).
  * Li and Voth (2021) C. Li and G. A. Voth, Journal of chemical theory and computation 17, 5759 (2021).
  * Voth (2006) G. A. Voth, Accounts of Chemical Research 39, 143 (2006).
  * Watkins  _et al._ (2019) L. C. Watkins, R. Liang, J. M. Swanson, W. F. DeGrado, and G. A. Voth, Journal of the American Chemical Society 141, 11667 (2019).
  * Lee  _et al._ (2018) M. Lee, C. Bai, M. Feliks, R. Alhadeff, and A. Warshel, Proceedings of the National Academy of Sciences 115, 10321 (2018).
  * Maag  _et al._ (2021) D. Maag, T. Mast, M. Elstner, Q. Cui, and T. Kubař, Proceedings of the National Academy of Sciences 118, e2024803118 (2021).
  * Liu  _et al._ (2021) Y. Liu, C. Li, M. Gupta, N. Verma, A. K. Johri, R. M. Stroud, and G. A. Voth, Proceedings of the National Academy of Sciences 118, e2101932118 (2021).
  * Li  _et al._ (2022) C. Li, Z. Yue, S. Newstead, and G. A. Voth, Biophysical Journal (2022).
  * Shavitt and Bartlett (2009) I. Shavitt and R. J. Bartlett,  _Many-body methods in chemistry and physics: MBPT and coupled-cluster theory_ (Cambridge university press, 2009).