<jupyter_start><jupyter_text>From Second Quantization to Equation-of-Motion Coupled-Cluster using SymPy Table of contents 1. [Introduction](Introduction)2. [Second Quantization](Second-Quantization)3. [Normal product](Normal-product)4. [Contraction](Contraction)5. [Wicks theorem](Wicks-theorem)6. [Particle-Hole formalism](Particle-Hole-formalism)6. [Hartree-Fock](Hartree-Fock)7. [Baker-Campbell-Hausdorff](Baker-Campbell-Hausdorff)8. [Coupled-Cluster](Coupled-Cluster)9. [Equation-of-motion Coupled-Cluster](Equation-of-motion-Coupled-Cluster)  IntroductionWe will have hands on tutorial for the derivation of EOM-CCSD amplitudes, one and two particle density matrix. I have developed a symbolic library using [SymPy](https://www.sympy.org/en/index.html) for deriving analytical expressions which can be easily extended to any operator. First, we will derive the fermionic algebra from first quantization and study their properties. $\newcommand{\ket}[1]{\left|{1}\right\rangle}\newcommand{\bra}[1]{\left\langle{1}\right|}\newcommand{\braket}[2]{\left\langle{1}\middle|{2}\right\rangle}$  Second QuantizationIn this section we will derive the fermionic algebra in second quantization representation from first quantization. Once the fermionic algebra is established, we will derive the relation between the its elements$$ \ket{k} = {a^{+}_{k}}\ket{vac} $$The opertor $a^{+}_{k}$ has created an electron on the vaccum state. $$\phi_{k}(1) \longleftrightarrow {a^{+}_{k}}\ket{vac} $$$\ket{vac}$ is an abstract vaccum state and $\braket{vac}{vac} = 1$If there are two electrons the wavefunction in first quantization is $$ \Phi(1,2) = \frac{1}{\sqrt{2}} \begin{vmatrix}\phi_{i}(1) & \phi_{k}(1) \\\phi_{i}(2) & \phi_{k}(2) \end{vmatrix} $$The same in second quantization follows as (keeping normalization included) $$\ket{ik} = {a^{+}_{i}a^{+}_{k}}\ket{vac} $$ As $\Phi(1,2)$ is an antisymmetric function therefore the creation operator in second quantization should respect this antisymmetry. $$\Phi(2,1) = - \Phi(1,2)  \longleftrightarrow {a^{+}_{k}a^{+}_{i}}\ket{vac} = - {a^{+}_{i}a^{+}_{k}}\ket{vac} $$One can see that the antisymmetry requirement is fulfilled by having the anticommutation relation between creation operators.$$a^{+}_{k}a^{+}_{i} + a^{+}_{i}a^{+}_{k} = \{a^{+}_{k} , a^{+}_{i} \} = 0 \tag{QED 1} $$ Now we will repeat the same excercice for anhilation operators. Consider first one-electron system by creating an electron on $i$th orbital using $a_{i}^{+}$. We can define the anhilator operator $a_{i}$, which act on the former state and returns back the system to $\ket{vac}$. $$a_{i} a_{i}^{+} \ket{vac} = \ket{vac}$$ When there is no electron in a orbital $i$, then it should result to zero. $$a_{i} \ket{vac} = 0$$In order to establish the anticommutation of anhilation operators we will utilize anticommutation relation of creation operators.$$ a_{i} a_{k} \ket{ki} = \ket{vac} \ \ \  \& \ \ \  a_{k} a_{i} \ket{ik} = \ket{vac} = -a_{k} a_{i} \ket{ki}$$Therefore, $$ a_{i}a_{k} + a_{k}a_{i}  = \{a_{k}, a_{i}\} = 0 \tag{QED 2}$$ Now we already have two anticommutation relation and will try to derive the third, between creation and anhilation operators. Consider two orbitals $i \ \& \ k$, where $i\neq k$$$a_{i}a_{k}^{+}\ket{i} = a_{i}\ket{ki} = - a_{i}\ket{ik} = - \ket{k}$$Reversing the action, $$a_{k}^{+} a_{i}\ket{i} = \ket{k}$$Comparing both $$a_{k}^{+} a_{i} + a_{i}a_{k}^{+} = \{ a_{i}, a_{k}^{+} \} = 0 $$But when $i = k$, then $$a_{i}^{+} a_{i} + a_{i}a_{i}^{+} = \{ a_{i}, a_{i}^{+} \} = 1 $$Therefore the final expression which captures both sitation is $$a_{i}^{+} a_{k} + a_{k}a_{i}^{+} = \{ a_{i}^{+}, a_{k} \} = \delta_{ik} \tag{QED 3}$$  Last, we need to establish the adjoint relation between anhilation and creation operator. Consider From first quantization $$\braket {\phi_{i}}{\phi_{k}} = \delta_{ik}$$Using second quantization on right hand side where $\dagger$ indicates hermitian conjugate $$\bra{vac}(a_{i}^{+})^{\dagger} a_{k}^{+}\ket{vac} = \delta_{ik}$$It is only true if $$(a_{i}^{+})^{\dagger} = a_{i}  \tag{QED 4}$$Hence, $a_{i}^{+} = a_{i}^{\dagger}$ which means creation operator for $i$th orbital is hermitian conjugate ofanhilation operator for the same orbital. From now on for creation operator we will use  $\dagger$ as superscript instead of $+$  Particle number operatorWe succesfully showed the transformation of basis from first quantization to second quantization by creating fermionic algebra. Now we will do the same for operators. Lets look at number operator for $N$ fermions in $N$ orthogonal basis. The wavefuntion can be written as $$\ket{\Psi_{N}} =  a^{\dagger}_{N}a^{\dagger}_{N-1} ... a^{\dagger}_{2} a^{\dagger}_{1}\ket{vac}$$Define number operator for $k$th orbital $$a^{\dagger}_{k}a_{k}  $$Act this opeartor on $\ket{\Psi_{N}}$$$a^{\dagger}_{k}a_{k} \ket{\Psi_{N}} = a^{\dagger}_{k}a_{k} (a^{\dagger}_{N}a^{\dagger}_{N-1} ... a^{\dagger}_{k} ...  a^{\dagger}_{2} a^{\dagger}_{1}\ket{vac})$$Using eq.(1) and eq.(3), both $a^{\dagger}_{k}a_{k}$ can move next to $a^{\dagger}_{k}$ inside $\Psi_{N}$$$a^{\dagger}_{k}a_{k} \ket{\Psi_{N}} = a^{\dagger}_{N}a^{\dagger}_{N-1} ...a^{\dagger}_{k}a_{k}a^{\dagger}_{k} ...  a^{\dagger}_{2} a^{\dagger}_{1}\ket{vac}  = a^{\dagger}_{N}a^{\dagger}_{N-1} ...a^{\dagger}_{k}...  a^{\dagger}_{2} a^{\dagger}_{1}\ket{vac} = +1\ket{\Psi_{N}}$$ Particle number operator for $N$ particles$$\sum_{i=1}^{N} a^{\dagger}_{i}a_{i} \ket{\Psi_{N}} = N \ket{\Psi_{N}} $$  Parity of permutation$$p q r s \xrightarrow[]{(12)(34)} q p s r$$ The total number of transposition required to decompose any permutation $P$ is known as its parity, $\mathcal{p}$If the total number is even then its called even parityIf the total number is odd then its called odd paritySign of permutation $P$ is defined as $$sgn(P) =  -1^{p}$$Hence,$$sgn(P_{odd}) =  -1 \hspace{2cm} sgn(P_{even}) =  1$$  Normal product Given the product of operators, the normal product is defined as rearrangment of operators such that all the creation operator are moved to left and all the anhilation operator to right, then multiply it by sign of permutation.Example: Consiper $Q = a_{p}^{\dagger}a_{q}a_{r}^{\dagger}$Now, the in order to write normal product of $Q$, one has to permute $a_{q}a_{r}^{\dagger}$ and multiply thw whole by -1 ( $\because$ odd parity permutation).$$n[a_{p}^{\dagger}a_{q}a_{r}^{\dagger}] = -a_{p}^{\dagger}a_{r}^{\dagger}a_{q} $$Excercise: $Q =  a_{s}a_{p}^{\dagger}a_{r}a_{q}^{\dagger}$1) $n[Q]$ 2) $n[n[Q]]$ (nomal product is idempotent operator) Solution: $Q =  a_{s}a_{p}^{\dagger}a_{r}a_{q}^{\dagger}$1) $n[Q]=-a_{q}^{\dagger}a_{p}^{\dagger}a_{r}a_{s}=a_{p}^{\dagger}a_{q}^{\dagger}a_{r}a_{s}$ 2) $n[n[Q]]=n[Q]$ Now, we will compute normal product using SymPy<jupyter_code>from secondquant_gen import AnnihilateFermion, CreateFermion, NO
from sympy import symbols

p,q,r,s = symbols('p,q,r,s')

Q = AnnihilateFermion(s)*CreateFermion(p)*AnnihilateFermion(r)*CreateFermion(q)

display(Q)

n = NO(Q)
nn = NO(n)

display(n,nn)<jupyter_output><empty_output><jupyter_text>Contraction Contraction between two operators is defined as itself minus its normal product. Symbolically ![con1](figures/con1.jpg)Example: Compute contraction of $Q = a_{p}a_{q}$![ex_con2](figures/ex_con1.jpg)Excercise: Compute contraction of the following 1) $a_{p}a_{q}^{\dagger}$2) $a_{p}^{\dagger}a_{q}$3) $a_{p}^{\dagger}a_{q}^{\dagger}$ Solution:![sol_con1](figures/sol_con1.jpg)Now, we will compute normal product using SymPy<jupyter_code>from secondquant_gen import contraction, F, Fd # F = AnnihilateFermion, Fd = CreateFermion
from sympy import symbols


