CCSD.tex

\documentclass{ctexart}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{fancyhdr}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage[a4paper]{geometry}
\usepackage{indentfirst}
\usepackage{siunitx}
\usepackage{mhchem}
\usepackage{subfigure}
\usepackage{float}
\usepackage{caption}
\usepackage{multirow}
\usepackage{pifont}
\usepackage{caption}
\usepackage{booktabs}
\usepackage{physics}
\usepackage{color}
\usepackage{url}
\usepackage{listings}
\usepackage[colorlinks=true, citecolor=blue, urlcolor=blue, linkcolor=blue, breaklinks=true, pdfpagelabels=false]{hyperref}
\usepackage[style=chem-acs, backend=biber]{biblatex}

\numberwithin{figure}{section}
\numberwithin{table}{section}

\addbibresource{ref.bib}

\newcommand{\me}{\mathrm{e}}
\newcommand{\mi}{\mathrm{i}}
\newcommand{\p}{\partial}
\newcommand{\mycite}[1]{\textsuperscript{\cite{#1}}}
\newcommand{\m}[1]{\mathrm{#1}}
\newcommand{\bu}{\bullet}

\lstset{breaklines=true, basicstyle=\ttfamily}

\DeclareMathOperator{\arcsinh}{arcsinh}
\DeclareMathOperator{\arccosh}{arccosh}
\DeclareMathOperator{\arctanh}{arctanh}
\DeclareMathOperator{\arccsch}{arccsch}

\geometry{top=2.54cm,bottom=2.54cm,left=2.54cm,right=2.54cm}
\newcommand\msi[2]{\SI[group-digits=false]{#1}{#2}}

\title{Coupled-Cluster Singles and Doubles}
\author{Peter W\thanks{Email: wtpeter@pku.edu.cn}}
\date{202501}

\pagestyle{fancy}
\fancyhf{}
\fancyhead[L]{}
\fancyhead[C]{}
\fancyhead[R]{}
\fancyfoot[C]{\thepage}


\begin{document}
\maketitle

这个笔记给出了CC方法的基本理论，并给出了CCSD方程的具体推导和Python代码实现。

\section{耦合簇方法}

这一部分参考了 Ref.~\cite{helgaker2000molecularelectronicstructure} 中的13.1、13.2和13.4节。

CC的波函数拟设是
\[\ket{\m{CC}}=\me^{\hat{T}}\ket{\m{HF}},\]
其中$\hat{T}$是激发算符
\[\hat{T}=\sum_{\mu} t_\mu \hat \tau_{\mu}=\hat{T}_1+\hat{T}_2+\cdots,\]
$\mu$指各种激发方式。
最常见的单激发和双激发算符是
\[\hat{T}_1=\sum_{ia} t_i^a a_a^\dagger a_i=\sum_{ia} t_i^a \hat\tau_i^a,\quad \hat{T}_2=\frac{1}{4}\sum_{ijab} t_{ij}^{ab}a_a^\dagger a_i a_b^\dagger a_j=\frac{1}{4}\sum_{ijab} t_{ij}^{ab}\hat\tau_{ij}^{ab},\]
其中$i,j,k,l,m$表示占据轨道，$a,b,c,d,e$表示空轨道。
在本文中使用的轨道均是自旋轨道，即GCCSD。
一般将$t$称为激发振幅。

注意到对任意的$\mu$和$\nu$激发算符，不论其阶数和轨道指标都有
\[[\hat\tau_\mu,\hat\tau_\nu]=0,\]
所以实际上
\[\me^{\hat{T}}=\prod_{\mu} (1+t_\mu \hat\tau_\mu),\]
就像Grassmann数一样。
所以CC波函数的形式也可以写成
\[\ket{\m{CC}}=\left[\prod_\mu (1+t_\mu \hat\tau_\mu)\right] \ket{\m{HF}},\]
而CI波函数的形式则是
\[\ket{\m{CI}}=\left(1+\sum_\mu C_\mu\hat\tau_\mu\right)\ket{\m{HF}}.\]
所以在截断的情况下，尽管二者的变分参数数量是一样的，
但是CC可以提供到FCI中所有的激发组态，而CI只有截断阶的激发组态。
在不截断的情况下，Full-CC和Full-CI是等价的。
在平衡位置附近，CCSD和FCI相比，得到的单激发和双激发组态权重只相差大约5\%，而四激发组态也可以恢复FCI大约84\%的结果，
但是几乎得不到三激发组态。
另外，随着截断阶数的增加（S、D、T、Q ……），CC几乎可以按照固定的比例逼近FCI结果，而CI则不行。

CC相比CI还具有大小一致性，并且这种一致性相比MPPT更好，是逐项一致而非逐阶一致的。

CC波函数/能量的优化远比CI更复杂。
对于CI而言，由于波函数的形式是线性组合，对能量求极小值完全等价于向各个组态的投影法，因此对波函数的优化在形式上非常简单。
但是CC波函数的形式是非线性的，投影法并不等价于对能量的极小化，而极小化又由于$\hat{T}$的非厄米（也非反厄米）而非常复杂，
所以实践上还是基本采用投影法，从结果上投影法得到的能量和极小化能量得到的能量非常接近，但是它不再是FCI能量的上界。

在使用投影法计算CC能量和系数时，简单来看是（假设）CC波函数满足
\[H\me^T\ket{\m{HF}}=E\me^{T} \ket{\m{HF}},\]
并将它投影到$\ket{\m{HF}}$和$\ket{\m{\mu}}=\hat \tau_\mu \ket{\m{HF}}$上得到
\[\mel{\m{HF}}{{H}\me^{{T}}}{\m{HF}}=E,\quad \mel{\mu}{H\me^T}{\m{HF}}=E\mel{\mu}{\me^{T}}{\m{HF}}.\]
而更常用的做法是考虑对相似变换后的哈密顿$\me^{-T}H\me^{T}$做投影，得到能量方程
\[\mel{\m{HF}}{\me^{-T}H\me^{T}}{\m{HF}}=E,\]
和振幅方程
\[\mel{\mu}{\me^{-T}H\me^T}{\m{HF}}=0.\]
这种做法的优势是可以使用BCH公式
\[\me^{-T}H\me^{T}=H+[H,T]+\frac{1}{2}[[H,T],T]+\cdots,\]
并利用Slater-Condon规则，大幅简化公式的复杂程度。

% 如果参考态是正则HF基态，那么可以将哈密顿改写为
% \[H=f+\Phi,\]
% 其中$f$是Fock算符
% \[f=\sum_{pq}\left[h_{pq}+\sum_i \mel{pi}{}{qi}\right]a_p^\dagger a_q=\sum_p \varepsilon_p a_p^\dagger a_p,\]
% $\Phi$是涨落算符
% \[\Phi=\frac{1}{2}\sum_{pqrs}\braket{pq}{rs}a_p^\dagger a_q^\dagger a_s a_r-\sum_i \mel{pi}{}{qi}a_p^\dagger a_q=\frac{1}{4}\sum_{pqrs}\mel{pq}{}{rs}a_p^\dagger a_q^\dagger a_s a_r-\sum_i \mel{pi}{}{qi}a_p^\dagger a_q.\]
% 可以看到
% \[[f,\tau_\mu]=\varepsilon_\mu t_\mu \tau_\mu,\]
% 其中$\varepsilon_\mu$是轨道能量差，对单双激发的情况下是
% \[\varepsilon_{ia}=\varepsilon_a-\varepsilon_i,\quad \varepsilon_{ijab}=\varepsilon_{a}+\varepsilon_b-\varepsilon_i-\varepsilon_j.\]
% 因此可以将投影方程写成
% \begin{align}
%     E&=E_{\m{HF}}+\mel{\m{HF}}{\Phi(T_2+\frac{1}{2}T_1^2)}{\m{HF}}, \label{energy1}\\
%     \varepsilon_\mu t_\mu&=-\mel{\mu}{\me^{-T}\Phi\me^T}{\m{HF}}, \label{amplitude1}
% \end{align}
% 其中$E_{\m{HF}}$是HF基态能量，第一个式子中没有高于二次的项是由于$\Phi$至多只有双体算符，
% 而没有$T_1$的一次项是因为$\mel{\m{HF}}{\Phi T_1}{\m{HF}}=0$（Slater-Condon规则）。

% 为了优化振幅，一般采用赝牛顿法，即直接将Eq.~(\ref{amplitude1})的右侧计算出来然后除以能量差，即在第$n$次迭代中采用
% \[\Delta t_\mu^{(n)}=-\varepsilon_\mu^{-1} \mel{\mu}{\me^{-T^{(n)}}\Phi\me^{T^{(n)}}}{\m{HF}},\]
% 其中$T^{(n)}$是由$t^{(n)}$计算得到的激发算符。


\section{产生湮灭算符期望值的计算}

这一部分参考了 Ref.~\cite{crawford2000introductioncoupled}中的55-63页。

为了计算一串产生湮灭算符在HF基态上的期望值，在此引入Wick定理。
可以先将算符分成准粒子的产生和湮灭算符，它们的定义是作用在某个态上产生一个准粒子或湮灭一个准粒子。
对于HF基态$\ket{\m{HF}}$来说，准粒子是电子或者空穴，这时准粒子产生算符是$a_a^\dagger$和$a_i$，
湮灭算符是$a_a$和$a_i^\dagger$，满足准粒子湮灭算符作用在HF基态上为零
\[a_a\ket{\m{HF}}=a_i^\dagger \ket{\m{HF}}=0.\]
接下来引入算符正规序的概念，它要求将算符串$ABC$通过对换排列成准粒子的产生算符在左侧，准粒子湮灭算符在右侧，将符号记为$N(ABC\cdots)$，
对费米子算符来说还需要在每个对换中加上一个负号。
下一个与正规序相关的概念是收缩，是指两个算符串本身减去它们的正规序，记做
\[A^{\bu}B^{\bu}=AB-N(AB).\]
根据收缩的定义，可以看到有只有两种不为零的收缩
\[a_a^{\bu} a_b^{\dagger \bu}=a_a a_b^\dagger - N(a_a a_b^\dagger)=a_a a_b^\dagger + a_b^\dagger a_a=\delta_{ab},\]
和
\[a_i^{\dagger \bu} a_j^\bu=a_i^\dagger a_j-N(a_i^\dagger a_j)=a_i^\dagger a_j+a_j a_i^\dagger=\delta_{ij}.\]
这里要注意，算符对在收缩之后会变成$c$-数（与其他所有符号都对易），就不再是算符了。

有了这些准备，可以引入Wick定理，它的表述是
\begin{align*}
    ABC\cdots XYZ&=N(ABC\cdot XYZ)\\
    &\quad +N(A^{\bu}B^{\bu}C\cdots XYZ)+N(A^{\bu}BC^{\bu}\cdots XYZ)+(\m{singles})\\
    &\quad +N(A^{\bu}B^{\bu}C^{\bu\bu}D^{\bu\bu}\cdot XYZ)+N(A^{\bu}B^{\bu}C^{\bu\bu}DE^{\bu\bu}\cdots XYZ)+(\m{doubles})\\
    &\quad +\cdots,
\end{align*}
即算符串可以写为所有可能的收缩的正规序之和，这些收缩包括单收缩（一对）、双收缩（两对）等等。
此外，还有推广的Wick定理，它是
\begin{align*}
    N(ABC\cdots)N(XYZ\cdots)&=N(ABC\cdots XYZ\cdots)\\
    &\quad +N(A^{\bu}BC\cdots X^{\bu}YZ\cdots)+N(A^{\bu}BC\cdots XY^{\bu}Z\cdots)+(\m{singles})\\
    &\quad +N(A^{\bu}B^{\bu\bu}C\cdots X^{\bu}Y^{\bu\bu}Z\cdots)+N(A^{\bu}B^{\bu\bu}C\cdots X^{\bu}YZ^{\bu\bu}\cdots)+(\m{doubles})\\
    &\quad +\cdots,
\end{align*}
即收缩要在两个不同组的正规序之间进行。
Wick定理的好处是在求算符串求在HF基态的期望值时，所有含有算符的正规序项全部为零，只有完全收缩项非零。

接下来可以将哈密顿写成正规序形式，在二次量子化下的哈密顿是
\[H=\sum_{pq}h_{pq}+\frac{1}{4}\sum_{pqrs}\mel{pq}{}{rs}a_p^\dagger a_q^\dagger a_s a_r,\]
其中$\mel{pq}{}{rs}$是反对称双电子积分，这里用$p,q,r,s$表示任意轨道。
如果用Wick定理将它变成正规序，可以得到
\[H=\sum_{pq}f_{pq}N(a_p^\dagger a_q)+\frac{1}{4}\sum_{pqrs}\mel{pq}{}{rs}N(a_p^\dagger a_q^\dagger a_s a_r)+E_{\m{HF}}=F_N+V_N+E_{\m{HF}},\]
其中下标$N$代表正规序，接下来将采用如下形式的正规序哈密顿量来进行CC的相关计算
\[H_N=H-\mel{\m{HF}}{H}{\m{HF}}=F_N+V_N.\]
注意一个算符的正规序是它减去它的期望值，这个结论也可以被推广。

同时，可以将激发算符也写成正规序$T_N$的形式，例如单激发算符
\[T_{1N}=\sum_{ia}t_i^a a_a^\dagger a_i=\sum_{ia} t_i^a N(a_a^\dagger a_i),\]
和双激发算符
\[T_{2N}=\sum_{ijab}t_{ij}^{ab}a_a^\dagger a_i a_b^\dagger a_j=\sum_{ijab}t_{ij}^{ab}a_a^\dagger a_b^\dagger a_ja_i=\sum_{ijab}t_{ij}^{ab}N(a_a^\dagger a_b^\dagger a_ja_i).\]
这里利用了正规序是算符串本身减去期望值（零），可以自然地加上正规序符号。


\section{CCSD方法}

CCSD方法即是将激发算符只保留单双激发，$T=T_1+T_2=T_{1N}+T_{2N}=T_N$。
而能量方程和振幅方程中的算符部分先写成
\[E_{\m{corr}}=\mel{\m{HF}}{\me^{-T_N}H_N\me^{T_N}}{\m{HF}},\quad 0=\mel{\mu}{\me^{-T_N}H_N\me^{T_N}}{\m{HF}},\]
其中能量方程扣除了常数的HF基态能量。
可以证明在BCH公式中只会截断到四次对易，即
\[\me^{-T_N}H_N\me^{T_N}=H_N+[H_N,T_N]+\frac{1}{2}[[H_N,T_N],T_N]+\frac{1}{6}[[[H_N,T_N],T_N],T_N]+\frac{1}{24}[[[[H_N,T_N],T_N],T_N],T_N].\]

利用SymPy可以推导上述公式，代码由 \url{https://github.com/sympy/sympy/blob/1.12/examples/intermediate/coupled_cluster.py} 修改而来，见 SymPyCC.ipynb。
其中使用的哈密顿量为正则轨道下的$H_N$（轨道能量为$\varepsilon_p$），在对$ijab$以外的指标使用Einstein求和约定、用$\{AB\cdots\}$表示正规序$N(AB\cdots)$的情况下为
\[H_N={f^{p}_{p}} \left\{{a^\dagger_{p}} a_{p}\right\} - \frac{{v^{pq}_{rs}} \left\{{a^\dagger_{p}} {a^\dagger_{q}} a_{r} a_{s}\right\}}{4},\quad f_p^p=\varepsilon_{p},\quad v_{rs}^{pq}=\mel{pq}{}{rs}.\]
使用符号计算可以得到能量公式
\[E_{\m{corr}}=- \frac{{t^{c}_{l}} {t^{d}_{k}} {v^{kl}_{cd}}}{2} + \frac{{t^{cd}_{kl}} {v^{kl}_{cd}}}{4}.\]
以及$T_1$振幅方程
\[0={f^{a}_{a}} {t^{a}_{i}} - {f^{i}_{i}} {t^{a}_{i}} - {t^{c}_{k}} {t^{d}_{i}} {t^{a}_{l}} {v^{kl}_{cd}} - {t^{c}_{k}} {t^{d}_{i}} {v^{ak}_{cd}} + {t^{c}_{k}} {t^{a}_{l}} {v^{kl}_{ic}} + {t^{c}_{k}} {t^{ad}_{il}} {v^{kl}_{cd}} + {t^{c}_{k}} {v^{ak}_{ic}} - \frac{{t^{c}_{i}} {t^{ad}_{kl}} {v^{kl}_{cd}}}{2} + \frac{{t^{a}_{l}} {t^{cd}_{ik}} {v^{kl}_{cd}}}{2} + \frac{{t^{cd}_{ik}} {v^{ak}_{cd}}}{2} - \frac{{t^{ac}_{kl}} {v^{kl}_{ic}}}{2}.\]
还有$T_2$振幅方程
\begin{align*}
0&={f^{a}_{a}} {t^{ab}_{ij}} P(ab) - {f^{i}_{i}} {t^{ab}_{ij}} P(ij) + {t^{c}_{k}} {t^{d}_{i}} {t^{ab}_{jl}} {v^{kl}_{cd}} P(ij) + {t^{c}_{k}} {t^{a}_{l}} {t^{bd}_{ij}} {v^{kl}_{cd}} P(ab) - {t^{c}_{k}} {t^{ad}_{ij}} {v^{bk}_{cd}} P(ab) + {t^{c}_{k}} {t^{ab}_{il}} {v^{kl}_{jc}} P(ij)\\
&\quad + {t^{c}_{i}} {t^{d}_{j}} {t^{a}_{k}} {t^{b}_{l}} {v^{kl}_{cd}} + {t^{c}_{i}} {t^{d}_{j}} {t^{a}_{k}} {v^{bk}_{cd}} P(ab) + \frac{{t^{c}_{i}} {t^{d}_{j}} {t^{ab}_{kl}} {v^{kl}_{cd}}}{2} + {t^{c}_{i}} {t^{d}_{j}} {v^{ab}_{cd}} - {t^{c}_{i}} {t^{a}_{k}} {t^{b}_{l}} {v^{kl}_{jc}} P(ij) - {t^{c}_{i}} {t^{a}_{k}} {t^{bd}_{jl}} {v^{kl}_{cd}} P(ab) P(ij)\\
&\quad - {t^{c}_{i}} {t^{a}_{k}} {v^{bk}_{jc}} P(ab) P(ij) - {t^{c}_{i}} {t^{ad}_{jk}} {v^{bk}_{cd}} P(ab) P(ij) - \frac{{t^{c}_{i}} {t^{ab}_{kl}} {v^{kl}_{jc}} P(ij)}{2} - {t^{c}_{i}} {v^{ab}_{jc}} P(ij) + \frac{{t^{a}_{k}} {t^{b}_{l}} {t^{cd}_{ij}} {v^{kl}_{cd}}}{2} + {t^{a}_{k}} {t^{b}_{l}} {v^{kl}_{ij}}\\
&\quad + \frac{{t^{a}_{k}} {t^{cd}_{ij}} {v^{bk}_{cd}} P(ab)}{2} + {t^{a}_{k}} {t^{bc}_{il}} {v^{kl}_{jc}} P(ab) P(ij) + {t^{a}_{k}} {v^{bk}_{ij}} P(ab) + \frac{{t^{cd}_{ij}} {t^{ab}_{kl}} {v^{kl}_{cd}}}{4} + \frac{{t^{cd}_{ij}} {v^{ab}_{cd}}}{2} + \frac{{t^{cd}_{jk}} {t^{ab}_{il}} {v^{kl}_{cd}} P(ij)}{2}\\
&\quad - \frac{{t^{ac}_{kl}} {t^{bd}_{ij}} {v^{kl}_{cd}} P(ab)}{2} + {t^{ac}_{ik}} {t^{bd}_{jl}} {v^{kl}_{cd}} P(ab) + {t^{ac}_{ik}} {v^{bk}_{jc}} P(ab) P(ij) + \frac{{t^{ab}_{kl}} {v^{kl}_{ij}}}{2} + {v^{ab}_{ij}},
\end{align*}
其中$P(pq)$是对称化操作符，即$f(p,q)P(pq)=f(p,q)-f(q,p)$。

接下来对上述方程进行一些简单的化简和符号变换，把求和显式写出，并将反对称化双电子积分回归到常见的$\mel{pq}{}{rs}$的形式。
首先能量方程变为
\[E_{\m{corr}}=\frac{1}{2}\sum_{ijab}t_i^a t_j^b \mel{ij}{}{ab}+\frac{1}{4}\sum_{ijab}t_{ij}^{ab}\mel{ij}{}{ab}.\]
$T_1$振幅方程变为
\begin{align*}
    D_i^a t_i^a&=\sum_{kc}t_k^c\mel{ak}{}{ic}-\sum_{kcd}t_k^ct_i^d\mel{ak}{}{cd}+\sum_{klc}t_k^ct_l^a \mel{kl}{}{ic}-\sum_{kcdl}t_k^c t_i^d t_l^a\mel{kl}{}{cd}+\frac{1}{2}\sum_{kcd}t_{ik}^{cd}\mel{ak}{}{cd}\\
    &\quad -\frac{1}{2}\sum_{klc}t_{kl}^{ac}\mel{kl}{}{ic} +\sum_{klcd}t_k^c t_{il}^{ad}\mel{kl}{}{cd}-\frac{1}{2}\sum_{klcd}t_i^ct_{kl}^{ad}\mel{kl}{}{cd}+\frac{1}{2}\sum_{klcd}t_l^a t_{ik}^{cd}\mel{kl}{}{cd},
\end{align*}
其中$D_i^a=\varepsilon_i-\varepsilon_a$。
$T_2$振幅方程变为
\begin{align*}
    D_{ij}^{ab}t_{ij}^{ab}&=\mel{ab}{}{ij}+\frac{1}{2}\sum_{kl}t_{kl}^{ab}\mel{kl}{}{ij}+\frac{1}{2}\sum_{cd}t_{ij}^{cd}\mel{ab}{}{cd}+P(ab)P(ij)\sum_{kc}t_{ik}^{ac}\mel{bk}{}{jc}\\
    &\quad +P(ab)\sum_{k} t_k^a\mel{bk}{}{ij}-P(ij)\sum_{j}t_i^c \mel{ab}{}{jc}+P(ab)P(ij)\sum_{klc}t_k^a t_{il}^{bc}\mel{kl}{}{jc}\\
    &\quad -P(ab)P(ij)\sum_{kcd}t_i^c t_{jk}^{ad}\mel{bk}{}{cd}-P(ab)\sum_{kcd}t_k^ct_{ij}^{ad}\mel{bk}{}{cd}+P(ij)\sum_{klc}t_k^ct_{il}^{ab}\mel{kl}{}{jc}\\
    &\quad +\sum_{kl}t_k^a t_l^b \mel{kl}{}{ij}+\sum_{cd}t_i^ct_j^d\mel{ab}{}{cd}-P(ab)P(ij)\sum_{kc}t_i^ct_k^a\mel{bk}{}{jc}\\
    &\quad +\frac{1}{4}\sum_{klcd}t_{ij}^{cd}t_{kl}^{ab}\mel{kl}{}{cd}-\frac{1}{2}P(ab)\sum_{klcd}t_{kl}^{ac}t_{ij}^{bd}\mel{kl}{}{cd}+\frac{1}{2}P(ij)\sum_{klcd}t_{jk}^{cd}t_{il}^{ab}\mel{kl}{}{cd}\\
    &\quad +\frac{1}{2}P(ab)P(ij)\sum_{klcd}t_{ik}^{ac}t_{jl}^{bd}\mel{kl}{}{cd}+\frac{1}{2}P(ab)\sum_{kcd}t_k^a t_{ij}^{cd} \mel{bk}{}{cd}-\frac{1}{2}P(ij)\sum_{klc}t_i^ct_{kl}^{ab}\mel{kl}{}{jc}\\
    &\quad -P(ij)\sum_{klc}t_i^c t_k^at_l^b \mel{kl}{}{jc}+P(ab)\sum_{kcd}t_i^c t_j^d t_k^a \mel{bk}{}{cd}+\frac{1}{2}\sum_{klcd}t_i^c t_j^d t_{kl}^{ab} \mel{kl}{}{cd}\\
    &\quad +\frac{1}{2}\sum_{klcd}t_k^a t_l^b t_{ij}^{cd}\mel{kl}{}{cd}+P(ij)\sum_{klcd}t_k^c t_i^d t_{jl}^{ab}\mel{kl}{}{cd}+P(ab)\sum_{klcd}t_k^c t_l^a t_{ij}^{bd} \mel{kl}{}{cd}\\
    &\quad -P(ab)P(ij)\sum_{klcd}t_i^c t_k^a t_{jl}^{bd}\mel{kl}{}{cd}+\sum_{klcd}t_i^c t_j^d t_k^a t_l^b \mel{kl}{}{cd},
\end{align*}
其中$D_{ij}^{ab}=\varepsilon_i+\varepsilon_j-\varepsilon_a-\varepsilon_b$。
这三个方程即是正则轨道下，
Ref~\cite{crawford2000introductioncoupled}中的[134]、[152]和[153]式，其中[152]式中第二行的最后一项应反号；
两个振幅方程也是Ref~\cite{stanton1991directproduct}中的(1)和(2)式，其中(2)式第一项应为复共轭。

在优化振幅时，一般采用赝牛顿法，即直接将上述方程的右侧计算出来然后除以能量差$D$，直到振幅和能量达到收敛标准。
大部分情况下可以使用MP2振幅作为CCSD振幅的初猜，即
\[t_i^a=0,\quad t_{ij}^{ab}=\frac{\mel{ab}{}{ij}}{\varepsilon_i+\varepsilon_j-\varepsilon_a-\varepsilon_b},\]
这相当于假设上述方程的第$-1$次迭代使用的振幅$t^{(-1)}$是全零的。


\section{DIIS}

CC方法也可以在上述的赝牛顿法基础上加入DIIS加速收敛，即取
\[t^{(n+1)}=\sum_k^M w_k t^{(k)},\]
其中$M$是子空间大小，权重$w_k$满足
\[\sum_k^M w_k=1.\]
一般取权重的方式是使迭代步长极小化，即
\[\Delta t^{\m{ave}}=\sum_{k}^M w_k \Delta t^{(k)}.\]
因此误差向量是$\vb*{e}_k=\Delta t^{(k)}$，可以简单地直接将$\Delta t_i^a$和$\Delta t_{ij}^{ab}$平铺并拼接成一个向量，
然后求解
\[\begin{pmatrix}
    B_{11} & B_{12} & \cdots & B_{1M} & -1\\
    B_{21} & B_{22} & \cdots & B_{2M} & -1\\
    \vdots & \vdots & \ddots & \vdots & \vdots\\
    B_{M1} & B_{M2} & \cdots & B_{MM} & -1\\
    -1 & -1 & \cdots & -1 & 0
\end{pmatrix}\begin{pmatrix}
    w_1\\
    w_2\\
    \vdots\\
    w_M\\
    \lambda
\end{pmatrix}=\begin{pmatrix}
    0\\
    0\\
    \vdots\\
    0\\
    -1
\end{pmatrix},\]
可以得到权重，其中$B_{ij}=\vb*{e}_i\cdot \vb*{e}_j$，$\lambda$是拉格朗日乘子。


\section{实现}

CCSD部分的代码都在CCSD.ipynb。
使用了PySCF进行HF参考态的计算并获得双电子积分，之后将RHF结果变为GHF的自旋轨道形式，
其中\texttt{kernel()}和\texttt{kernel\_diis()}可以进行CCSD和CCSD+DIIS的计算。



\printbibliography[title={参考文献}]

\end{document}