p,q = symbols('p q')

C1 = contraction(F(p),F(q))
C2 = contraction(F(p),Fd(q))
C3 = contraction(Fd(p),F(q))
C4 = contraction(Fd(p),Fd(q))

display(C1,C2,C3,C4)<jupyter_output><empty_output><jupyter_text>Normal products with ContractionsGiven normal product of operators with contraction of operators inside the normal product. Then take all the contractions out in pairs and same order, leave uncontracted operators inside normal product and multiply sign of permutation. Example: ![NPC_ex1](figures/NPC_ex1.jpg) Example: ![NPC_ex2](figures/NPC_ex2.jpg)  Wicks Theorem Given a product of operators then it can be rewritten as normal product of the same plus normal product with all possible contractions inside.  Example: ![wick1](figures/wick1.jpg)Now, we will show use of wicks theorem using SymPy<jupyter_code>from secondquant_gen import wicks, F, Fd
from sympy import symbols, Dummy


p,q,r,s,t,u = symbols('p q r s t u')

E1 = wicks(F(p)*Fd(q)*Fd(r))
display(E1)
E2 = wicks(F(p)*F(q)*Fd(r)*F(s)*Fd(t)*Fd(u))
display(E2)<jupyter_output><empty_output><jupyter_text>Expectation value with $\ket{vac}$Let $Q = M_{1}M_{2}M_{3} ... M_{n-1}M_{n}$ product of $n$ operators ![prop_wick1](figures/prop_wick1.jpg)From this one can easily see the importance of Wick's theorem. Lets do some checks using SymPy Exercise:$Q = a_{p}a_{q}^{\dagger}a_{r}^{\dagger}$$Q = a_{p}a_{q}a_{r}^{\dagger}a_{s}a_{t}^{\dagger}a_{u}^{\dagger}$Compute $\bra{vac}Q\ket{vac}$?Solution:1. Conver $Q$ to normal order using wicks theorem.2. Only terms which are fully contracted will result in non-zero<jupyter_code>from secondquant_gen import wicks, F, Fd
from sympy import symbols


p,q,r,s,t,u = symbols('p q r s t u')

E1 = wicks(F(p)*Fd(q)*Fd(r),keep_only_fully_contracted=True)
display(E1)
E2 = wicks(F(p)*F(q)*Fd(r)*F(s)*Fd(t)*Fd(u),keep_only_fully_contracted=True)
display(E2)<jupyter_output><empty_output><jupyter_text>Particle Hole formalism This procedure makes the evaluation of certain matrix elements easier like $$\bra{\Phi_{0}}\hat{H}\ket{\Phi_{0}}$$$$\bra{\Phi_{i}^{a}}\hat{H}\ket{\Phi_{0}}$$$$\bra{\Phi_{ij}^{ab}}\hat{H}\ket{\Phi_{0}}$$where, $\Phi_{0}$ is Hartree-Fock determinant. Given a Hartree-Fock solution we get occupied and virtual orbitals. Lets assign symbols $(i,j,k,...)$ to occupied and $(a,b,c,...)$ to virtual. We will define new set of operators $$ b_{p} =\left\{	\begin{array}{ll}		a_{p}^{\dagger}  & \mbox{if } p = i \\		a_{p}            & \mbox{if } p = a	\end{array}\right.$$$$ b_{p}^\dagger =\left\{	\begin{array}{ll}		a_{p}            & \mbox{if } p = i \\		a_{p}^{\dagger}  & \mbox{if } p = a	\end{array}\right.$$This makes the HF determinant a Fermi vaccum$$ b_{p}\ket{\Phi_{0}} = 0$$This follows $$ \{ b_{p},b_{q} \} = 0 \hspace{1cm} \{ b_{p}^{\dagger},b_{q}^{\dagger} \} = 0 \hspace{1cm} \{ b_{p},b_{q}^\dagger \}  = \delta_{pq}  $$  Normal product of PHFIt is defined similarly, as rearrangment of operators ($b_{p}$ \& $b_{q}^{\dagger}$) such that all the creation operator are moved to left and all the anhilation operator to right, then multiply it by sign of permutation.Example: Consiper $R = b_{p}^{\dagger}b_{q}b_{r}^{\dagger}$Now, the in order to write normal product of $R$, one has to permute $b_{q}b_{r}^{\dagger}$ and multiply thw whole by -1 ( $\because$ odd parity permutation).$$N[b_{p}^{\dagger}b_{q}b_{r}^{\dagger}] = -b_{p}^{\dagger}b_{r}^{\dagger}b_{q} $$The thing one has to focus is this definition shows different result that $n[ ]$.Lets do an excercise to explain this. We are going to compute normal order of $a_{p}^{\dagger}a_{q}a_{r}^{\dagger}$ for $p=i,q=j,r=a$$$N[a_{i}^{\dagger}a_{j}a_{a}^{\dagger}] = N[b_{i}b_{j}^{\dagger}b_{a}^{\dagger}] = b_{j}^{\dagger}b_{a}^{\dagger}b_{i} = -b_{a}^{\dagger}b_{j}^{\dagger}b_{i}  =  -a_{a}^{\dagger}a_{j}a_{i}^{\dagger}  $$Also compute $$n[a_{i}^{\dagger}a_{j}a_{a}^{\dagger}] = -a_{i}^{\dagger}a_{a}^{\dagger}a_{j} = a_{a}^{\dagger}a_{i}^{\dagger}a_{j} = a_{a}^{\dagger}\delta_{ij} -  a_{a}^{\dagger}a_{j}^{\dagger}a_{i} $$<jupyter_code>from sympy.physics.secondquant import NO, F, Fd # F = AnnihilateFermion, Fd = CreateFermion
from sympy import symbols


i,j = symbols('i j', below_fermi=True)
a = symbols('a', above_fermi=True)

display(NO(Fd(i)*F(j)*Fd(a)))<jupyter_output><empty_output><jupyter_text>Contractions in PHF Same way they are defined but the operators are from PHF ![con_phf](figures/con_phf.jpg)Similarlry the contraction inside normal product, Wicks theorem are defined.  Expectation value with $\ket{\Phi_{0}}$Let $R = M_{1}M_{2}M_{3} ... M_{n-1}M_{n}$ product of $n$ operators ![phf_prop_wick1](figures/phf_prop_wick1.jpg)From this one can easily see the importance of Wick's theorem. Lets do some checks using SymPy and derive Hartree-Fock energy  Electronic Hamiltonian$$H_{el} = \sum_{pq} h_{pq} a_{p}^{\dagger} a_{q} + \frac{1}{4}\sum_{pqrs} v^{pq}_{rs} a_{p}^{\dagger}a_{q}^{\dagger}a_{s}a_{r}$$where, $$h_{pq} = \bra{\phi_{p} (x)}(-\frac{1}{2}\nabla^{2} -\sum_{\alpha} \frac{Z_{\alpha}}{r_{\alpha x}} )\ket{\phi_{q} (x)}$$$$v_{pqrs} = \bra{\phi_{p} \phi_{q}}\ket{\phi_{r} \phi_{s}} = \braket{\phi_{p} \phi_{q}}{\phi_{r}\phi_{s}} - \braket{\phi_{p} \phi_{q}}{\phi_{s}\phi_{r}}$$ Hartree-FockLets find Hartree-Fock energy<jupyter_code>from sympy.physics.secondquant import wicks, F, Fd, NO  # F = AnnihilateFermion, Fd = CreateFermion
from sympy.physics.secondquant import substitute_dummies, contraction
from sympy.physics.secondquant import AntiSymmetricTensor, evaluate_deltas
from sympy import symbols, Rational, Dummy
    
p, q, r, s = symbols('p,q,r,s', cls=Dummy)                                      

h = AntiSymmetricTensor('h', (p,), (q,))                                        
pq = Fd(p)*F(q)                                                                 
v = AntiSymmetricTensor('v', (p, q), (r, s))                                    
pqsr = Fd(p)*Fd(q)*F(s)*F(r)

H = h*pq + Rational(1, 4)*v*pqsr                                                                                                                                              
display(H) # Without summation                                                                    
                                                                                
eq = wicks(H,keep_only_fully_contracted=True)                                   
index_rule = {'below':  'ijklmno','above':  'abcde', 'general': 'pqrs'}
eq = substitute_dummies(eq, new_indices=True, pretty_indices=index_rule)
display(eq)                                                                     

eq = evaluate_deltas(eq)                                                                
display(eq)<jupyter_output><empty_output><jupyter_text>Normal ordered Hamiltonian$$H_{N} = \sum_{pq} h_{pq} N[a_{p}^{\dagger} a_{q}] + \frac{1}{4}\sum_{pqrs} v^{pq}_{rs} N[a_{p}^{\dagger}a_{q}^{\dagger}a_{s}a_{r}]$$One can show that $$H_{el} = H_{N} + E_{HF}$$From now on we will use $H_{N}$ or computing matrix elements, also $H_{N}$ measures the correlation effect  Baker-Campbell-Hausdorff$$e^{-B} A  e^{B} = A + B + [A,B] + \frac{1}{2!}[[A,B],B] + \frac{1}{3!}[[[A,B],B],B] +  ... $$if $A$ is electronic Hamiltonian or normal ordered Hamiltonian and $B$ any excitation operator, then due to the stucture of Hamiltonian the series truncate after fourth expansion. This is the reason we call coupled cluster as coupled cluster.   Coupled-ClusterCoupled Cluster Singles and Doubles (CCSD)$$\Psi_{CCSD} = e^{T} \ket{\Phi_{0}}$$where, $$ T_{1} = \sum_{ia}t_{i}^{a}a^{\dagger}i$$$$ T_{2} = \frac{1}{4}\sum_{ijab}t_{ij}^{ab}a^{\dagger}b^{\dagger}ji$$Equations which need to be solved$$ E_{CCSD}^{cor} = \bra{\Phi_{0}}e^{-T} H_{N} e^{T} \ket{\Phi_{0}}   \tag{Correltion energy}$$  $$ 0 = \bra{\Phi_{i}^{a}}e^{-T} H_{N} e^{T} \ket{\Phi_{0}}   \tag{$T_{1}$ amplitude}$$  $$ 0 = \bra{\Phi_{ij}^{ab}}e^{-T} H_{N} e^{T} \ket{\Phi_{0}}   \tag{$T_{2}$ amplitude}$$  Lets see the power of SymPy now,<jupyter_code>from sympy.physics.secondquant import wicks, F, Fd, NO  # F = AnnihilateFermion, Fd = CreateFermion
from sympy.physics.secondquant import substitute_dummies, contraction
from sympy.physics.secondquant import AntiSymmetricTensor, evaluate_deltas
from sympy.physics.secondquant import PermutationOperator, simplify_index_permutations
from sympy import symbols, Rational, Dummy
import BCH  
    
p, q, r, s = symbols('p,q,r,s', cls=Dummy) 

i,j = symbols('i,j' , below_fermi=True)                                     
a,b = symbols('a,b' , above_fermi=True)

f = AntiSymmetricTensor('f', (p,), (q,))                                        
pq = NO(Fd(p)*F(q))                                                               
v = AntiSymmetricTensor('v', (p, q), (r, s))                                    
pqsr = NO(Fd(p)*Fd(q)*F(s)*F(r))

H = f*pq + Rational(1, 4)*v*pqsr                                                                                                                                              

ccsd = BCH.level(H,"SD")

eq = wicks(ccsd,simplify_kronecker_deltas=True,keep_only_fully_contracted=True)                                   
index_rule = {'below':  'ijklmno','above':  'abcdef', 'general': 'pqrs'}
e_ccsd = substitute_dummies(eq, new_indices=True, pretty_indices=index_rule)                                                        
display(e_ccsd) 
eq = wicks(Fd(i)*F(a)*ccsd,simplify_kronecker_deltas=True,keep_only_fully_contracted=True)   
index_rule = {'below':  'jklmno','above':  'bcdef', 'general': 'pqrs'}
t1 = (substitute_dummies(eq, new_indices=True, pretty_indices=index_rule))        
display(t1) 
eq = wicks(Fd(i)*Fd(j)*F(b)*F(a)*ccsd,simplify_kronecker_deltas=True,keep_only_fully_contracted=True)   
index_rule = {'below':  'klmno','above':  'cdef', 'general': 'pqrs'}
t2 = substitute_dummies(eq, new_indices=True, pretty_indices=index_rule)    
P = PermutList = [PermutationOperator(i,j),PermutationOperator(a,b)]
t2 = simplify_index_permutations(t2,PermutList)
display(t2)<jupyter_output><empty_output><jupyter_text>Equation-of-motion Coupled-ClusterEOM-CC allows to compute excited and open-shell character electronic states. There are many flavor depending on the target state from reference. Energy and amplitude equations of EOM-CC$$\bar{H} = e^{-T} H_{N} e^{T}$$$$\bra{\Phi_{ij..}^{ab..}}\bar{H}\hat{R}\ket{\Phi_{0}} = E_{EOM} \bra{\Phi_{ij..}^{ab..}}\hat{R}\ket{\Phi_{0}}$$$$\bra{\Phi_{ij..}^{ab..}}[\bar{H}-E_{cc},\hat{R}]\ket{\Phi_{0}} = \Delta E_{EOM} \bra{\Phi_{ij..}^{ab..}}\hat{R}\ket{\Phi_{0}}$$Last equation is the one which we will derive using SymPy Example - EOM-IP-CCSD using genralized Davidson method$$\begin{pmatrix}                                                                   \bar{H}_{SS} - E_{cc} & \bar{H}_{SD} \\                                         \bar{H}_{DS}          &  \bar{H}_{DD} - E_{cc}                                \end{pmatrix}  \begin{pmatrix}                                                                   R_{1}\\                                                                         R_{2}                                                                         \end{pmatrix}=\omega                                                                     \begin{pmatrix}                                                                   R_{1}\\                                                                         R_{2}                                                                         \end{pmatrix}$$                                                               and          $$\begin{pmatrix}                                                                   L_{1} & &  L_{2}                                                              \end{pmatrix}       \begin{pmatrix}                                                                   \bar{H}_{SS}-E_{cc} & \bar{H}_{SD} \\                                           \bar{H}_{DS} & \bar{H}_{DD}-E_{cc}                                            \end{pmatrix}                                                                   =                                                                               \omega                                                                          \begin{pmatrix}                                                                   L_{1} & & L_{2},                                                              \end{pmatrix}    $$                                                          Using matrix equations Need right $\sigma$ amplitudes of EOM-IP-CCSD amplitudes and for this we need $$\bra{\Phi_{i}}[\bar{H}-E_{cc},\hat{R_{1}}]\ket{\Phi_{0}} = ((\bar{H}_{SS} - E_{cc}) R_{1})$$$$\bra{\Phi_{ij}^{a}}[\bar{H},\hat{R_{1}}]\ket{\Phi_{0}} = (\bar{H}_{DS} R_{1}) $$$$\bra{\Phi_{i}}[\bar{H},\hat{R_{2}}]\ket{\Phi_{0}} = (\bar{H}_{SD} R_{2})$$$$\bra{\Phi_{ij}^{a}}[\bar{H}-E_{cc},\hat{R_{2}}]\ket{\Phi_{0}} =((\bar{H}_{DD} - E_{cc}) R_{2}) $$Therefore, trial vectors are defined as                                                             $$\sigma_{1}=((\bar{H}_{SS}-E_{cc})R_{1})+(\bar{H}_{SD}R_{2})$$$$\sigma_{2}=(\bar{H}_{DS}R_{1}) +((\bar{H}_{DD}-E_{cc})R_{2})$$Similarly, left vectors are defined. Lets use power of SymPy for right $\sigma$ amplitudes of EOM-IP-CCSD<jupyter_code>import EOM,SIGMA
flavor1 = "IP"
R0_f1 = EOM.R0(flavor1)
R1_f1 = EOM.R1(flavor1)
R2_f1 = EOM.R2(flavor1)
Rf1  = R0_f1 + R1_f1 + R2_f1
                                                                            
SIGMA.RVECTORS(R0_f1,R1_f1,R2_f1,flavor1)

L0_f1 = EOM.L0(flavor1)
L1_f1 = EOM.L1(flavor1)
L2_f1 = EOM.L2(flavor1)
Lf1  = L0_f1 + L1_f1 + L2_f1
                                                                            
SIGMA.LVECTORS(L0_f1,L1_f1,L2_f1,flavor1)<jupyter_output><empty_output><jupyter_text>Properties  One particle density matrix (OPDM)$$\gamma^I_{pq}= \bra{\Psi_{I}}p^{\dagger}q\ket{\Psi_{I}}$$ One particle transition density matrix (OPTDM)$$\gamma^{IJ}_{pq}= \bra{\Psi_{I}}p^{\dagger}q\ket{\Psi_{J}}$$<jupyter_code>import EOM, DM, TDM 
flavor1 = "IP"
R0_f1 = EOM.R0(flavor1)
R1_f1 = EOM.R1(flavor1)
R2_f1 = EOM.R2(flavor1)
Rf1  = R0_f1 + R1_f1 + R2_f1
                                                                                
L0_f1 = EOM.L0(flavor1)
L1_f1 = EOM.L1(flavor1)
L2_f1 = EOM.L2(flavor1)
Lf1  = L0_f1 + L1_f1 + L2_f1
DM.OPDM(Lf1,Rf1,flavor1)
flavor2 = "CCSD"                                                                
R0_f2 = EOM.R0(flavor2)                                                    
R1_f2 = EOM.R1(flavor2)                                                    
R2_f2 = EOM.R2(flavor2)                                                    
Rf2  = R0_f2 + R1_f2 + R1_f2                                                   
                                                                                
L0_f2 = EOM.L0(flavor2)                                                    
L1_f2 = EOM.L1(flavor2)                                                    
L2_f2 = EOM.L2(flavor2)                                                    
Lf2  = L0_f2 + L1_f2 + L2_f2                                                

TDM.OPTDM(Lf1,Rf1,Lf2,Rf2,flavor1,flavor2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17446043165467626, 'fraction_numerical': 0.014607825934374452, 'mean_word_length': 2.6804972554084596, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 117, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '660604', 'license': 'MIT', 'n_tokens_mistral': 7765, 'n_tokens_neox': 7386, 'n_words': 1776, 'path': 'SQ2EOM.ipynb', 'repo_name': 'sgulania/SQ2EOM', 'size': '22796'}