From 04704ba6256c7d1e250b6adc1e1dec831b6cb555 Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Thu, 9 May 2024 16:12:36 -0400
Subject: [PATCH 1/7] Linear response summetric

---
 pycc/ccresponse.py | 190 +++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 183 insertions(+), 7 deletions(-)

diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py
index 9f51250..0f34b70 100644
--- a/pycc/ccresponse.py
+++ b/pycc/ccresponse.py
@@ -296,6 +296,182 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis
                     check.append(polar)
 
 
+    def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
+        """
+        	Calculate the CC linear response function for polarizability at field-frequency omega(w1).
+
+        	The linear response function, <<A;B(w1)>> generally reuires the following perturbed wave functions and frequencies:
+
+        	Parameters
+        	----------
+        	pertkey_a: string
+        		String identifying the one-electron perturbation, A along a cartesian axis
+
+        	Return
+        	------
+        	polar: float
+        	     A value of the chosen linear response function corresponding to compute polariazabiltity in a specified cartesian diresction.
+        """
+
+        contract = self.ccwfn.contract
+
+        # Defining the l1 and l2
+        l1 = self.cclambda.l1
+        l2 = self.cclambda.l2
+        hbar = self.hbar
+        L = self.ccwfn.H.L
+        o = self.ccwfn.o
+        v = self.ccwfn.v
+
+        # Please refer to eqn 45 of [Crawford: 10.1002/wcms.1406].
+        # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
+        # <<A;B>> = <0|(1 + L^(0)) { [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)] + [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] } |0>
+
+        # <0|(1 + L^(0)) [\bar{A}^(0), X^(1)(B)] |0>
+        # <0| [\bar{A}^(0), X^(1)(B)] |0> + <0| L^(0) [\bar{A}^(0), X^(1)(B)] |0>
+        #           First term                      Second term
+
+        polar = 0.0
+        pertbar_A = self.pertbar[pertkey_a]
+        pertbar_B = self.pertbar[pertkey_b]
+        Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
+
+        # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
+        polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
+        # Second term (contributions from the L1 and X1)
+        temp = contract("ia, ic -> ac", l1, X1_B)
+        polar += contract("ac, ac -> ", temp, pertbar_A.Avv)
+        temp = contract("ia, ka -> ik", l1, X1_B)
+        polar -= contract("ik, ki -> ", temp, pertbar_A.Aoo)
+        # Second term (contributions from the L1 and X2)
+        temp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
+        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_B)
+        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_B)
+        # Second term (contributions from the L2 and X1)
+        temp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
+        polar += contract("ia, ia -> ", temp, X1_B)
+        temp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
+        polar -= 0.5 * contract("ak, ka -> ", temp, X1_B)
+        temp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
+        polar -= 0.5 * contract("bk, kb -> ", temp, X1_B)
+        # Second term (contributions from the L2 and X2)
+        temp = contract("ijab, kjab -> ik", l2, X2_B)
+        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_A.Aoo)
+        temp = contract("ijab, kiba-> jk", l2, X2_B)
+        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_A.Aoo)
+        temp = contract("ijab, ijac -> bc", l2, X2_B)
+        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_A.Avv)
+        temp = contract("ijab, ijcb -> ac", l2, X2_B)
+        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_A.Avv)
+
+        # <0|(1 + L^(0)) [\bar{B}^(1), X^(1)(B)] |0>
+        # <0| [\bar{B}^(1), X^(1)(B)] |0> + <0| L^(0) [\bar{B}^(0), X^(1)(B)] |0>
+        #           Third term                      Fourth term
+
+        # Third term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
+        polar += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
+        # Fourth term (contributions from the L1 and X1)
+        temp = contract("ia, ic -> ac", l1, X1_A)
+        polar += contract("ac, ac -> ", temp, pertbar_B.Avv)
+        temp = contract("ia, ka -> ik", l1, X1_A)
+        polar -= contract("ik, ki -> ", temp, pertbar_B.Aoo)
+        # Fourth term (contributions from the L1 and X2)
+        temp = contract("ia, jb -> ijab", l1, pertbar_B.Aov)
+        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_A)
+        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_A)
+        # Fourth term (contributions from the L2 and X1)
+        temp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
+        polar += contract("ia, ia -> ", temp, X1_A)
+        temp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
+        polar -= 0.5 * contract("ak, ka -> ", temp, X1_A)
+        temp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
+        polar -= 0.5 * contract("bk, kb -> ", temp, X1_A)
+        # Fourth term (contributions from the L2 and X2)
+        temp = contract("ijab, kjab -> ik", l2, X2_A)
+        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_B.Aoo)
+        temp = contract("ijab, kiba-> jk", l2, X2_A)
+        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_B.Aoo)
+        temp = contract("ijab, ijac -> bc", l2, X2_A)
+        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_B.Avv)
+        temp = contract("ijab, ijcb -> ac", l2, X2_A)
+        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
+
+        # <0|(1 + L^(0))[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
+        # <0|[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0> + <0| L^(0) [[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
+        #               Fifth term                          Sixth term
+
+        # Fifth term (Contribution from T1, there's no contribution from T2)
+        temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_B)
+        polar += 2.0 * contract("jb, jb -> ", temp, X1_A)
+
+        # Sixth term (contributions from the L1 and X1, also, no contributing term from L1 and X2)
+        temp = contract("ja, jb -> ab", hbar.Hov, X1_A)
+        temp = contract("ab, ia -> ib", temp, X1_B)
+        polar += -2.0 * contract("ib, ib -> ", l1, temp)
+        temp = contract("ijak, jc -> icak", L[o, o, v, o], X1_A)
+        temp = contract("icak, ia -> kc", temp, X1_B)
+        polar += - 2.0 * contract("kc, kc -> ", temp, l1)
+        temp = contract("cjab, jb -> ca", L[v, o, v, v], X1_A)
+        temp = contract("ca, ka -> kc", temp, X1_B)
+        polar += 2.0 * contract("kc, kc -> ", temp, l1)
+
+        # Sixth term (contributions from the L2 and X1)
+        temp = contract("jdka, jc -> cdka", hbar.Hovov, X1_A)
+        temp = contract("cdka, la -> klcd", temp, X1_B)
+        polar += -1.0 * contract("klcd, klcd -> ", temp, l2)
+        temp = contract("cjka, jd -> cdka", hbar.Hovov.swapaxes(0, 1), X1_A)
+        temp = contract("cdka, la -> cdkl", temp, X1_B)
+        polar += -1.0 * contract("cdkl, klcd -> ", temp, l2)
+
+
+        # Sixth term (contributions from the L2 and X2)
+        # I don't know if I should use the L term or the H_bar term for Hoovv
+        # There is a permutationally
+        temp = contract("klab, mnab -> klmn", L[o, o, v, v], X2_B)
+        temp = contract("klmn, klfe -> mnef", temp, X2_A)
+        polar += contract("mnef, mnef -> ", temp, l2)
+        temp = contract("klab, nkaf -> lbfn", L[o, o, v, v], X2_B)
+        temp = contract("lbfn, mlbe -> mnef", temp, X2_A)
+        polar -= contract("mnef, mnef -> ", temp, l2)
+        temp = contract("klab, nlae -> kbne", L[o, o, v, v], X2_B)
+        temp = contract("kbne, mkbf -> mnef", temp, X2_A)
+        polar += contract("mnef, mnef -> ", temp, l2)
+
+        # Also, we can consider terms like L1, X1, X2 == L1, X2, X1 or L2, X1, X2
+        # L1, X1, X2
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        temp = contract("jc, jldc -> ld", temp, X2_A)
+        polar += -1.0 * contract("ld, ld -> ", temp, l1)
+        temp = contract("ijac, ijda -> dc", L[o, o, v, v], X2_B)
+        temp = contract("dc, lc -> ld", temp, X1_A)
+        polar += contract("ld, ld -> ", temp, l1)
+        temp = contract("ijac, ilac-> jl", L[o, o, v, v], X2_B)
+        temp = contract("jl, jd -> ld", temp, X1_A)
+        polar += -1.0 * contract("ld, ld -> ", temp, l1)
+
+        # L2, X1, X2
+        temp = contract("dkab, ma -> dkmb", hbar.Hvovv, X1_B)
+        temp = contract("dkmb, lkbe -> lmde", temp, X2_A)
+        polar += 0.5 * contract("lmde, lmde -> ", temp, l2)
+        temp = contract("dkab, mkae -> dmbe", hbar.Hvovv, X2_B)
+        temp = contract("dmbe, lb -> dmel", temp, X1_A)
+        polar -= 0.5 * contract("dmel, lmde -> ", temp, l2)
+        temp = contract("dkab, mlab -> mldk", hbar.Hvovv, X2_B)
+        temp = contract("mldk, ke -> dmel", temp, X1_A)
+        polar += 0.5 * contract("dmel, lmde -> ", temp, l2)
+
+        temp = contract("kjla, klde -> ajde", hbar.Hooov, X2_B)
+        temp = contract("ajde, ia -> ijde", temp, X1_A)
+        polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
+        temp = contract("klja, ikae -> ijle", hbar.Hooov, X2_B)
+        temp = contract("ijle, ld -> ijde", temp, X1_A)
+        polar += 0.5 * contract("ijde, ijde -> ", temp, l2)
+        temp = contract("klja, ilad -> ijkd", hbar.Hooov, X2_B)
+        temp = contract("ijkd, ke -> ijde", temp, X1_A)
+        polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
+
+        return -1.0 * polar
+
     def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         """
 	Calculate the CC linear response function for polarizability at field-frequency omega(w1).
@@ -328,14 +504,14 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         polar2 = 0
         pertbar_A = self.pertbar[pertkey_a]
         Avvoo = pertbar_A.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
-        # <0|Y1(B) * A_bar|0>
-        polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
-        # <0|Y2(B) * A_bar|0>
-        polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
-        polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
+        # # <0|Y1(B) * A_bar|0>
+        # polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
+        # # <0|Y2(B) * A_bar|0>
+        # polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
+        # polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
         # <0|[A_bar, X(B)]|0>
         polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
-        # <0|L1(0) [A_bar, X2(B)]|0>
+        # <0|L1(0) [A_bar, X1(B)]|0>
         tmp = contract("ia, ic -> ac", l1, X1_B)
         polar2 += contract("ac, ac -> ", tmp, pertbar_A.Avv)
         tmp = contract("ia, ka -> ik", l1, X1_B)
@@ -351,7 +527,7 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         polar2 -= 0.5 * contract("ak, ka -> ", tmp, X1_B)
         tmp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
         polar2 -= 0.5 * contract("bk, kb -> ", tmp, X1_B)
-        # <0|L2(0)[A_bar, X1(B)]|0>
+        # <0|L2(0)[A_bar, X2(B)]|0>
         tmp = contract("ijab, kjab -> ik", l2, X2_B)
         polar2 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_A.Aoo)
         tmp = contract("ijab, kiba-> jk", l2, X2_B)

From c384f965b9695e07a9fba293ea2d4fbbef529caf Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Wed, 19 Jun 2024 11:13:31 -0400
Subject: [PATCH 2/7] Implementing LR symmetric

---
 pycc/ccresponse.py            | 336 +++++++++++++++++++++-------------
 pycc/tests/test_036_lr.py     | 179 ++++++++++++------
 pycc/tests/test_036_lr_sym.py | 141 ++++++++++++++
 3 files changed, 467 insertions(+), 189 deletions(-)
 create mode 100644 pycc/tests/test_036_lr_sym.py

diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py
index 0f34b70..dd2a625 100644
--- a/pycc/ccresponse.py
+++ b/pycc/ccresponse.py
@@ -9,6 +9,7 @@
 import time
 from .utils import helper_diis
 from .cclambda import cclambda
+from .hamiltonian import Hamiltonian
 
 class ccresponse(object):
     """
@@ -314,14 +315,14 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         """
 
         contract = self.ccwfn.contract
+        o = self.ccwfn.o
+        v = self.ccwfn.v
 
         # Defining the l1 and l2
         l1 = self.cclambda.l1
         l2 = self.cclambda.l2
         hbar = self.hbar
         L = self.ccwfn.H.L
-        o = self.ccwfn.o
-        v = self.ccwfn.v
 
         # Please refer to eqn 45 of [Crawford: 10.1002/wcms.1406].
         # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
@@ -336,141 +337,216 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         pertbar_B = self.pertbar[pertkey_b]
         Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
 
-        # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
-        polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
-        # Second term (contributions from the L1 and X1)
-        temp = contract("ia, ic -> ac", l1, X1_B)
-        polar += contract("ac, ac -> ", temp, pertbar_A.Avv)
-        temp = contract("ia, ka -> ik", l1, X1_B)
-        polar -= contract("ik, ki -> ", temp, pertbar_A.Aoo)
-        # Second term (contributions from the L1 and X2)
-        temp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
-        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_B)
-        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_B)
-        # Second term (contributions from the L2 and X1)
-        temp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
-        polar += contract("ia, ia -> ", temp, X1_B)
-        temp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
-        polar -= 0.5 * contract("ak, ka -> ", temp, X1_B)
-        temp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
-        polar -= 0.5 * contract("bk, kb -> ", temp, X1_B)
-        # Second term (contributions from the L2 and X2)
-        temp = contract("ijab, kjab -> ik", l2, X2_B)
-        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_A.Aoo)
-        temp = contract("ijab, kiba-> jk", l2, X2_B)
-        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_A.Aoo)
-        temp = contract("ijab, ijac -> bc", l2, X2_B)
-        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_A.Avv)
-        temp = contract("ijab, ijcb -> ac", l2, X2_B)
-        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_A.Avv)
-
-        # <0|(1 + L^(0)) [\bar{B}^(1), X^(1)(B)] |0>
-        # <0| [\bar{B}^(1), X^(1)(B)] |0> + <0| L^(0) [\bar{B}^(0), X^(1)(B)] |0>
-        #           Third term                      Fourth term
-
-        # Third term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
-        polar += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
-        # Fourth term (contributions from the L1 and X1)
-        temp = contract("ia, ic -> ac", l1, X1_A)
-        polar += contract("ac, ac -> ", temp, pertbar_B.Avv)
-        temp = contract("ia, ka -> ik", l1, X1_A)
-        polar -= contract("ik, ki -> ", temp, pertbar_B.Aoo)
-        # Fourth term (contributions from the L1 and X2)
-        temp = contract("ia, jb -> ijab", l1, pertbar_B.Aov)
-        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_A)
-        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_A)
-        # Fourth term (contributions from the L2 and X1)
-        temp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
-        polar += contract("ia, ia -> ", temp, X1_A)
-        temp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
-        polar -= 0.5 * contract("ak, ka -> ", temp, X1_A)
-        temp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
-        polar -= 0.5 * contract("bk, kb -> ", temp, X1_A)
-        # Fourth term (contributions from the L2 and X2)
-        temp = contract("ijab, kjab -> ik", l2, X2_A)
-        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_B.Aoo)
-        temp = contract("ijab, kiba-> jk", l2, X2_A)
-        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_B.Aoo)
-        temp = contract("ijab, ijac -> bc", l2, X2_A)
-        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_B.Avv)
-        temp = contract("ijab, ijcb -> ac", l2, X2_A)
-        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
+        # # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
+        # polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
+        # # Second term (contributions from the L1 and X1)
+        # temp = contract("ia, ic -> ac", l1, X1_B)
+        # polar += contract("ac, ac -> ", temp, pertbar_A.Avv)
+        # temp = contract("ia, ka -> ik", l1, X1_B)
+        # polar -= contract("ik, ki -> ", temp, pertbar_A.Aoo)
+        # # Second term (contributions from the L1 and X2)
+        # temp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
+        # polar += 2.0 * contract("ijab, ijab -> ", temp, X2_B)
+        # polar += -1.0 * contract("ijab, ijba -> ", temp, X2_B)
+        # # Second term (contributions from the L2 and X1)
+        # temp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
+        # polar += contract("ia, ia -> ", temp, X1_B)
+        # temp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
+        # polar -= 0.5 * contract("ak, ka -> ", temp, X1_B)
+        # temp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
+        # polar -= 0.5 * contract("bk, kb -> ", temp, X1_B)
+        # # Second term (contributions from the L2 and X2)
+        # temp = contract("ijab, kjab -> ik", l2, X2_B)
+        # polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_A.Aoo)
+        # temp = contract("ijab, kiba-> jk", l2, X2_B)
+        # polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_A.Aoo)
+        # temp = contract("ijab, ijac -> bc", l2, X2_B)
+        # polar += 0.5 * contract("bc, bc -> ", temp, pertbar_A.Avv)
+        # temp = contract("ijab, ijcb -> ac", l2, X2_B)
+        # polar += 0.5 * contract("ac, ac -> ", temp, pertbar_A.Avv)
+        #
+        # # <0|(1 + L^(0)) [\bar{B}^(1), X^(1)(B)] |0>
+        # # <0| [\bar{B}^(1), X^(1)(B)] |0> + <0| L^(0) [\bar{B}^(0), X^(1)(B)] |0>
+        # #           Third term                      Fourth term
+        #
+        # # Third term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
+        # polar += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
+        # # Fourth term (contributions from the L1 and X1)
+        # temp = contract("ia, ic -> ac", l1, X1_A)
+        # polar += contract("ac, ac -> ", temp, pertbar_B.Avv)
+        # temp = contract("ia, ka -> ik", l1, X1_A)
+        # polar -= contract("ik, ki -> ", temp, pertbar_B.Aoo)
+        # # Fourth term (contributions from the L1 and X2)
+        # temp = contract("ia, jb -> ijab", l1, pertbar_B.Aov)
+        # polar += 2.0 * contract("ijab, ijab -> ", temp, X2_A)
+        # polar += -1.0 * contract("ijab, ijba -> ", temp, X2_A)
+        # # Fourth term (contributions from the L2 and X1)
+        # temp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
+        # polar += contract("ia, ia -> ", temp, X1_A)
+        # temp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
+        # polar -= 0.5 * contract("ak, ka -> ", temp, X1_A)
+        # temp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
+        # polar -= 0.5 * contract("bk, kb -> ", temp, X1_A)
+        # # Fourth term (contributions from the L2 and X2)
+        # temp = contract("ijab, kjab -> ik", l2, X2_A)
+        # polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_B.Aoo)
+        # temp = contract("ijab, kiba-> jk", l2, X2_A)
+        # polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_B.Aoo)
+        # temp = contract("ijab, ijac -> bc", l2, X2_A)
+        # polar += 0.5 * contract("bc, bc -> ", temp, pertbar_B.Avv)
+        # temp = contract("ijab, ijcb -> ac", l2, X2_A)
+        # polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
 
         # <0|(1 + L^(0))[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
         # <0|[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0> + <0| L^(0) [[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
         #               Fifth term                          Sixth term
 
-        # Fifth term (Contribution from T1, there's no contribution from T2)
-        temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_B)
-        polar += 2.0 * contract("jb, jb -> ", temp, X1_A)
+        # Expanding the permutational operator and implementing it explicitly
+        temp = contract("ijab, ia -> jb", L[o, o, v, v].copy(), X1_A)
+        polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
+
+        temp = contract("je, ja -> ea", X1_A, hbar.Hov)
+        temp = contract("ea, ma -> me", temp, X1_B)
+        polar -= contract("me, me -> ", temp, l1)
+        temp = contract("ma, me -> ea", X1_A, hbar.Hov)
+        temp = contract("ea, je -> ja", temp, X1_B)
+        polar -= contract("ja, me -> ", temp, l1)
+
+        temp = contract("ijam, je -> ieam", hbar.Hooov.swapaxes(2, 3), X1_A)
+        temp = contract("ieam, ia -> me", temp, X1_B)
+        polar -= contract("me, me -> ", temp, l1)
+        temp = contract("ijem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
+        temp = contract("ajem, je -> ma", temp, X1_B)
+        polar -= contract("ma, ma -> ", temp, l1)
+
+        temp = contract("ejab, jb -> ea", hbar.Hvovv, X1_A)
+        temp = contract("ea, ma -> me", temp, X1_B)
+        polar += contract("me, me -> ", temp, l1)
+        temp = contract("emab, jb -> emaj", hbar.Hvovv, X1_A)
+        temp = contract("emaj, ma -> je", temp, X1_B)
+        polar += contract("je, je -> ", temp, l1)
+
+        temp = contract("jfma, je -> efma", hbar.Hovov, X1_A)
+        temp = contract("efma, na -> mnef", temp, X1_B)
+        polar -= contract("mnef, mnef -> ", temp, l2)
+        temp = contract("nfma, na -> mf", hbar.Hovov, X1_A)
+        temp = contract("mf, je -> mjef", temp, X1_B)
+        polar -= contract("mjef, mjef -> ", temp, l2)
+
+        temp = contract("ejma, jf -> efma", hbar.Hovov.swapaxes(0, 1), X1_A)
+        temp = contract("efma, na -> mnef", temp, X1_B)
+        polar -= contract("mnef, mnef -> ", temp, l2)
+        temp = contract("enmf, na -> amfe", hbar.Hovov.swapaxes(0, 1), X1_A)
+        temp = contract("amfe, jf -> mjae", temp, X1_B)
+        polar -= contract("mjae, mjae -> ", temp, l2)
+
+        temp = contract("abef, nb -> anef", hbar.Hvvvv, X1_A)
+        temp = contract("anef, ma -> mnef", temp, X1_B)
+        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
+        temp = contract("abef, ma -> amef", hbar.Hvvvv, X1_A)
+        temp = contract("amef, nb -> mnef", temp, X1_B)
+        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
+
+        temp = contract("mnij, jf -> mnif", hbar.Hoooo, X1_A)
+        temp = contract("mnif, ie -> mnef", temp, X1_B)
+        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
+        temp = contract("mnij, ie -> mnje", hbar.Hoooo, X1_A)
+        temp = contract("mnje, jf -> mnef", temp, X1_B)
+        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
+
+
 
-        # Sixth term (contributions from the L1 and X1, also, no contributing term from L1 and X2)
-        temp = contract("ja, jb -> ab", hbar.Hov, X1_A)
-        temp = contract("ab, ia -> ib", temp, X1_B)
-        polar += -2.0 * contract("ib, ib -> ", l1, temp)
-        temp = contract("ijak, jc -> icak", L[o, o, v, o], X1_A)
-        temp = contract("icak, ia -> kc", temp, X1_B)
-        polar += - 2.0 * contract("kc, kc -> ", temp, l1)
-        temp = contract("cjab, jb -> ca", L[v, o, v, v], X1_A)
-        temp = contract("ca, ka -> kc", temp, X1_B)
-        polar += 2.0 * contract("kc, kc -> ", temp, l1)
-
-        # Sixth term (contributions from the L2 and X1)
-        temp = contract("jdka, jc -> cdka", hbar.Hovov, X1_A)
-        temp = contract("cdka, la -> klcd", temp, X1_B)
-        polar += -1.0 * contract("klcd, klcd -> ", temp, l2)
-        temp = contract("cjka, jd -> cdka", hbar.Hovov.swapaxes(0, 1), X1_A)
-        temp = contract("cdka, la -> cdkl", temp, X1_B)
-        polar += -1.0 * contract("cdkl, klcd -> ", temp, l2)
 
 
+
+
+
+
+
+
+
+
+        # Begin HX_1Y_1
+        # Fifth term (Contribution from T1, there's no contribution from T2)
+        # temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)  # L[o, o, v, v] or hbar.Hoovv
+        # polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
+
+        # Sixth term (contributions from the L1 and X1, also, no contributing term from L1 and X2)
+        # temp = contract("ja, jb -> ab", hbar.Hov, X1_A)
+        # temp = contract("ab, ia -> ib", temp, X1_B)
+        # polar -= 2.0 * contract("ib, ib -> ", l1, temp)
+        # temp = contract("ijka, jc -> icak", hbar.Hooov, X1_A)    # L[o, o, v, o] or hbar.Hoovo
+        # temp = contract("icak, ia -> kc", temp, X1_B)
+        # polar += 2.0 * contract("kc, kc -> ", temp, l1)
+        # temp = contract("cjab, jb -> ca", hbar.Hvovv, X1_A)      ## L[v, o, v, v] or hbar.Hvovv
+        # temp = contract("ca, ka -> kc", temp, X1_B)
+        # polar += 2.0 * contract("kc, kc -> ", temp, l1)
+        #
+        # # Sixth term (contributions from the L2 and X1)
+        # temp = contract("je, jfma -> amef", X1_A, hbar.Hovov)
+        # temp = contract("amef, na -> mnef", temp, X1_B)
+        # polar -= contract("mnef, mnef -> ", temp, l2)
+        # temp = contract("je, fjma -> amef", X1_A, hbar.Hovov.swapaxes(0, 1))
+        # temp = contract("amef, na -> mnef", temp, X1_B)
+        # polar -= contract("mnef, mnef -> ", temp, l2)
+        # temp = contract("nb, abef -> anef", X1_A, hbar.Hvvvv)
+        # temp = contract("anef, ma -> mnef", temp, X1_B)
+        # polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
+        # temp = contract("jf, mnij -> mnif", X1_A, hbar.Hoooo)
+        # temp = contract("mnif, ie -> mnef", temp, X1_B)
+        # polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
+        # End HX_1Y_1
+
+
+        # Begin HX_{2}Y_{2}
+        # There is no contribution from L1 and X2
         # Sixth term (contributions from the L2 and X2)
         # I don't know if I should use the L term or the H_bar term for Hoovv
         # There is a permutationally
-        temp = contract("klab, mnab -> klmn", L[o, o, v, v], X2_B)
-        temp = contract("klmn, klfe -> mnef", temp, X2_A)
-        polar += contract("mnef, mnef -> ", temp, l2)
-        temp = contract("klab, nkaf -> lbfn", L[o, o, v, v], X2_B)
-        temp = contract("lbfn, mlbe -> mnef", temp, X2_A)
-        polar -= contract("mnef, mnef -> ", temp, l2)
-        temp = contract("klab, nlae -> kbne", L[o, o, v, v], X2_B)
-        temp = contract("kbne, mkbf -> mnef", temp, X2_A)
-        polar += contract("mnef, mnef -> ", temp, l2)
-
-        # Also, we can consider terms like L1, X1, X2 == L1, X2, X1 or L2, X1, X2
-        # L1, X1, X2
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
-        temp = contract("jc, jldc -> ld", temp, X2_A)
-        polar += -1.0 * contract("ld, ld -> ", temp, l1)
-        temp = contract("ijac, ijda -> dc", L[o, o, v, v], X2_B)
-        temp = contract("dc, lc -> ld", temp, X1_A)
-        polar += contract("ld, ld -> ", temp, l1)
-        temp = contract("ijac, ilac-> jl", L[o, o, v, v], X2_B)
-        temp = contract("jl, jd -> ld", temp, X1_A)
-        polar += -1.0 * contract("ld, ld -> ", temp, l1)
-
-        # L2, X1, X2
-        temp = contract("dkab, ma -> dkmb", hbar.Hvovv, X1_B)
-        temp = contract("dkmb, lkbe -> lmde", temp, X2_A)
-        polar += 0.5 * contract("lmde, lmde -> ", temp, l2)
-        temp = contract("dkab, mkae -> dmbe", hbar.Hvovv, X2_B)
-        temp = contract("dmbe, lb -> dmel", temp, X1_A)
-        polar -= 0.5 * contract("dmel, lmde -> ", temp, l2)
-        temp = contract("dkab, mlab -> mldk", hbar.Hvovv, X2_B)
-        temp = contract("mldk, ke -> dmel", temp, X1_A)
-        polar += 0.5 * contract("dmel, lmde -> ", temp, l2)
-
-        temp = contract("kjla, klde -> ajde", hbar.Hooov, X2_B)
-        temp = contract("ajde, ia -> ijde", temp, X1_A)
-        polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
-        temp = contract("klja, ikae -> ijle", hbar.Hooov, X2_B)
-        temp = contract("ijle, ld -> ijde", temp, X1_A)
-        polar += 0.5 * contract("ijde, ijde -> ", temp, l2)
-        temp = contract("klja, ilad -> ijkd", hbar.Hooov, X2_B)
-        temp = contract("ijkd, ke -> ijde", temp, X1_A)
-        polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
-
-        return -1.0 * polar
+        temp = contract("klab, mnab -> klmn", L[o, o, v, v], X2_A)
+        temp = contract("klmn, klfe -> mnef", temp, X2_B)
+        polar += 1/16 * contract("mnef, mnef -> ", temp, l2)
+        # temp = contract("klab, nkaf -> lbfn", L[o, o, v, v], X2_A)
+        # temp = contract("lbfn, mlbe -> mnef", temp, X2_B)
+        # polar -= contract("mnef, mnef -> ", temp, l2)
+        # temp = contract("klab, nlae -> kbne", L[o, o, v, v], X2_A)
+        # temp = contract("kbne, mkbf -> mnef", temp, X2_B)
+        # polar += contract("mnef, mnef -> ", temp, l2)
+        #
+        # # Also, we can consider terms like L1, X1, X2 == L1, X2, X1 or L2, X1, X2
+        # # L1, X1, X2
+        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        # temp = contract("jc, jldc -> ld", temp, X2_A)
+        # polar += -1.0 * contract("ld, ld -> ", temp, l1)
+        # temp = contract("ijac, ijda -> dc", L[o, o, v, v], X2_B)
+        # temp = contract("dc, lc -> ld", temp, X1_A)
+        # polar += contract("ld, ld -> ", temp, l1)
+        # temp = contract("ijac, ilac-> jl", L[o, o, v, v], X2_B)
+        # temp = contract("jl, jd -> ld", temp, X1_A)
+        # polar += -1.0 * contract("ld, ld -> ", temp, l1)
+        #
+        # # L2, X1, X2
+        # temp = contract("dkab, ma -> dkmb", hbar.Hvovv, X1_B)
+        # temp = contract("dkmb, lkbe -> lmde", temp, X2_A)
+        # polar += 0.5 * contract("lmde, lmde -> ", temp, l2)
+        # temp = contract("dkab, mkae -> dmbe", hbar.Hvovv, X2_B)
+        # temp = contract("dmbe, lb -> dmel", temp, X1_A)
+        # polar -= 0.5 * contract("dmel, lmde -> ", temp, l2)
+        # temp = contract("dkab, mlab -> mldk", hbar.Hvovv, X2_B)
+        # temp = contract("mldk, ke -> dmel", temp, X1_A)
+        # polar += 0.5 * contract("dmel, lmde -> ", temp, l2)
+        #
+        # temp = contract("kjla, klde -> ajde", hbar.Hooov, X2_B)
+        # temp = contract("ajde, ia -> ijde", temp, X1_A)
+        # polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
+        # temp = contract("klja, ikae -> ijle", hbar.Hooov, X2_B)
+        # temp = contract("ijle, ld -> ijde", temp, X1_A)
+        # polar += 0.5 * contract("ijde, ijde -> ", temp, l2)
+        # temp = contract("klja, ilad -> ijkd", hbar.Hooov, X2_B)
+        # temp = contract("ijkd, ke -> ijde", temp, X1_A)
+        # polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
+
+        return polar#-1.0 * polar
 
     def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         """
@@ -505,10 +581,10 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         pertbar_A = self.pertbar[pertkey_a]
         Avvoo = pertbar_A.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
         # # <0|Y1(B) * A_bar|0>
-        # polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
-        # # <0|Y2(B) * A_bar|0>
-        # polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
-        # polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
+        polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
+        # <0|Y2(B) * A_bar|0>
+        polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
+        polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
         # <0|[A_bar, X(B)]|0>
         polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
         # <0|L1(0) [A_bar, X1(B)]|0>
diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py
index 3b2bf78..a555240 100644
--- a/pycc/tests/test_036_lr.py
+++ b/pycc/tests/test_036_lr.py
@@ -5,74 +5,135 @@
 # Import package, test suite, and other packages as needed
 import psi4
 import numpy as np
+import sys
+sys.path.append("/Users/jattakumi/pycc/pycc/")
 import pytest
 import pycc
-from ..data.molecules import *
+from data.molecules import *
 
 def test_linresp():
+    h2o2 = """
+    O -0.182400 -0.692195 -0.031109
+    O 0.182400 0.692195 -0.031109
+    H 0.533952 -1.077444 0.493728
+    H -0.533952 1.077444 0.493728
+    symmetry c1
+    """
+
+    # H2_4 chain
+    h2_4 = """
+    H 0.000000 0.000000 0.000000
+    H 0.750000 0.000000 0.000000
+    H 0.000000 1.500000 0.000000
+    H 0.375000 1.500000 -0.649520
+    H 0.000000 3.000000 0.000000
+    H -0.375000 3.000000 -0.649520
+    H 0.000000 4.500000 -0.000000
+    H -0.750000 4.500000 -0.000000
+    symmetry c1
+    """
+
+    # H2_7 chain
+    h2_7 = """
+    H 0.000000 0.000000 0.000000
+    H 0.750000 0.000000 0.000000
+    H 0.000000 1.500000 0.000000
+    H 0.375000 1.500000 -0.649520
+    H 0.000000 3.000000 0.000000
+    H -0.375000 3.000000 -0.649520
+    H 0.000000 4.500000 -0.000000
+    H -0.750000 4.500000 -0.000000
+    H 0.000000 6.000000 -0.000000
+    H -0.375000 6.000000 0.649520
+    H 0.000000 7.500000 -0.000000
+    H 0.375000 7.500000 -0.649520
+    H 0.000000 9.000000 -0.000000
+    H 0.750000 9.000000 0.000000
+    symmetry c1
+    """
+
+    threshold = [1e-03, 1e-04, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10, 1e-11, 1e-12, 0]
+    lmo = 'PNO++'
+    geom = "h2_4"
+    psi4.core.clean_options()
     psi4.set_memory('2 GiB')
-    psi4.core.set_output_file('output.dat', False)
+    psi4.set_output_file('output.dat', False)
     psi4.set_options({'basis': 'aug-cc-pvdz',
                       'scf_type': 'pk',
+                      'mp2_type': 'conv',
+                      'freeze_core': 'true',
                       'e_convergence': 1e-12,
                       'd_convergence': 1e-12,
-                      'r_convergence': 1e-12
+                      'r_convergence': 1e-12,
     })
-    mol = psi4.geometry(moldict["H2O"])
+
+    mol = psi4.geometry(h2_4)
     rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
 
-    e_conv = 1e-12
-    r_conv = 1e-12
-
-    cc = pycc.ccwfn(rhf_wfn)
-    ecc = cc.solve_cc(e_conv, r_conv)
-    hbar = pycc.cchbar(cc)
-    cclambda = pycc.cclambda(cc, hbar)
-    lecc = cclambda.solve_lambda(e_conv, r_conv)
-    density = pycc.ccdensity(cc, cclambda)
-
-    resp = pycc.ccresponse(density)
-
-    omega1 = 0.0656
-
-    # Creating dictionaries
-    # X_1 = X(-omega); X_2 = X(omega)
-    # Y_1 = Y(-omega); Y_2 = Y(omega)
-    X_1 = {}
-    X_2 = {}
-    Y_1 = {}
-    Y_2 = {}
-
-    for axis in range(0, 3):
-        string = "MU_" + resp.cart[axis]
-
-        A = resp.pertbar[string]
-
-        X_2[string] = resp.solve_right(A, omega1)
-        Y_2[string] = resp.solve_left(A, omega1)
-        X_1[string] = resp.solve_right(A, -omega1)
-        Y_1[string] = resp.solve_left(A, -omega1)
-
-    # Grabbing X, Y and declaring the matrix space for LR
-    polar_AB = np.zeros((3,3))
-
-    for a in range(0, 3):
-        string_a = "MU_" + resp.cart[a]
-        for b in range(0,3):
-            string_b = "MU_" + resp.cart[b]
-            Y1_B, Y2_B, _ = Y_2[string_b]
-            X1_B, X2_B, _ = X_2[string_b]
-            polar_AB[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
-
-    polar_AB_avg = np.average([polar_AB[0,0], polar_AB[1,1], polar_AB[2,2]])
-
-    #validating from psi4
-    polar_XX = 9.92992070420665
-    polar_YY = 13.443740151331559
-    polar_ZZ = 11.342765745046526
-    polar_avg = 11.572142200333
-
-    assert(abs(polar_AB[0,0] - polar_XX) < 1e-8)
-    assert(abs(polar_AB[1,1] - polar_YY) < 1e-8)
-    assert(abs(polar_AB[2,2] - polar_ZZ) < 1e-8)
-    assert(abs(polar_AB_avg - polar_avg) < 1e-8)
+    for t in threshold:
+        e_conv = 1e-12
+        r_conv = 1e-12
+
+        cc = pycc.ccwfn(rhf_wfn, local_mos = 'BOYS', local= lmo, local_cutoff = t, filter=True)
+        ecc = cc.solve_cc(e_conv, r_conv)
+        hbar = pycc.cchbar(cc)
+        cclambda = pycc.cclambda(cc, hbar)
+        lecc = cclambda.solve_lambda(e_conv, r_conv)
+        density = pycc.ccdensity(cc, cclambda)
+
+        resp = pycc.ccresponse(density)
+
+        omega1 = 0.0656
+
+        # Creating dictionaries
+        # X_1 = X(-omega); X_2 = X(omega)
+        # Y_1 = Y(-omega); Y_2 = Y(omega)
+        X_1 = {}
+        X_2 = {}
+        Y_1 = {}
+        Y_2 = {}
+
+        for axis in range(0, 3):
+            string = "MU_" + resp.cart[axis]
+
+            A = resp.pertbar[string]
+
+            X_2[string] = resp.solve_right(A, omega1)
+            Y_2[string] = resp.solve_left(A, omega1)
+            X_1[string] = resp.solve_right(A, -omega1)
+            Y_1[string] = resp.solve_left(A, -omega1)
+
+        # Grabbing X, Y and declaring the matrix space for LR
+        polar_AB = np.zeros((3,3))
+
+        for a in range(0, 3):
+            string_a = "MU_" + resp.cart[a]
+            for b in range(0,3):
+                string_b = "MU_" + resp.cart[b]
+                Y1_B, Y2_B, _ = Y_2[string_b]
+                X1_B, X2_B, _ = X_2[string_b]
+                polar_AB[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
+
+        print("Dynamic Polarizability Tensor @ w=0.0656 a.u.:")
+        print(polar_AB)
+        print("Average Dynamic Polarizability:")
+        polar_AB_avg = np.average([polar_AB[0,0], polar_AB[1,1], polar_AB[2,2]])
+        print(polar_AB_avg)
+
+        B_avg = str(geom) + "_Dynpolar_" + str(lmo) + ".txt"
+
+        with open(B_avg, 'a') as f1:
+            f1.write(str(polar_AB_avg) + ' ' + str(t) + '\n')
+
+        del cc, hbar, cclambda, density, resp
+
+    # #validating from psi4
+    # polar_XX = 9.92992070420665
+    # polar_YY = 13.443740151331559
+    # polar_ZZ = 11.342765745046526
+    # polar_avg = 11.572142200333
+    #
+    # assert(abs(polar_AB[0,0] - polar_XX) < 1e-8)
+    # assert(abs(polar_AB[1,1] - polar_YY) < 1e-8)
+    # assert(abs(polar_AB[2,2] - polar_ZZ) < 1e-8)
+    # assert(abs(polar_AB_avg - polar_avg) < 1e-8)
diff --git a/pycc/tests/test_036_lr_sym.py b/pycc/tests/test_036_lr_sym.py
new file mode 100644
index 0000000..c0354d3
--- /dev/null
+++ b/pycc/tests/test_036_lr_sym.py
@@ -0,0 +1,141 @@
+"""
+Test CCSD linear response functions.
+"""
+import sys
+sys.path.append("/Users/jattakumi/pycc/")
+import numpy as np
+# Import package, test suite, and other packages as needed
+import psi4
+from opt_einsum import contract
+
+import pycc
+#from ccwfn import ccwfn
+#from cchbar import cchbar
+#from cclambda import cclambda
+#from ccdensity import ccdensity
+#from ccresponse import ccresponse
+
+geom = """
+O
+H 1 1.8084679
+H 1 1.8084679 2 104.5
+units bohr
+symmetry c1 
+no_reorient
+"""
+
+hf = """
+F  0.000000000000   0.000000000000  0.000000000000
+H  0.000000000000   0.000000000000  -1.732800000000
+units bohr
+no_reorient
+symmetry c1
+"""
+
+hof = """
+          O          -0.947809457408    -0.132934425181     0.000000000000
+          H          -1.513924046286     1.610489987673     0.000000000000
+          F           0.878279174340     0.026485523618     0.000000000000
+unit bohr
+no_reorient
+symmetry c1
+"""
+psi4.core.clean_options()
+psi4.set_memory('2 GiB')
+psi4.core.set_output_file('output.dat', False)
+psi4.set_options({'basis': 'cc-pVDZ',
+                  'scf_type': 'pk',
+                  'mp2_type': 'conv',
+                  'freeze_core': 'true',
+                  'e_convergence': 1e-12,
+                  'd_convergence': 1e-12,
+                  'r_convergence': 1e-12})
+mol = psi4.geometry(hof)
+rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+e_conv = 1e-12
+r_conv = 1e-12
+
+cc = pycc.ccwfn(rhf_wfn)
+ecc = cc.solve_cc(e_conv, r_conv)
+hbar = pycc.cchbar(cc)
+cclambda = pycc.cclambda(cc, hbar)
+lecc = cclambda.solve_lambda(e_conv, r_conv)
+density = pycc.ccdensity(cc, cclambda)
+
+resp = pycc.ccresponse(density)
+
+omega1 = 0.0
+# omega1 = 0.0656
+
+
+#A = resp.pertbar.Aoo
+
+#resp.linresp(omega1)
+# Creating dictionaries
+# X_1 = X(-omega); X_2 = X(omega)
+# Y_1 = Y(-omega); Y_2 = Y(omega)
+# X_neg = X1(-omega) , X2(-omega)
+# X_pos, Y_neg, Y_pos
+# X_1 = {}
+# X_2 = {}
+# Y_1 = {}
+# Y_2 = {}
+X_A = {}
+X_B = {}
+
+for axis in range(0, 3):
+    string = "MU_" + resp.cart[axis]
+
+    A = resp.pertbar[string]
+    # B = resp.pertbar[string]
+    #print("Aoo",A)
+
+    # A -> -omega
+    # B -> +omega
+    X_A[string] = resp.solve_right(A, omega1)
+    X_B[string] = resp.solve_right(A, -omega1)
+
+    # X_2[string] = resp.solve_right(A, omega1)
+    # Y_2[string] = resp.solve_left(A, omega1)
+    # X_1[string] = resp.solve_right(A, -omega1)
+    # Y_1[string] = resp.solve_left(A, -omega1)
+
+#resp.polar(omega1)
+# Grabbing X, Y and declaring the matrix space for LR
+# polar_AB_pos = np.zeros((3,3))
+polar_AB_neg = np.zeros((3,3))
+# polar_AB_aveg = np.zeros((3,3))
+
+
+for a in range(0, 3):
+    string_a = "MU_" + resp.cart[a]
+    X1_A, X2_A, _ = X_A[string_a]
+    for b in range(0, 3):
+        string_b = "MU_" + resp.cart[b]
+        # Y1_B, Y2_B, _ = Y_2[string_b]
+        # X1_B, X2_B, _ = X_2[string_b]
+        X_1B, X_2B, _ = X_B[string_b]
+        # polar_AB_pos[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
+        polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
+        # polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
+
+# print(polar_AB_pos[0,0])
+# print(polar_AB_pos[0,1])
+# print(polar_AB_pos[0,2])
+# print(polar_AB_pos[1,0])
+# print(polar_AB_pos[1,1])
+# print(polar_AB_pos[1,2])
+# print(polar_AB_pos[2,0])
+# print(polar_AB_pos[2,1])
+# print(polar_AB_pos[2,2])
+# print("Comparing asymmetric versus symmetric")
+print(polar_AB_neg[0,0])
+print(polar_AB_neg[0,1])
+print(polar_AB_neg[0,2])
+print(polar_AB_neg[1,0])
+print(polar_AB_neg[1,1])
+print(polar_AB_neg[1,2])
+print(polar_AB_neg[2,0])
+print(polar_AB_neg[2,1])
+print(polar_AB_neg[2,2])

From fa2f26cadf8fca46011de3f6cb417b744157880b Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Mon, 2 Sep 2024 15:22:47 -0400
Subject: [PATCH 3/7] Still working on the LR symmetric

---
 pycc/ccresponse.py            | 368 ++++++++++++++++++++++++----------
 pycc/tests/test_036_lr_sym.py |   2 +
 2 files changed, 268 insertions(+), 102 deletions(-)

diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py
index dd2a625..e0e26f7 100644
--- a/pycc/ccresponse.py
+++ b/pycc/ccresponse.py
@@ -323,6 +323,9 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         l2 = self.cclambda.l2
         hbar = self.hbar
         L = self.ccwfn.H.L
+        F = self.ccwfn.H.F
+        t2 = self.ccwfn.t2
+        ERI = self.ccwfn.H.ERI
 
         # Please refer to eqn 45 of [Crawford: 10.1002/wcms.1406].
         # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
@@ -402,58 +405,84 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         #               Fifth term                          Sixth term
 
         # Expanding the permutational operator and implementing it explicitly
-        temp = contract("ijab, ia -> jb", L[o, o, v, v].copy(), X1_A)
-        polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
-
-        temp = contract("je, ja -> ea", X1_A, hbar.Hov)
-        temp = contract("ea, ma -> me", temp, X1_B)
-        polar -= contract("me, me -> ", temp, l1)
-        temp = contract("ma, me -> ea", X1_A, hbar.Hov)
-        temp = contract("ea, je -> ja", temp, X1_B)
-        polar -= contract("ja, me -> ", temp, l1)
-
-        temp = contract("ijam, je -> ieam", hbar.Hooov.swapaxes(2, 3), X1_A)
-        temp = contract("ieam, ia -> me", temp, X1_B)
-        polar -= contract("me, me -> ", temp, l1)
-        temp = contract("ijem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
-        temp = contract("ajem, je -> ma", temp, X1_B)
-        polar -= contract("ma, ma -> ", temp, l1)
-
-        temp = contract("ejab, jb -> ea", hbar.Hvovv, X1_A)
-        temp = contract("ea, ma -> me", temp, X1_B)
-        polar += contract("me, me -> ", temp, l1)
-        temp = contract("emab, jb -> emaj", hbar.Hvovv, X1_A)
-        temp = contract("emaj, ma -> je", temp, X1_B)
-        polar += contract("je, je -> ", temp, l1)
-
-        temp = contract("jfma, je -> efma", hbar.Hovov, X1_A)
-        temp = contract("efma, na -> mnef", temp, X1_B)
-        polar -= contract("mnef, mnef -> ", temp, l2)
-        temp = contract("nfma, na -> mf", hbar.Hovov, X1_A)
-        temp = contract("mf, je -> mjef", temp, X1_B)
-        polar -= contract("mjef, mjef -> ", temp, l2)
-
-        temp = contract("ejma, jf -> efma", hbar.Hovov.swapaxes(0, 1), X1_A)
-        temp = contract("efma, na -> mnef", temp, X1_B)
-        polar -= contract("mnef, mnef -> ", temp, l2)
-        temp = contract("enmf, na -> amfe", hbar.Hovov.swapaxes(0, 1), X1_A)
-        temp = contract("amfe, jf -> mjae", temp, X1_B)
-        polar -= contract("mjae, mjae -> ", temp, l2)
-
-        temp = contract("abef, nb -> anef", hbar.Hvvvv, X1_A)
-        temp = contract("anef, ma -> mnef", temp, X1_B)
-        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
-        temp = contract("abef, ma -> amef", hbar.Hvvvv, X1_A)
-        temp = contract("amef, nb -> mnef", temp, X1_B)
-        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
-
-        temp = contract("mnij, jf -> mnif", hbar.Hoooo, X1_A)
-        temp = contract("mnif, ie -> mnef", temp, X1_B)
-        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
-        temp = contract("mnij, ie -> mnje", hbar.Hoooo, X1_A)
-        temp = contract("mnje, jf -> mnef", temp, X1_B)
-        polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
-
+        # temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)
+        # polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
+        # # #
+        # temp = contract("je, ja -> ea", X1_A, hbar.Hov)
+        # temp = contract("ea, ma -> me", temp, X1_B)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("je, ja -> ea", X1_B, hbar.Hov)
+        # temp = contract("ea, ma -> me", temp, X1_A)
+        # polar -= contract("me, me -> ", temp, l1)
+        # #
+        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
+        # temp = contract("ieam, ia -> me", temp, X1_B)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_B)
+        # temp = contract("ieam, ia -> me", temp, X1_A)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
+        # temp = contract("ajem, je -> ma", temp, X1_B)
+        # polar += contract("ma, ma -> ", temp, l1)
+        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
+        # temp = contract("ajem, je -> ma", temp, X1_A)
+        # polar += contract("ma, ma -> ", temp, l1)
+        # # #
+        # temp = contract("ejab, jb -> ea", (2.0 * hbar.Hvovv), X1_A)
+        # temp = contract("ea, ma -> me", temp, X1_B)
+        # polar += contract("me, me -> ", temp, l1)
+        # temp = contract("ejab, jb -> ea", hbar.Hvovv.swapaxes(2, 3), X1_A)
+        # temp = contract("ea, ma -> me", temp, X1_B)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("ejba, ja -> be", (2.0 * hbar.Hvovv), X1_B)
+        # temp = contract("be, mb -> me", temp, X1_A)
+        # polar += contract("me, me -> ", temp, l1)
+        # temp = contract("ejba, ja -> be", hbar.Hvovv.swapaxes(2, 3), X1_B)
+        # temp = contract("be, mb -> me", temp, X1_A)
+        # polar -= contract("me, me -> ", temp, l1)
+        #
+        # temp = contract("jfma, je -> amef", hbar.Hovov, X1_A)
+        # temp = contract("amef, na -> mnef", temp, X1_B)
+        # polar -= contract("mnef, mnef -> ", temp, l2)
+        #
+        # temp = contract("jfma, je -> amef", hbar.Hovov, X1_B)
+        # temp = contract("amef, na -> mnef", temp, X1_A)
+        # polar -= contract("mnef, mnef -> ", temp, l2)
+        #
+        # temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_A)
+        # temp = contract("mnej, jc -> mnec", temp, X1_B)
+        # polar -= contract("mnec, mnec -> ", temp, l2)
+        #
+        # temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_B)
+        # temp = contract("mnej, jc -> mnec", temp, X1_A)
+        # polar -= contract("mnec, mnec -> ", temp, l2)
+        #
+        # temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_A)
+        # temp = contract("abej, ie -> ijab", temp, X1_B)
+        # polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
+        #
+        # temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_B)
+        # temp = contract("abej, ie -> ijab", temp, X1_A)
+        # polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
+        #
+        # temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_A)
+        # temp = contract("anij, nb -> ijab", temp, X1_B)
+        # polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
+        #
+        # temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_B)
+        # temp = contract("anij, nb -> ijab", temp, X1_A)
+        # polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
+        #
+        # Goo = contract('mjab,ijab->mi', t2, l2)
+        # Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
+        # r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
+        # r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
+        # r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
+        # r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
+        # polar += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)     #Gvv
+        # # print("Gvv\n", polar_Gvv)
+        # polar += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)     #Goo
+        # # print("Goo\n", polar_Goo)
 
 
 
@@ -475,12 +504,32 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         # temp = contract("ja, jb -> ab", hbar.Hov, X1_A)
         # temp = contract("ab, ia -> ib", temp, X1_B)
         # polar -= 2.0 * contract("ib, ib -> ", l1, temp)
-        # temp = contract("ijka, jc -> icak", hbar.Hooov, X1_A)    # L[o, o, v, o] or hbar.Hoovo
-        # temp = contract("icak, ia -> kc", temp, X1_B)
-        # polar += 2.0 * contract("kc, kc -> ", temp, l1)
-        # temp = contract("cjab, jb -> ca", hbar.Hvovv, X1_A)      ## L[v, o, v, v] or hbar.Hvovv
-        # temp = contract("ca, ka -> kc", temp, X1_B)
-        # polar += 2.0 * contract("kc, kc -> ", temp, l1)
+        #
+        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
+        # temp = contract("ieam, ia -> me", temp, X1_B)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_B)
+        # temp = contract("ieam, ia -> me", temp, X1_A)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
+        # temp = contract("ajem, je -> ma", temp, X1_B)
+        # polar += contract("ma, ma -> ", temp, l1)
+        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
+        # temp = contract("ajem, je -> ma", temp, X1_A)
+        # polar += contract("ma, ma -> ", temp, l1)
+
+        # temp = contract("ejab, jb -> ea", (hbar.Hvovv), X1_A)
+        # temp = contract("ea, ma -> me", temp, X1_B)
+        # polar += contract("me, me -> ", temp, l1)
+        # temp = contract("ejab, ja -> be", hbar.Hvovv, X1_A)
+        # temp = contract("be, mb -> me", temp, X1_B)
+        # polar -= contract("me, me -> ", temp, l1)
+        # temp = contract("ejba, ma -> jmbe", (2 * hbar.Hvvvo.swapaxes(1, 3)), X1_B)
+        # temp = contract("jmbe, jb -> me", temp, X1_A)
+        # polar += contract("me, me -> ", temp, l1)
+        # temp = contract("jeba, ma -> jmbe", (hbar.Hvovv.swapaxes(0, 3).swapaxes(2, 3)), X1_B)
+        # temp = contract("jmbe, jb -> me", temp, X1_A)
+        # polar -= contract("me, me -> ", temp, l1)
         #
         # # Sixth term (contributions from the L2 and X1)
         # temp = contract("je, jfma -> amef", X1_A, hbar.Hovov)
@@ -497,54 +546,169 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         # polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
         # End HX_1Y_1
 
-
-        # Begin HX_{2}Y_{2}
-        # There is no contribution from L1 and X2
-        # Sixth term (contributions from the L2 and X2)
-        # I don't know if I should use the L term or the H_bar term for Hoovv
-        # There is a permutationally
-        temp = contract("klab, mnab -> klmn", L[o, o, v, v], X2_A)
-        temp = contract("klmn, klfe -> mnef", temp, X2_B)
-        polar += 1/16 * contract("mnef, mnef -> ", temp, l2)
-        # temp = contract("klab, nkaf -> lbfn", L[o, o, v, v], X2_A)
-        # temp = contract("lbfn, mlbe -> mnef", temp, X2_B)
-        # polar -= contract("mnef, mnef -> ", temp, l2)
-        # temp = contract("klab, nlae -> kbne", L[o, o, v, v], X2_A)
-        # temp = contract("kbne, mkbf -> mnef", temp, X2_B)
-        # polar += contract("mnef, mnef -> ", temp, l2)
+        # Begin HX_2Y_2
+        # temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
+        # temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
+        # polar += 2 * contract("jlbd, jlbd -> ", temp, l2)
+        #
+        # temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
+        # temp = contract("bc, klcd -> klbd", temp, X2_B)
+        # polar -= contract("klbd, klbd -> ", temp, l2)
+        # temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
+        # temp = contract("jk, jlbd -> klbd", temp, X2_B)
+        # polar -= contract("klbd, klbd -> ", temp, l2)
+        # temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A)
+        # temp = contract("ikab, ilad -> klbd", temp, X2_B)
+        # polar -= contract("klbd, klbd -> ", temp, l2)
+        #
+        # temp = contract("ikad, klcd -> ilac", L[o, o, v, v], X2_B)
+        # temp = contract("ilac, ijab -> jlbc", temp, X2_A)
+        # polar -= contract("jlbc, jlbc -> ", temp, l2)
+        # temp = contract("ikad, ikac -> cd", L[o, o, v, v], X2_B)
+        # temp = contract("cd, jlbd -> jlbc", temp, X2_A)
+        # polar -= contract("jlbc, jlbc -> ", temp, l2)
+        # temp = contract("ikad, ilad -> kl", L[o, o, v, v], X2_B)
+        # temp = contract("kl, jkbc -> jlbc", temp, X2_A)
+        # polar -= contract("jlbc, jlbc -> ", temp, l2)
+        #
+        #
+        # temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A)
+        # temp = contract("klij, klcd -> ijcd", temp, X2_B)
+        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        #
+        # temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_B)
+        # temp = contract("klij, klcd -> ijcd", temp, X2_A)
+        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        #
+        # temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A)
+        # temp = contract("ilbc, jlbd -> ijcd", temp, X2_B)
+        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
         #
-        # # Also, we can consider terms like L1, X1, X2 == L1, X2, X1 or L2, X1, X2
-        # # L1, X1, X2
+        # temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_B)
+        # temp = contract("ilbc, jlbd -> ijcd", temp, X2_A)
+        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        #
+        # temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A)
+        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
+        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        #
+        # temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_B)
+        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
+        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        # End HX_{2}Y_{2}
+
+
+        # Begin L_{1}HX_{1}Y_{2}
+        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
+        # temp = contract("jc, jkbc -> kb", temp, X2_B)
+        # polar -= contract("kb, kb", temp, l1)
+        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        # temp = contract("jc, jkbc -> kb", temp, X2_A)
+        # polar -= contract("kb, kb", temp, l1)
+        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
+        # temp = contract("jc, jkcb -> kb", temp, X2_B)
+        # polar += 2.0 * contract("kb, kb", temp, l1)
         # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
-        # temp = contract("jc, jldc -> ld", temp, X2_A)
-        # polar += -1.0 * contract("ld, ld -> ", temp, l1)
-        # temp = contract("ijac, ijda -> dc", L[o, o, v, v], X2_B)
-        # temp = contract("dc, lc -> ld", temp, X1_A)
-        # polar += contract("ld, ld -> ", temp, l1)
-        # temp = contract("ijac, ilac-> jl", L[o, o, v, v], X2_B)
-        # temp = contract("jl, jd -> ld", temp, X1_A)
-        # polar += -1.0 * contract("ld, ld -> ", temp, l1)
+        # temp = contract("jc, jkcb -> kb", temp, X2_A)
+        # polar += 2.0 * contract("kb, kb", temp, l1)
+        #
+        # temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_B)
+        # temp = contract("jk, jb -> kb", temp, X1_A)
+        # polar -= contract("kb, kb -> ", temp, l1)
+        # temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
+        # temp = contract("jk, jb -> kb", temp, X1_B)
+        # polar -= contract("kb, kb -> ", temp, l1)
+        #
+        # temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_B)
+        # temp = contract("bc, kc -> kb", temp, X1_A)
+        # polar -= contract("kb, kb -> ", temp, l1)
+        # temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
+        # temp = contract("bc, kc -> kb", temp, X1_B)
+        # polar -= contract("kb, kb -> ", temp, l1)
+        # End L_{1}HX_{1}Y_{2}
+
+        # Begin L_{2}HX_{1}Y_{2}
+        # temp = contract("kdab, ia -> ikbd", hbar.Hvovv.swapaxes(0, 1), X1_A)
+        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
+        # polar -= contract("ijcd, ijcd -> ", temp, l2)
+        # temp = contract("kdab, ia -> ikbd", hbar.Hvovv.swapaxes(0, 1), X1_B)
+        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
+        # polar -= contract("ijcd, ijcd -> ", temp, l2)
+
+        # temp = contract("kdab, jb -> jkad", hbar.Hvovv.swapaxes(0, 1), X1_A)
+        # temp = contract("jkad, ikac -> ijcd", temp, X2_B)
+        # polar -= contract("ijcd, ijcd -> ", temp, l2)
+        # temp = contract("kdab, jb -> jkad", hbar.Hvovv.swapaxes(0, 1), X1_B)
+        # temp = contract("jkad, ikac -> ijcd", temp, X2_A)
+        # polar -= contract("ijcd, ijcd -> ", temp, l2)
+        #
+        # temp = contract("kdab, kc -> abcd", hbar.Hvovv.swapaxes(0, 1), X1_A)
+        # temp = contract("abcd, ijab -> ijcd", temp, X2_B)
+        # polar -= contract("ijcd, ijcd -> ", temp, l2)
+        # temp = contract("kdab, kc -> abcd", hbar.Hvovv.swapaxes(0, 1), X1_B)
+        # temp = contract("abcd, ijab -> ijcd", temp, X2_A)
+        # polar -= contract("ijcd, ijcd -> ", temp, l2)
+
+        # temp = contract("kjal, ia -> ijkl", hbar.Hooov.swapaxes(2, 3), X1_A)
+        # temp = contract("ijkl, jkbc -> ilbc", temp, X2_B)
+        # polar += contract("ilbc, ilbc -> ", temp, l2)
+        # temp = contract("kjal, ia -> ijkl", hbar.Hooov.swapaxes(2, 3), X1_B)
+        # temp = contract("ijkl, jkbc -> ilbc", temp, X2_A)
+        # polar += contract("ilbc, ilbc -> ", temp, l2)
+        #
+        # temp = contract("kjal, jb -> klab", hbar.Hooov.swapaxes(2, 3), X1_A)
+        # temp = contract("klab, ikac -> ilbc", temp, X2_B)
+        # polar += contract("ilbc, ilbc -> ", temp, l2)
+        # temp = contract("kjal, jb -> klab", hbar.Hooov.swapaxes(2, 3), X1_B)
+        # temp = contract("klab, ikac -> ilbc", temp, X2_A)
+        # polar += contract("ilbc, ilbc -> ", temp, l2)
+        #
+        # temp = contract("kjal, kc -> jlac", hbar.Hooov.swapaxes(2, 3), X1_A)
+        # temp = contract("jlac, ijab -> ilbc", temp, X2_B)
+        # polar += contract("ilbc, ilbc -> ", temp, l2)
+        # temp = contract("kjal, kc -> jlac", hbar.Hooov.swapaxes(2, 3), X1_B)
+        # temp = contract("jlac, ijab -> ilbc", temp, X2_A)
+        # polar += contract("ilbc, ilbc -> ", temp, l2)
+
+        # # temp = contract("ja, ia -> ij", hbar.Hov, X1_A)
+        # # temp = contract("ij, jkbc -> ikbc", temp, X2_B)
+        # # polar -= contract("ikbc, ikbc -> ", temp, l2)
+        # # temp = contract("ib, ia -> ab", hbar.Hov, X1_A)
+        # # temp = contract("ab, jkbc -> jkac", temp, X2_B)
+        # # polar -= contract("jkac, jkac -> ", temp, l2)
+        # # temp = contract("ja, ia -> ij", hbar.Hov, X1_B)
+        # # temp = contract("ij, jkbc -> ikbc", temp, X2_A)
+        # # polar -= contract("ikbc, ikbc -> ", temp, l2)
+        # # temp = contract("ib, ia -> ab", hbar.Hov, X1_B)
+        # # temp = contract("ab, jkbc -> jkac", temp, X2_A)
+        # # polar -= contract("jkac, jkac -> ", temp, l2)
         #
-        # # L2, X1, X2
-        # temp = contract("dkab, ma -> dkmb", hbar.Hvovv, X1_B)
-        # temp = contract("dkmb, lkbe -> lmde", temp, X2_A)
-        # polar += 0.5 * contract("lmde, lmde -> ", temp, l2)
-        # temp = contract("dkab, mkae -> dmbe", hbar.Hvovv, X2_B)
-        # temp = contract("dmbe, lb -> dmel", temp, X1_A)
-        # polar -= 0.5 * contract("dmel, lmde -> ", temp, l2)
-        # temp = contract("dkab, mlab -> mldk", hbar.Hvovv, X2_B)
-        # temp = contract("mldk, ke -> dmel", temp, X1_A)
-        # polar += 0.5 * contract("dmel, lmde -> ", temp, l2)
+        # temp = contract("ijam, ia -> jm", (2 * hbar.Hooov.swapaxes(2, 3)), X1_A)
+        # temp = contract("jm, jkbc -> mkbc", temp, X2_B)
+        # polar -= contract("mkbc, mkbc -> ", temp, l2)
+        # temp = contract("ijam, ia -> jm", (2 * hbar.Hooov.swapaxes(2, 3)), X1_B)
+        # temp = contract("jm, jkbc -> mkbc", temp, X2_A)
+        # polar -= contract("mkbc, mkbc -> ", temp, l2)
+        # temp = contract("ijma, ia -> jm", hbar.Hooov, X1_A)
+        # temp = contract("jm, jkbc -> mkbc", temp, X2_B)
+        # polar += contract("mkbc, mkbc -> ", temp, l2)
+        # temp = contract("ijma, ia -> jm", hbar.Hooov, X1_B)
+        # temp = contract("jm, jkbc -> mkbc", temp, X2_A)
+        # polar += contract("mkbc, mkbc -> ", temp, l2)
         #
-        # temp = contract("kjla, klde -> ajde", hbar.Hooov, X2_B)
-        # temp = contract("ajde, ia -> ijde", temp, X1_A)
-        # polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
-        # temp = contract("klja, ikae -> ijle", hbar.Hooov, X2_B)
-        # temp = contract("ijle, ld -> ijde", temp, X1_A)
-        # polar += 0.5 * contract("ijde, ijde -> ", temp, l2)
-        # temp = contract("klja, ilad -> ijkd", hbar.Hooov, X2_B)
-        # temp = contract("ijkd, ke -> ijde", temp, X1_A)
-        # polar -= 0.5 * contract("ijde, ijde -> ", temp, l2)
+        # temp = contract("djab, ia -> ijbd", (2 * hbar.Hvovv), X1_A)
+        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_B)
+        # polar += contract("ikcd, ikcd -> ", temp, l2)
+        # temp = contract("djab, ia -> ijbd", (2 * hbar.Hvovv), X1_B)
+        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_A)
+        # polar += contract("ikcd, ikcd -> ", temp, l2)
+        # temp = contract("djba, ia -> ijbd", hbar.Hvovv, X1_A)
+        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_B)
+        # polar -= contract("ikcd, ikcd -> ", temp, l2)
+        # temp = contract("djba, ia -> ijbd", hbar.Hvovv, X1_B)
+        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_A)
+        # polar -= contract("ikcd, ikcd -> ", temp, l2)
+        # End L_{2}HX_{1}Y_{2}
 
         return polar#-1.0 * polar
 
diff --git a/pycc/tests/test_036_lr_sym.py b/pycc/tests/test_036_lr_sym.py
index c0354d3..5084765 100644
--- a/pycc/tests/test_036_lr_sym.py
+++ b/pycc/tests/test_036_lr_sym.py
@@ -112,6 +112,8 @@
     string_a = "MU_" + resp.cart[a]
     X1_A, X2_A, _ = X_A[string_a]
     for b in range(0, 3):
+        # string_a = "MU_" + resp.cart[a]
+        # X1_A, X2_A, _ = X_A[string_a]
         string_b = "MU_" + resp.cart[b]
         # Y1_B, Y2_B, _ = Y_2[string_b]
         # X1_B, X2_B, _ = X_2[string_b]

From e6ca86232d1924d50e3eadc0a814edae92d77bbf Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Thu, 5 Sep 2024 17:56:47 -0400
Subject: [PATCH 4/7] Linear response symmetric

---
 pycc/ccresponse.py            | 689 ++++++++++++++++------------------
 pycc/data/molecules.py        |   9 +
 pycc/tests/test_036_lr_sym.py |  44 +--
 3 files changed, 365 insertions(+), 377 deletions(-)

diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py
index e0e26f7..5700ffe 100644
--- a/pycc/ccresponse.py
+++ b/pycc/ccresponse.py
@@ -340,374 +340,353 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         pertbar_B = self.pertbar[pertkey_b]
         Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
 
-        # # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
-        # polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
-        # # Second term (contributions from the L1 and X1)
-        # temp = contract("ia, ic -> ac", l1, X1_B)
-        # polar += contract("ac, ac -> ", temp, pertbar_A.Avv)
-        # temp = contract("ia, ka -> ik", l1, X1_B)
-        # polar -= contract("ik, ki -> ", temp, pertbar_A.Aoo)
-        # # Second term (contributions from the L1 and X2)
-        # temp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
-        # polar += 2.0 * contract("ijab, ijab -> ", temp, X2_B)
-        # polar += -1.0 * contract("ijab, ijba -> ", temp, X2_B)
-        # # Second term (contributions from the L2 and X1)
-        # temp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
-        # polar += contract("ia, ia -> ", temp, X1_B)
-        # temp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
-        # polar -= 0.5 * contract("ak, ka -> ", temp, X1_B)
-        # temp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
-        # polar -= 0.5 * contract("bk, kb -> ", temp, X1_B)
-        # # Second term (contributions from the L2 and X2)
-        # temp = contract("ijab, kjab -> ik", l2, X2_B)
-        # polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_A.Aoo)
-        # temp = contract("ijab, kiba-> jk", l2, X2_B)
-        # polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_A.Aoo)
-        # temp = contract("ijab, ijac -> bc", l2, X2_B)
-        # polar += 0.5 * contract("bc, bc -> ", temp, pertbar_A.Avv)
-        # temp = contract("ijab, ijcb -> ac", l2, X2_B)
-        # polar += 0.5 * contract("ac, ac -> ", temp, pertbar_A.Avv)
-        #
-        # # <0|(1 + L^(0)) [\bar{B}^(1), X^(1)(B)] |0>
-        # # <0| [\bar{B}^(1), X^(1)(B)] |0> + <0| L^(0) [\bar{B}^(0), X^(1)(B)] |0>
-        # #           Third term                      Fourth term
-        #
-        # # Third term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
-        # polar += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
-        # # Fourth term (contributions from the L1 and X1)
-        # temp = contract("ia, ic -> ac", l1, X1_A)
-        # polar += contract("ac, ac -> ", temp, pertbar_B.Avv)
-        # temp = contract("ia, ka -> ik", l1, X1_A)
-        # polar -= contract("ik, ki -> ", temp, pertbar_B.Aoo)
-        # # Fourth term (contributions from the L1 and X2)
-        # temp = contract("ia, jb -> ijab", l1, pertbar_B.Aov)
-        # polar += 2.0 * contract("ijab, ijab -> ", temp, X2_A)
-        # polar += -1.0 * contract("ijab, ijba -> ", temp, X2_A)
-        # # Fourth term (contributions from the L2 and X1)
-        # temp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
-        # polar += contract("ia, ia -> ", temp, X1_A)
-        # temp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
-        # polar -= 0.5 * contract("ak, ka -> ", temp, X1_A)
-        # temp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
-        # polar -= 0.5 * contract("bk, kb -> ", temp, X1_A)
-        # # Fourth term (contributions from the L2 and X2)
-        # temp = contract("ijab, kjab -> ik", l2, X2_A)
-        # polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_B.Aoo)
-        # temp = contract("ijab, kiba-> jk", l2, X2_A)
-        # polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_B.Aoo)
-        # temp = contract("ijab, ijac -> bc", l2, X2_A)
-        # polar += 0.5 * contract("bc, bc -> ", temp, pertbar_B.Avv)
-        # temp = contract("ijab, ijcb -> ac", l2, X2_A)
-        # polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
+        # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
+        polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
+        # Second term (contributions from the L1 and X1)
+        temp = contract("ia, ic -> ac", l1, X1_B)
+        polar += contract("ac, ac -> ", temp, pertbar_A.Avv)
+        temp = contract("ia, ka -> ik", l1, X1_B)
+        polar -= contract("ik, ki -> ", temp, pertbar_A.Aoo)
+        # Second term (contributions from the L1 and X2)
+        temp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
+        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_B)
+        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_B)
+        # Second term (contributions from the L2 and X1)
+        temp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
+        polar += contract("ia, ia -> ", temp, X1_B)
+        temp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
+        polar -= 0.5 * contract("ak, ka -> ", temp, X1_B)
+        temp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
+        polar -= 0.5 * contract("bk, kb -> ", temp, X1_B)
+        # Second term (contributions from the L2 and X2)
+        temp = contract("ijab, kjab -> ik", l2, X2_B)
+        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_A.Aoo)
+        temp = contract("ijab, kiba-> jk", l2, X2_B)
+        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_A.Aoo)
+        temp = contract("ijab, ijac -> bc", l2, X2_B)
+        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_A.Avv)
+        temp = contract("ijab, ijcb -> ac", l2, X2_B)
+        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_A.Avv)
+
+        # <0|(1 + L^(0)) [\bar{B}^(1), X^(1)(B)] |0>
+        # <0| [\bar{B}^(1), X^(1)(B)] |0> + <0| L^(0) [\bar{B}^(0), X^(1)(B)] |0>
+        #           Third term                      Fourth term
+
+        # Third term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
+        polar += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
+        # Fourth term (contributions from the L1 and X1)
+        temp = contract("ia, ic -> ac", l1, X1_A)
+        polar += contract("ac, ac -> ", temp, pertbar_B.Avv)
+        temp = contract("ia, ka -> ik", l1, X1_A)
+        polar -= contract("ik, ki -> ", temp, pertbar_B.Aoo)
+        # Fourth term (contributions from the L1 and X2)
+        temp = contract("ia, jb -> ijab", l1, pertbar_B.Aov)
+        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_A)
+        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_A)
+        # Fourth term (contributions from the L2 and X1)
+        temp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
+        polar += contract("ia, ia -> ", temp, X1_A)
+        temp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
+        polar -= 0.5 * contract("ak, ka -> ", temp, X1_A)
+        temp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
+        polar -= 0.5 * contract("bk, kb -> ", temp, X1_A)
+        # Fourth term (contributions from the L2 and X2)
+        temp = contract("ijab, kjab -> ik", l2, X2_A)
+        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_B.Aoo)
+        temp = contract("ijab, kiba-> jk", l2, X2_A)
+        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_B.Aoo)
+        temp = contract("ijab, ijac -> bc", l2, X2_A)
+        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_B.Avv)
+        temp = contract("ijab, ijcb -> ac", l2, X2_A)
+        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
 
         # <0|(1 + L^(0))[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
         # <0|[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0> + <0| L^(0) [[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
         #               Fifth term                          Sixth term
 
         # Expanding the permutational operator and implementing it explicitly
-        # temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)
-        # polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
-        # # #
-        # temp = contract("je, ja -> ea", X1_A, hbar.Hov)
-        # temp = contract("ea, ma -> me", temp, X1_B)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("je, ja -> ea", X1_B, hbar.Hov)
-        # temp = contract("ea, ma -> me", temp, X1_A)
-        # polar -= contract("me, me -> ", temp, l1)
-        # #
-        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
-        # temp = contract("ieam, ia -> me", temp, X1_B)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_B)
-        # temp = contract("ieam, ia -> me", temp, X1_A)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
-        # temp = contract("ajem, je -> ma", temp, X1_B)
-        # polar += contract("ma, ma -> ", temp, l1)
-        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
-        # temp = contract("ajem, je -> ma", temp, X1_A)
-        # polar += contract("ma, ma -> ", temp, l1)
+        temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)
+        polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
         # # #
-        # temp = contract("ejab, jb -> ea", (2.0 * hbar.Hvovv), X1_A)
-        # temp = contract("ea, ma -> me", temp, X1_B)
-        # polar += contract("me, me -> ", temp, l1)
-        # temp = contract("ejab, jb -> ea", hbar.Hvovv.swapaxes(2, 3), X1_A)
-        # temp = contract("ea, ma -> me", temp, X1_B)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("ejba, ja -> be", (2.0 * hbar.Hvovv), X1_B)
-        # temp = contract("be, mb -> me", temp, X1_A)
-        # polar += contract("me, me -> ", temp, l1)
-        # temp = contract("ejba, ja -> be", hbar.Hvovv.swapaxes(2, 3), X1_B)
-        # temp = contract("be, mb -> me", temp, X1_A)
-        # polar -= contract("me, me -> ", temp, l1)
-        #
-        # temp = contract("jfma, je -> amef", hbar.Hovov, X1_A)
-        # temp = contract("amef, na -> mnef", temp, X1_B)
-        # polar -= contract("mnef, mnef -> ", temp, l2)
-        #
-        # temp = contract("jfma, je -> amef", hbar.Hovov, X1_B)
-        # temp = contract("amef, na -> mnef", temp, X1_A)
-        # polar -= contract("mnef, mnef -> ", temp, l2)
-        #
-        # temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_A)
-        # temp = contract("mnej, jc -> mnec", temp, X1_B)
-        # polar -= contract("mnec, mnec -> ", temp, l2)
-        #
-        # temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_B)
-        # temp = contract("mnej, jc -> mnec", temp, X1_A)
-        # polar -= contract("mnec, mnec -> ", temp, l2)
-        #
-        # temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_A)
-        # temp = contract("abej, ie -> ijab", temp, X1_B)
-        # polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
-        #
-        # temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_B)
-        # temp = contract("abej, ie -> ijab", temp, X1_A)
-        # polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
-        #
-        # temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_A)
-        # temp = contract("anij, nb -> ijab", temp, X1_B)
-        # polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
-        #
-        # temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_B)
-        # temp = contract("anij, nb -> ijab", temp, X1_A)
-        # polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
+        temp = contract("je, ja -> ea", X1_A, hbar.Hov)
+        temp = contract("ea, ma -> me", temp, X1_B)
+        polar -= contract("me, me -> ", temp, l1)
+        temp = contract("je, ja -> ea", X1_B, hbar.Hov)
+        temp = contract("ea, ma -> me", temp, X1_A)
+        polar -= contract("me, me -> ", temp, l1)
         #
-        # Goo = contract('mjab,ijab->mi', t2, l2)
-        # Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
-        # r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
-        # r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
-        # r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
-        # r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
-        # polar += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)     #Gvv
-        # # print("Gvv\n", polar_Gvv)
-        # polar += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)     #Goo
-        # # print("Goo\n", polar_Goo)
-
-
-
-
-
-
-
-
-
-
-
-
-        # Begin HX_1Y_1
-        # Fifth term (Contribution from T1, there's no contribution from T2)
-        # temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)  # L[o, o, v, v] or hbar.Hoovv
-        # polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
-
-        # Sixth term (contributions from the L1 and X1, also, no contributing term from L1 and X2)
-        # temp = contract("ja, jb -> ab", hbar.Hov, X1_A)
-        # temp = contract("ab, ia -> ib", temp, X1_B)
-        # polar -= 2.0 * contract("ib, ib -> ", l1, temp)
-        #
-        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
-        # temp = contract("ieam, ia -> me", temp, X1_B)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_B)
-        # temp = contract("ieam, ia -> me", temp, X1_A)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
-        # temp = contract("ajem, je -> ma", temp, X1_B)
-        # polar += contract("ma, ma -> ", temp, l1)
-        # temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
-        # temp = contract("ajem, je -> ma", temp, X1_A)
-        # polar += contract("ma, ma -> ", temp, l1)
-
-        # temp = contract("ejab, jb -> ea", (hbar.Hvovv), X1_A)
-        # temp = contract("ea, ma -> me", temp, X1_B)
-        # polar += contract("me, me -> ", temp, l1)
-        # temp = contract("ejab, ja -> be", hbar.Hvovv, X1_A)
-        # temp = contract("be, mb -> me", temp, X1_B)
-        # polar -= contract("me, me -> ", temp, l1)
-        # temp = contract("ejba, ma -> jmbe", (2 * hbar.Hvvvo.swapaxes(1, 3)), X1_B)
-        # temp = contract("jmbe, jb -> me", temp, X1_A)
-        # polar += contract("me, me -> ", temp, l1)
-        # temp = contract("jeba, ma -> jmbe", (hbar.Hvovv.swapaxes(0, 3).swapaxes(2, 3)), X1_B)
-        # temp = contract("jmbe, jb -> me", temp, X1_A)
-        # polar -= contract("me, me -> ", temp, l1)
-        #
-        # # Sixth term (contributions from the L2 and X1)
-        # temp = contract("je, jfma -> amef", X1_A, hbar.Hovov)
-        # temp = contract("amef, na -> mnef", temp, X1_B)
-        # polar -= contract("mnef, mnef -> ", temp, l2)
-        # temp = contract("je, fjma -> amef", X1_A, hbar.Hovov.swapaxes(0, 1))
-        # temp = contract("amef, na -> mnef", temp, X1_B)
-        # polar -= contract("mnef, mnef -> ", temp, l2)
-        # temp = contract("nb, abef -> anef", X1_A, hbar.Hvvvv)
-        # temp = contract("anef, ma -> mnef", temp, X1_B)
-        # polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
-        # temp = contract("jf, mnij -> mnif", X1_A, hbar.Hoooo)
-        # temp = contract("mnif, ie -> mnef", temp, X1_B)
-        # polar += 0.5 * contract("mnef, mnef -> ", temp, l2)
-        # End HX_1Y_1
-
-        # Begin HX_2Y_2
-        # temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
-        # temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
-        # polar += 2 * contract("jlbd, jlbd -> ", temp, l2)
-        #
-        # temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
-        # temp = contract("bc, klcd -> klbd", temp, X2_B)
-        # polar -= contract("klbd, klbd -> ", temp, l2)
-        # temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
-        # temp = contract("jk, jlbd -> klbd", temp, X2_B)
-        # polar -= contract("klbd, klbd -> ", temp, l2)
-        # temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A)
-        # temp = contract("ikab, ilad -> klbd", temp, X2_B)
-        # polar -= contract("klbd, klbd -> ", temp, l2)
-        #
-        # temp = contract("ikad, klcd -> ilac", L[o, o, v, v], X2_B)
-        # temp = contract("ilac, ijab -> jlbc", temp, X2_A)
-        # polar -= contract("jlbc, jlbc -> ", temp, l2)
-        # temp = contract("ikad, ikac -> cd", L[o, o, v, v], X2_B)
-        # temp = contract("cd, jlbd -> jlbc", temp, X2_A)
-        # polar -= contract("jlbc, jlbc -> ", temp, l2)
-        # temp = contract("ikad, ilad -> kl", L[o, o, v, v], X2_B)
-        # temp = contract("kl, jkbc -> jlbc", temp, X2_A)
-        # polar -= contract("jlbc, jlbc -> ", temp, l2)
-        #
-        #
-        # temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A)
-        # temp = contract("klij, klcd -> ijcd", temp, X2_B)
-        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        #
-        # temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_B)
-        # temp = contract("klij, klcd -> ijcd", temp, X2_A)
-        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        #
-        # temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A)
-        # temp = contract("ilbc, jlbd -> ijcd", temp, X2_B)
-        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        #
-        # temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_B)
-        # temp = contract("ilbc, jlbd -> ijcd", temp, X2_A)
-        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        #
-        # temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A)
-        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
-        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        #
-        # temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_B)
-        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
-        # polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        # End HX_{2}Y_{2}
-
+        temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
+        temp = contract("ieam, ia -> me", temp, X1_B)
+        polar -= contract("me, me -> ", temp, l1)
+        temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_B)
+        temp = contract("ieam, ia -> me", temp, X1_A)
+        polar -= contract("me, me -> ", temp, l1)
+        temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
+        temp = contract("ajem, je -> ma", temp, X1_B)
+        polar += contract("ma, ma -> ", temp, l1)
+        temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
+        temp = contract("ajem, je -> ma", temp, X1_A)
+        polar += contract("ma, ma -> ", temp, l1)
+        # #
+        temp = contract("ejab, jb -> ea", (2.0 * hbar.Hvovv), X1_A)
+        temp = contract("ea, ma -> me", temp, X1_B)
+        polar += contract("me, me -> ", temp, l1)
+        temp = contract("ejab, jb -> ea", hbar.Hvovv.swapaxes(2, 3), X1_A)
+        temp = contract("ea, ma -> me", temp, X1_B)
+        polar -= contract("me, me -> ", temp, l1)
+        temp = contract("ejba, ja -> be", (2.0 * hbar.Hvovv), X1_B)
+        temp = contract("be, mb -> me", temp, X1_A)
+        polar += contract("me, me -> ", temp, l1)
+        temp = contract("ejba, ja -> be", hbar.Hvovv.swapaxes(2, 3), X1_B)
+        temp = contract("be, mb -> me", temp, X1_A)
+        polar -= contract("me, me -> ", temp, l1)
+
+        temp = contract("jfma, je -> amef", hbar.Hovov, X1_A)
+        temp = contract("amef, na -> mnef", temp, X1_B)
+        polar -= contract("mnef, mnef -> ", temp, l2)
+
+        temp = contract("jfma, je -> amef", hbar.Hovov, X1_B)
+        temp = contract("amef, na -> mnef", temp, X1_A)
+        polar -= contract("mnef, mnef -> ", temp, l2)
+
+        temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_A)
+        temp = contract("mnej, jc -> mnec", temp, X1_B)
+        polar -= contract("mnec, mnec -> ", temp, l2)
+
+        temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_B)
+        temp = contract("mnej, jc -> mnec", temp, X1_A)
+        polar -= contract("mnec, mnec -> ", temp, l2)
+
+        temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_A)
+        temp = contract("abej, ie -> ijab", temp, X1_B)
+        polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
+
+        temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_B)
+        temp = contract("abej, ie -> ijab", temp, X1_A)
+        polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
+
+        temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_A)
+        temp = contract("anij, nb -> ijab", temp, X1_B)
+        polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
+
+        temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_B)
+        temp = contract("anij, nb -> ijab", temp, X1_A)
+        polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
+
+        Goo = contract('mjab,ijab->mi', t2, l2)
+        Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
+        r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
+        r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
+        r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
+        r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
+        polar += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)     #Gvv
+        # print("Gvv\n", polar_Gvv)
+        polar += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)     #Goo
+        # print("Goo\n", polar_Goo)
+
+        # # Begin HX_2Y_2
+        temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
+        temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
+        polar += 2 * contract("jlbd, jlbd -> ", temp, l2)
+
+        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
+        temp = contract("bc, klcd -> klbd", temp, X2_B)
+        polar -= contract("klbd, klbd -> ", temp, l2)
+        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
+        temp = contract("jk, jlbd -> klbd", temp, X2_B)
+        polar -= contract("klbd, klbd -> ", temp, l2)
+        temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A)
+        temp = contract("ikab, ilad -> klbd", temp, X2_B)
+        polar -= contract("klbd, klbd -> ", temp, l2)
+
+        temp = contract("ikad, klcd -> ilac", L[o, o, v, v], X2_B)
+        temp = contract("ilac, ijab -> jlbc", temp, X2_A)
+        polar -= contract("jlbc, jlbc -> ", temp, l2)
+        temp = contract("ikad, ikac -> cd", L[o, o, v, v], X2_B)
+        temp = contract("cd, jlbd -> jlbc", temp, X2_A)
+        polar -= contract("jlbc, jlbc -> ", temp, l2)
+        temp = contract("ikad, ilad -> kl", L[o, o, v, v], X2_B)
+        temp = contract("kl, jkbc -> jlbc", temp, X2_A)
+        polar -= contract("jlbc, jlbc -> ", temp, l2)
+
+
+        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A)
+        temp = contract("klij, klcd -> ijcd", temp, X2_B)
+        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+
+        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_B)
+        temp = contract("klij, klcd -> ijcd", temp, X2_A)
+        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+
+        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A)
+        temp = contract("ilbc, jlbd -> ijcd", temp, X2_B)
+        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+
+        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_B)
+        temp = contract("ilbc, jlbd -> ijcd", temp, X2_A)
+        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+
+        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A)
+        temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
+        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+
+        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_B)
+        temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
+        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        # # End HX_{2}Y_{2}
 
         # Begin L_{1}HX_{1}Y_{2}
-        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
-        # temp = contract("jc, jkbc -> kb", temp, X2_B)
-        # polar -= contract("kb, kb", temp, l1)
-        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
-        # temp = contract("jc, jkbc -> kb", temp, X2_A)
-        # polar -= contract("kb, kb", temp, l1)
-        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
-        # temp = contract("jc, jkcb -> kb", temp, X2_B)
-        # polar += 2.0 * contract("kb, kb", temp, l1)
-        # temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
-        # temp = contract("jc, jkcb -> kb", temp, X2_A)
-        # polar += 2.0 * contract("kb, kb", temp, l1)
-        #
-        # temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_B)
-        # temp = contract("jk, jb -> kb", temp, X1_A)
-        # polar -= contract("kb, kb -> ", temp, l1)
-        # temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
-        # temp = contract("jk, jb -> kb", temp, X1_B)
-        # polar -= contract("kb, kb -> ", temp, l1)
-        #
-        # temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_B)
-        # temp = contract("bc, kc -> kb", temp, X1_A)
-        # polar -= contract("kb, kb -> ", temp, l1)
-        # temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
-        # temp = contract("bc, kc -> kb", temp, X1_B)
-        # polar -= contract("kb, kb -> ", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
+        temp = contract("jc, jkbc -> kb", temp, X2_B)
+        polar -= contract("kb, kb", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        temp = contract("jc, jkbc -> kb", temp, X2_A)
+        polar -= contract("kb, kb", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
+        temp = contract("jc, jkcb -> kb", temp, X2_B)
+        polar += 2.0 * contract("kb, kb", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        temp = contract("jc, jkcb -> kb", temp, X2_A)
+        polar += 2.0 * contract("kb, kb", temp, l1)
+
+        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_B)
+        temp = contract("jk, jb -> kb", temp, X1_A)
+        polar -= contract("kb, kb -> ", temp, l1)
+        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
+        temp = contract("jk, jb -> kb", temp, X1_B)
+        polar -= contract("kb, kb -> ", temp, l1)
+
+        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_B)
+        temp = contract("bc, kc -> kb", temp, X1_A)
+        polar -= contract("kb, kb -> ", temp, l1)
+        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
+        temp = contract("bc, kc -> kb", temp, X1_B)
+        polar -= contract("kb, kb -> ", temp, l1)
         # End L_{1}HX_{1}Y_{2}
 
         # Begin L_{2}HX_{1}Y_{2}
-        # temp = contract("kdab, ia -> ikbd", hbar.Hvovv.swapaxes(0, 1), X1_A)
-        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
-        # polar -= contract("ijcd, ijcd -> ", temp, l2)
-        # temp = contract("kdab, ia -> ikbd", hbar.Hvovv.swapaxes(0, 1), X1_B)
-        # temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
-        # polar -= contract("ijcd, ijcd -> ", temp, l2)
-
-        # temp = contract("kdab, jb -> jkad", hbar.Hvovv.swapaxes(0, 1), X1_A)
-        # temp = contract("jkad, ikac -> ijcd", temp, X2_B)
-        # polar -= contract("ijcd, ijcd -> ", temp, l2)
-        # temp = contract("kdab, jb -> jkad", hbar.Hvovv.swapaxes(0, 1), X1_B)
-        # temp = contract("jkad, ikac -> ijcd", temp, X2_A)
-        # polar -= contract("ijcd, ijcd -> ", temp, l2)
-        #
-        # temp = contract("kdab, kc -> abcd", hbar.Hvovv.swapaxes(0, 1), X1_A)
-        # temp = contract("abcd, ijab -> ijcd", temp, X2_B)
-        # polar -= contract("ijcd, ijcd -> ", temp, l2)
-        # temp = contract("kdab, kc -> abcd", hbar.Hvovv.swapaxes(0, 1), X1_B)
-        # temp = contract("abcd, ijab -> ijcd", temp, X2_A)
-        # polar -= contract("ijcd, ijcd -> ", temp, l2)
-
-        # temp = contract("kjal, ia -> ijkl", hbar.Hooov.swapaxes(2, 3), X1_A)
-        # temp = contract("ijkl, jkbc -> ilbc", temp, X2_B)
-        # polar += contract("ilbc, ilbc -> ", temp, l2)
-        # temp = contract("kjal, ia -> ijkl", hbar.Hooov.swapaxes(2, 3), X1_B)
-        # temp = contract("ijkl, jkbc -> ilbc", temp, X2_A)
-        # polar += contract("ilbc, ilbc -> ", temp, l2)
-        #
-        # temp = contract("kjal, jb -> klab", hbar.Hooov.swapaxes(2, 3), X1_A)
-        # temp = contract("klab, ikac -> ilbc", temp, X2_B)
-        # polar += contract("ilbc, ilbc -> ", temp, l2)
-        # temp = contract("kjal, jb -> klab", hbar.Hooov.swapaxes(2, 3), X1_B)
-        # temp = contract("klab, ikac -> ilbc", temp, X2_A)
-        # polar += contract("ilbc, ilbc -> ", temp, l2)
-        #
-        # temp = contract("kjal, kc -> jlac", hbar.Hooov.swapaxes(2, 3), X1_A)
-        # temp = contract("jlac, ijab -> ilbc", temp, X2_B)
-        # polar += contract("ilbc, ilbc -> ", temp, l2)
-        # temp = contract("kjal, kc -> jlac", hbar.Hooov.swapaxes(2, 3), X1_B)
-        # temp = contract("jlac, ijab -> ilbc", temp, X2_A)
-        # polar += contract("ilbc, ilbc -> ", temp, l2)
-
-        # # temp = contract("ja, ia -> ij", hbar.Hov, X1_A)
-        # # temp = contract("ij, jkbc -> ikbc", temp, X2_B)
-        # # polar -= contract("ikbc, ikbc -> ", temp, l2)
-        # # temp = contract("ib, ia -> ab", hbar.Hov, X1_A)
-        # # temp = contract("ab, jkbc -> jkac", temp, X2_B)
-        # # polar -= contract("jkac, jkac -> ", temp, l2)
-        # # temp = contract("ja, ia -> ij", hbar.Hov, X1_B)
-        # # temp = contract("ij, jkbc -> ikbc", temp, X2_A)
-        # # polar -= contract("ikbc, ikbc -> ", temp, l2)
-        # # temp = contract("ib, ia -> ab", hbar.Hov, X1_B)
-        # # temp = contract("ab, jkbc -> jkac", temp, X2_A)
-        # # polar -= contract("jkac, jkac -> ", temp, l2)
-        #
-        # temp = contract("ijam, ia -> jm", (2 * hbar.Hooov.swapaxes(2, 3)), X1_A)
-        # temp = contract("jm, jkbc -> mkbc", temp, X2_B)
-        # polar -= contract("mkbc, mkbc -> ", temp, l2)
-        # temp = contract("ijam, ia -> jm", (2 * hbar.Hooov.swapaxes(2, 3)), X1_B)
-        # temp = contract("jm, jkbc -> mkbc", temp, X2_A)
-        # polar -= contract("mkbc, mkbc -> ", temp, l2)
-        # temp = contract("ijma, ia -> jm", hbar.Hooov, X1_A)
-        # temp = contract("jm, jkbc -> mkbc", temp, X2_B)
-        # polar += contract("mkbc, mkbc -> ", temp, l2)
-        # temp = contract("ijma, ia -> jm", hbar.Hooov, X1_B)
-        # temp = contract("jm, jkbc -> mkbc", temp, X2_A)
-        # polar += contract("mkbc, mkbc -> ", temp, l2)
-        #
-        # temp = contract("djab, ia -> ijbd", (2 * hbar.Hvovv), X1_A)
-        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_B)
-        # polar += contract("ikcd, ikcd -> ", temp, l2)
-        # temp = contract("djab, ia -> ijbd", (2 * hbar.Hvovv), X1_B)
-        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_A)
-        # polar += contract("ikcd, ikcd -> ", temp, l2)
-        # temp = contract("djba, ia -> ijbd", hbar.Hvovv, X1_A)
-        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_B)
-        # polar -= contract("ikcd, ikcd -> ", temp, l2)
-        # temp = contract("djba, ia -> ijbd", hbar.Hvovv, X1_B)
-        # temp = contract("ijbd, jkbc -> ikcd", temp, X2_A)
-        # polar -= contract("ikcd, ikcd -> ", temp, l2)
+        # <O|L2(A)[[Hbar(0),X2(B)],X1(C)]]|0>
+        tmp = contract("klcd,lkdb->cb", X2_B, l2)
+        tmp = contract("jb,cb->jc", X1_A, tmp)
+        polar -= contract("jc,jc->", tmp, hbar.Hov)
+
+        tmp = contract("klcd,ljdc->kj", X2_B, l2)
+        tmp = contract("kj,jb->kb", tmp, X1_A)
+        polar -= contract("kb,kb->", tmp, hbar.Hov)
+
+        # down
+        tmp = contract('lkda,klcd->ac', l2, X2_B)
+        tmp2 = contract('jb,ajcb->ac', X1_A, hbar.Hvovv)
+        polar += 2.0 * contract('ac,ac->', tmp, tmp2)
+
+        tmp = contract('lkda,klcd->ac', l2, X2_B)
+        tmp2 = contract('jb,ajbc->ac', X1_A, hbar.Hvovv)
+        polar -= contract('ac,ac->', tmp, tmp2)
+
+        tmp = contract('jb,ljda->lbda', X1_A, l2)
+
+        # swapaxes
+        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_B, hbar.Hvovv)
+        tmp2 -= contract('klcd,akcb->ldab', X2_B, hbar.Hvovv)
+        polar += contract('lbda,ldab->', tmp, tmp2)
+
+        tmp = contract('ia,fkba->fkbi', X1_A, hbar.Hvovv)
+        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
+        polar -= contract('jkbc,kjbc->', X2_B, tmp)
+
+        tmp = contract('ia,fjac->fjic', X1_A, hbar.Hvovv)
+        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
+        polar -= contract('jkbc,jkbc->', X2_B, tmp)
+
+        tmp = contract('ia,jkfa->jkfi', X1_A, l2)
+        tmp2 = contract('jkbc,fibc->jkfi', X2_B, hbar.Hvovv)
+        polar -= contract('jkfi,jkfi->', tmp2, tmp)
+
+        tmp = contract('jb,kjib->ki', X1_A, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
+        polar -= 2.0 * contract('ki,ki->', tmp, tmp2)
+
+        tmp = contract('jb,jkib->ki', X1_A, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
+        polar += contract('ki,ki->', tmp, tmp2)
+
+        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_B)
+        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_B)
+        tmp = contract('jild,jb->bild', tmp, X1_A)
+        polar -= contract('bild,ilbd->', tmp, l2)
+
+        tmp = contract('ia,jkna->jkni', X1_A, hbar.Hooov)
+        tmp2 = contract('jkbc,nibc->jkni', X2_B, l2)
+        polar += contract('jkni,jkni->', tmp2, tmp)
+
+        tmp = contract('ia,nkab->nkib', X1_A, l2)
+        tmp = contract('jkbc,nkib->jnic', X2_B, tmp)
+        polar += contract('jnic,ijnc->', tmp, hbar.Hooov)
+
+        tmp = contract('ia,nkba->nkbi', X1_A, l2)
+        tmp = contract('jkbc,nkbi->jnci', X2_B, tmp)
+        polar += contract('jnci,jinc->', tmp, hbar.Hooov)
+
+        tmp = contract("klcd,lkdb->cb", X2_A, l2)
+        tmp = contract("jb,cb->jc", X1_B, tmp)
+        polar -= contract("jc,jc->", tmp, hbar.Hov)
+
+        tmp = contract("klcd,ljdc->kj", X2_A, l2)
+        tmp = contract("kj,jb->kb", tmp, X1_B)
+        polar -= contract("kb,kb->", tmp, hbar.Hov)
+
+        # down
+        tmp = contract('lkda,klcd->ac', l2, X2_A)
+        tmp2 = contract('jb,ajcb->ac', X1_B, hbar.Hvovv)
+        polar += 2.0 * contract('ac,ac->', tmp, tmp2)
+
+        tmp = contract('lkda,klcd->ac', l2, X2_A)
+        tmp2 = contract('jb,ajbc->ac', X1_B, hbar.Hvovv)
+        polar -= contract('ac,ac->', tmp, tmp2)
+
+        tmp = contract('jb,ljda->lbda', X1_B, l2)
+
+        # swapaxes
+        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_A, hbar.Hvovv)
+        tmp2 -= contract('klcd,akcb->ldab', X2_A, hbar.Hvovv)
+        polar += contract('lbda,ldab->', tmp, tmp2)
+
+        tmp = contract('ia,fkba->fkbi', X1_B, hbar.Hvovv)
+        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
+        polar -= contract('jkbc,kjbc->', X2_A, tmp)
+
+        tmp = contract('ia,fjac->fjic', X1_B, hbar.Hvovv)
+        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
+        polar -= contract('jkbc,jkbc->', X2_A, tmp)
+
+        tmp = contract('ia,jkfa->jkfi', X1_B, l2)
+        tmp2 = contract('jkbc,fibc->jkfi', X2_A, hbar.Hvovv)
+        polar -= contract('jkfi,jkfi->', tmp2, tmp)
+
+        tmp = contract('jb,kjib->ki', X1_B, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
+        polar -= 2.0 * contract('ki,ki->', tmp, tmp2)
+
+        tmp = contract('jb,jkib->ki', X1_B, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
+        polar += contract('ki,ki->', tmp, tmp2)
+
+        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_A)
+        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_A)
+        tmp = contract('jild,jb->bild', tmp, X1_B)
+        polar -= contract('bild,ilbd->', tmp, l2)
+
+        tmp = contract('ia,jkna->jkni', X1_B, hbar.Hooov)
+        tmp2 = contract('jkbc,nibc->jkni', X2_A, l2)
+        polar += contract('jkni,jkni->', tmp2, tmp)
+
+        tmp = contract('ia,nkab->nkib', X1_B, l2)
+        tmp = contract('jkbc,nkib->jnic', X2_A, tmp)
+        polar += contract('jnic,ijnc->', tmp, hbar.Hooov)
+
+        tmp = contract('ia,nkba->nkbi', X1_B, l2)
+        tmp = contract('jkbc,nkbi->jnci', X2_A, tmp)
+        polar += contract('jnci,jinc->', tmp, hbar.Hooov)
         # End L_{2}HX_{1}Y_{2}
 
         return polar#-1.0 * polar
diff --git a/pycc/data/molecules.py b/pycc/data/molecules.py
index 0ea22a1..d8a563c 100644
--- a/pycc/data/molecules.py
+++ b/pycc/data/molecules.py
@@ -265,6 +265,14 @@
 symmetry c1
 """
 
+h2o2 = """
+    O -0.182400 -0.692195 -0.031109
+    O 0.182400 0.692195 -0.031109
+    H 0.533952 -1.077444 0.493728
+    H -0.533952 1.077444 0.493728
+    symmetry c1
+    """
+
 moldict = {}
 moldict["He"] = he
 moldict["Be"] = be
@@ -285,3 +293,4 @@
 moldict["(S)-dimethylallene"] = sdma
 moldict["(S)-2-chloropropionitrile"] = s2cpn
 moldict["(R)-methylthiirane"] = rmt
+moldict["hydrogenperoxide"] = h2o2
diff --git a/pycc/tests/test_036_lr_sym.py b/pycc/tests/test_036_lr_sym.py
index 5084765..e04da79 100644
--- a/pycc/tests/test_036_lr_sym.py
+++ b/pycc/tests/test_036_lr_sym.py
@@ -77,10 +77,10 @@
 # Y_1 = Y(-omega); Y_2 = Y(omega)
 # X_neg = X1(-omega) , X2(-omega)
 # X_pos, Y_neg, Y_pos
-# X_1 = {}
-# X_2 = {}
-# Y_1 = {}
-# Y_2 = {}
+X_1 = {}
+X_2 = {}
+Y_1 = {}
+Y_2 = {}
 X_A = {}
 X_B = {}
 
@@ -96,14 +96,14 @@
     X_A[string] = resp.solve_right(A, omega1)
     X_B[string] = resp.solve_right(A, -omega1)
 
-    # X_2[string] = resp.solve_right(A, omega1)
-    # Y_2[string] = resp.solve_left(A, omega1)
-    # X_1[string] = resp.solve_right(A, -omega1)
-    # Y_1[string] = resp.solve_left(A, -omega1)
+    X_2[string] = resp.solve_right(A, omega1)
+    Y_2[string] = resp.solve_left(A, omega1)
+    X_1[string] = resp.solve_right(A, -omega1)
+    Y_1[string] = resp.solve_left(A, -omega1)
 
 #resp.polar(omega1)
 # Grabbing X, Y and declaring the matrix space for LR
-# polar_AB_pos = np.zeros((3,3))
+polar_AB_pos = np.zeros((3,3))
 polar_AB_neg = np.zeros((3,3))
 # polar_AB_aveg = np.zeros((3,3))
 
@@ -115,23 +115,23 @@
         # string_a = "MU_" + resp.cart[a]
         # X1_A, X2_A, _ = X_A[string_a]
         string_b = "MU_" + resp.cart[b]
-        # Y1_B, Y2_B, _ = Y_2[string_b]
-        # X1_B, X2_B, _ = X_2[string_b]
+        Y1_B, Y2_B, _ = Y_2[string_b]
+        X1_B, X2_B, _ = X_2[string_b]
         X_1B, X_2B, _ = X_B[string_b]
-        # polar_AB_pos[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
+        polar_AB_pos[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
         polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
         # polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
 
-# print(polar_AB_pos[0,0])
-# print(polar_AB_pos[0,1])
-# print(polar_AB_pos[0,2])
-# print(polar_AB_pos[1,0])
-# print(polar_AB_pos[1,1])
-# print(polar_AB_pos[1,2])
-# print(polar_AB_pos[2,0])
-# print(polar_AB_pos[2,1])
-# print(polar_AB_pos[2,2])
-# print("Comparing asymmetric versus symmetric")
+print(polar_AB_pos[0,0])
+print(polar_AB_pos[0,1])
+print(polar_AB_pos[0,2])
+print(polar_AB_pos[1,0])
+print(polar_AB_pos[1,1])
+print(polar_AB_pos[1,2])
+print(polar_AB_pos[2,0])
+print(polar_AB_pos[2,1])
+print(polar_AB_pos[2,2])
+print("Comparing asymmetric versus symmetric")
 print(polar_AB_neg[0,0])
 print(polar_AB_neg[0,1])
 print(polar_AB_neg[0,2])

From c40478c33a189d71d086dbab08a9893d95e88e7a Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Sat, 14 Sep 2024 23:08:00 +0200
Subject: [PATCH 5/7] Linear response symmetric

---
 pycc/ccresponse.py            | 380 ++++++++++++++++++++++++++++++++--
 pycc/tests/test_036_lr_sym.py |  22 +-
 2 files changed, 372 insertions(+), 30 deletions(-)

diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py
index 5700ffe..24197fd 100644
--- a/pycc/ccresponse.py
+++ b/pycc/ccresponse.py
@@ -296,6 +296,347 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis
                     X1[X_key], X2[X_key], polar = self.solve_right(self.pertbar[pertkey], -omega, e_conv, r_conv, maxiter, max_diis, start_diis)
                     check.append(polar)
 
+    def sym_linresp(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
+        contract = self.ccwfn.contract
+        o = self.ccwfn.o
+        v = self.ccwfn.v
+
+        l1 = self.cclambda.l1
+        l2 = self.cclambda.l2
+        hbar = self.hbar
+        L = self.ccwfn.H.L
+        t2 = self.ccwfn.t2
+        ERI = self.ccwfn.H.ERI
+
+
+    def sym_linresp(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
+        """
+        	Calculate the CC linear response function for polarizability at field-frequency omega(w1).
+
+        	The linear response function, <<A;B(w1)>> generally reuires the following perturbed wave functions and frequencies:
+
+        	Parameters
+        	----------
+        	pertkey_a: string
+        		String identifying the one-electron perturbation, A along a cartesian axis
+
+        	Return
+        	------
+        	polar: float
+        	     A value of the chosen linear response function corresponding to compute polariazabiltity in a specified cartesian diresction.
+        """
+
+        contract = self.ccwfn.contract
+        o = self.ccwfn.o
+        v = self.ccwfn.v
+
+        # Defining the l1 and l2
+        l1 = self.cclambda.l1
+        l2 = self.cclambda.l2
+        hbar = self.hbar
+        L = self.ccwfn.H.L
+        F = self.ccwfn.H.F
+        t2 = self.ccwfn.t2
+        ERI = self.ccwfn.H.ERI
+
+        # Please refer to eqn 45 of [Crawford: 10.1002/wcms.1406].
+        # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
+        # <<A;B>> = <0|(1 + L^(0)) { [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)] + [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] } |0>
+
+        # <0|(1 + L^(0)) [\bar{A}^(0), X^(1)(B)] |0>
+        # <0| [\bar{A}^(0), X^(1)(B)] |0> + <0| L^(0) [\bar{A}^(0), X^(1)(B)] |0>
+        #           First term                      Second term
+
+        polar1 = 0.0
+        pertbar_A = self.pertbar[pertkey_a]
+        pertbar_B = self.pertbar[pertkey_b]
+        Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
+
+
+        polar1 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
+
+        temp_ov = contract('ab, ib -> ia', pertbar_A.Avv, X1_B)
+        temp_ov -= contract('ji, ja -> ia', pertbar_A.Aoo, X1_B)
+        polar1 += contract('ia, ia', temp_ov, l1)
+
+        temp_ov = -1.0 * contract('ja, ijab -> ib', pertbar_A.Aov, X2_B)
+        temp_ov += 2.0 * contract('ja, jiab -> ib', pertbar_A.Aov, X2_B)
+        polar1 += contract("ib, ib -> ", temp_ov, l1)
+
+        # Still need to get the derivation right
+        tmp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
+        polar1 += contract("ia, ia -> ", tmp, X1_B)
+        tmp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
+        polar1 -= 0.5 * contract("ak, ka -> ", tmp, X1_B)
+        tmp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
+        polar1 -= 0.5 * contract("bk, kb -> ", tmp, X1_B)
+
+        tmp = contract("ijab, kjab -> ik", l2, X2_B)
+        polar1 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_A.Aoo)
+        tmp = contract("ijab, kiba-> jk", l2, X2_B)
+        polar1 -= 0.5 * contract("jk, kj -> ", tmp, pertbar_A.Aoo)
+        tmp = contract("ijab, ijac -> bc", l2, X2_B)
+        polar1 += 0.5 * contract("bc, bc -> ", tmp, pertbar_A.Avv)
+        tmp = contract("ijab, ijcb -> ac", l2, X2_B)
+        polar1 += 0.5 * contract("ac, ac -> ", tmp, pertbar_A.Avv)
+
+        polar1 += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
+
+        temp_ov = contract('ab, ib -> ia', pertbar_B.Avv, X1_A)
+        temp_ov -= contract('ji, ja -> ia', pertbar_B.Aoo, X1_A)
+        polar1 += contract('ia, ia', temp_ov, l1)
+
+        temp_ov = -1.0 * contract('ja, ijab -> ib', pertbar_B.Aov, X2_A)
+        temp_ov += 2.0 * contract('ja, jiab -> ib', pertbar_B.Aov, X2_A)
+        polar1 += contract("ib, ib -> ", temp_ov, l1)
+
+        # Still need to get the derivation right
+        tmp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
+        polar1 += contract("ia, ia -> ", tmp, X1_A)
+        tmp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
+        polar1 -= 0.5 * contract("ak, ka -> ", tmp, X1_A)
+        tmp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
+        polar1 -= 0.5 * contract("bk, kb -> ", tmp, X1_A)
+
+        tmp = contract("ijab, kjab -> ik", l2, X2_A)
+        polar1 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_B.Aoo)
+        tmp = contract("ijab, kiba-> jk", l2, X2_A)
+        polar1 -= 0.5 * contract("jk, kj -> ", tmp, pertbar_B.Aoo)
+        tmp = contract("ijab, ijac -> bc", l2, X2_A)
+        polar1 += 0.5 * contract("bc, bc -> ", tmp, pertbar_B.Avv)
+        tmp = contract("ijab, ijcb -> ac", l2, X2_A)
+        polar1 += 0.5 * contract("ac, ac -> ", tmp, pertbar_B.Avv)
+
+        polar1 += 2.0 * contract('ijab, ia, jb', L[o, o, v, v], X1_B, X1_A)
+
+        tmp = contract('ja, ia -> ij', hbar.Hov, X1_A)
+        tmp = contract('ij, jb -> ib', tmp, X1_B)
+        polar1 -= contract('ib, ib', tmp, l1)
+        tmp = contract('ja, ia -> ij', hbar.Hov, X1_B)
+        tmp = contract('ij, jb -> ib', tmp, X1_A)
+        polar1 -= contract('ib, ib', tmp, l1)
+
+        tmp = contract('jika, ia -> jk', 2 * hbar.Hooov, X1_B)
+        tmp = contract('jk, jb -> kb', tmp, X1_A)
+        polar1 -= contract('kb, kb', tmp, l1)
+        tmp = contract('jika, ia -> jk', hbar.Hooov.swapaxes(0, 1), X1_B)
+        tmp = contract('jk, jb -> kb', tmp, X1_A)
+        polar1 += contract('kb, kb', tmp, l1)
+        tmp = contract('jika, ia -> jk', 2 * hbar.Hooov, X1_A)
+        tmp = contract('jk, jb -> kb', tmp, X1_B)
+        polar1 -= contract('kb, kb', tmp, l1)
+        tmp = contract('jika, ia -> jk', hbar.Hooov.swapaxes(0, 1), X1_A)
+        tmp = contract('jk, jb -> kb', tmp, X1_B)
+        polar1 += contract('kb, kb', tmp, l1)
+
+        tmp = contract('cjab, jb -> ac', 2 * hbar.Hvovv, X1_A)
+        tmp = contract('ac, ia -> ic', tmp, X1_B)
+        polar1 += contract('ic, ic -> ', tmp, l1)
+        tmp = contract('cjab, jb -> ac', hbar.Hvovv.swapaxes(2, 3), X1_A)
+        tmp = contract('ac, ia -> ic', tmp, X1_B)
+        polar1 -= contract('ic, ic -> ', tmp, l1)
+        tmp = contract('cjab, jb -> ac', 2 * hbar.Hvovv, X1_B)
+        tmp = contract('ac, ia -> ic', tmp, X1_A)
+        polar1 += contract('ic, ic -> ', tmp, l1)
+        tmp = contract('cjab, jb -> ac', hbar.Hvovv.swapaxes(2, 3), X1_B)
+        tmp = contract('ac, ia -> ic', tmp, X1_A)
+        polar1 -= contract('ic, ic -> ', tmp, l1)
+
+        temp = contract("jcka, ia -> ijkc", hbar.Hovov, X1_A)
+        temp = contract("ijkc, jb -> kibc", temp, X1_B)
+        polar1 -= contract("kibc, kibc -> ", temp, l2)
+        temp = contract("jcka, ia -> ijkc", hbar.Hovov, X1_B)
+        temp = contract("ijkc, jb -> kibc", temp, X1_A)
+        polar1 -= contract("kibc, kibc -> ", temp, l2)
+        temp = contract("jcak, ia -> ijkc", hbar.Hovvo, X1_A)
+        temp = contract("ijkc, jb -> kicb", temp, X1_B)
+        polar1 -= contract("kicb, kicb -> ", temp, l2)
+        temp = contract("jcak, ia -> ijkc", hbar.Hovvo, X1_B)
+        temp = contract("ijkc, jb -> kicb", temp, X1_A)
+        polar1 -= contract("kicb, kicb -> ", temp, l2)
+        temp = contract('cdab, ia -> ibcd', hbar.Hvvvv, X1_B)
+        temp = contract('ibcd, jb -> ijcd', temp, X1_A)
+        polar1 += 0.5 * contract('ijcd, ijcd', temp, l2)
+        temp = contract('cdab, ia -> ibcd', hbar.Hvvvv, X1_A)
+        temp = contract('ibcd, jb -> ijcd', temp, X1_B)
+        polar1 += 0.5 * contract('ijcd, ijcd', temp, l2)
+        temp = contract('ijkl, ia -> klaj', hbar.Hoooo, X1_B)
+        temp = contract('klaj, jb -> klab', temp, X1_A)
+        polar1 += 0.5 * contract('klab, klab', temp, l2)
+        temp = contract('ijkl, ia -> klaj', hbar.Hoooo, X1_A)
+        temp = contract('klaj, jb -> klab', temp, X1_B)
+        polar1 += 0.5 * contract('klab, klab', temp, l2)
+        # I might re-index the three-body term for consistency
+        Goo = contract('mjab,ijab->mi', t2, l2)
+        Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
+        r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
+        r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
+        r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
+        r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
+        polar1 += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)     #Gvv
+        polar1 += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)     #Goo
+
+        temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
+        temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
+        polar1 += 2 * contract("jlbd, jlbd -> ", temp, l2)
+        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
+        temp = contract("bc, klcd -> klbd", temp, X2_B)
+        polar1 -= contract("klbd, klbd -> ", temp, l2)
+        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
+        temp = contract("jk, jlbd -> klbd", temp, X2_B)
+        polar1 -= contract("klbd, klbd -> ", temp, l2)
+        temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A)
+        temp = contract("ikab, ilad -> klbd", temp, X2_B)
+        polar1 -= contract("klbd, klbd -> ", temp, l2)
+        temp = contract("ikad, klcd -> ilac", L[o, o, v, v], X2_B)
+        temp = contract("ilac, ijab -> jlbc", temp, X2_A)
+        polar1 -= contract("jlbc, jlbc -> ", temp, l2)
+        temp = contract("ikad, ikac -> cd", L[o, o, v, v], X2_B)
+        temp = contract("cd, jlbd -> jlbc", temp, X2_A)
+        polar1 -= contract("jlbc, jlbc -> ", temp, l2)
+        temp = contract("ikad, ilad -> kl", L[o, o, v, v], X2_B)
+        temp = contract("kl, jkbc -> jlbc", temp, X2_A)
+        polar1 -= contract("jlbc, jlbc -> ", temp, l2)
+        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A)
+        temp = contract("klij, klcd -> ijcd", temp, X2_B)
+        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_B)
+        temp = contract("klij, klcd -> ijcd", temp, X2_A)
+        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A)
+        temp = contract("ilbc, jlbd -> ijcd", temp, X2_B)
+        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_B)
+        temp = contract("ilbc, jlbd -> ijcd", temp, X2_A)
+        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A)
+        temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
+        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_B)
+        temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
+        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
+        temp = contract("jc, jkbc -> kb", temp, X2_B)
+        polar1 -= contract("kb, kb", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        temp = contract("jc, jkbc -> kb", temp, X2_A)
+        polar1 -= contract("kb, kb", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
+        temp = contract("jc, jkcb -> kb", temp, X2_B)
+        polar1 += 2.0 * contract("kb, kb", temp, l1)
+        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
+        temp = contract("jc, jkcb -> kb", temp, X2_A)
+        polar1 += 2.0 * contract("kb, kb", temp, l1)
+        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_B)
+        temp = contract("jk, jb -> kb", temp, X1_A)
+        polar1 -= contract("kb, kb -> ", temp, l1)
+        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
+        temp = contract("jk, jb -> kb", temp, X1_B)
+        polar1 -= contract("kb, kb -> ", temp, l1)
+        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_B)
+        temp = contract("bc, kc -> kb", temp, X1_A)
+        polar1 -= contract("kb, kb -> ", temp, l1)
+        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
+        temp = contract("bc, kc -> kb", temp, X1_B)
+        polar1 -= contract("kb, kb -> ", temp, l1)
+
+        # <O|L2(A)[[Hbar(0),X2(B)],X1(C)]]|0>
+        tmp = contract("klcd,lkdb->cb", X2_B, l2)
+        tmp = contract("jb,cb->jc", X1_A, tmp)
+        polar1 -= contract("jc,jc->", tmp, hbar.Hov)
+        tmp = contract("klcd,ljdc->kj", X2_B, l2)
+        tmp = contract("kj,jb->kb", tmp, X1_A)
+        polar1 -= contract("kb,kb->", tmp, hbar.Hov)
+        # down
+        tmp = contract('lkda,klcd->ac', l2, X2_B)
+        tmp2 = contract('jb,ajcb->ac', X1_A, hbar.Hvovv)
+        polar1 += 2.0 * contract('ac,ac->', tmp, tmp2)
+        tmp = contract('lkda,klcd->ac', l2, X2_B)
+        tmp2 = contract('jb,ajbc->ac', X1_A, hbar.Hvovv)
+        polar1 -= contract('ac,ac->', tmp, tmp2)
+        tmp = contract('jb,ljda->lbda', X1_A, l2)
+        # swapaxes
+        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_B, hbar.Hvovv)
+        tmp2 -= contract('klcd,akcb->ldab', X2_B, hbar.Hvovv)
+        polar1 += contract('lbda,ldab->', tmp, tmp2)
+        tmp = contract('ia,fkba->fkbi', X1_A, hbar.Hvovv)
+        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
+        polar1 -= contract('jkbc,kjbc->', X2_B, tmp)
+        tmp = contract('ia,fjac->fjic', X1_A, hbar.Hvovv)
+        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
+        polar1 -= contract('jkbc,jkbc->', X2_B, tmp)
+        tmp = contract('ia,jkfa->jkfi', X1_A, l2)
+        tmp2 = contract('jkbc,fibc->jkfi', X2_B, hbar.Hvovv)
+        polar1 -= contract('jkfi,jkfi->', tmp2, tmp)
+        tmp = contract('jb,kjib->ki', X1_A, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
+        polar1 -= 2.0 * contract('ki,ki->', tmp, tmp2)
+        tmp = contract('jb,jkib->ki', X1_A, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
+        polar1 += contract('ki,ki->', tmp, tmp2)
+        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_B)
+        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_B)
+        tmp = contract('jild,jb->bild', tmp, X1_A)
+        polar1 -= contract('bild,ilbd->', tmp, l2)
+        tmp = contract('ia,jkna->jkni', X1_A, hbar.Hooov)
+        tmp2 = contract('jkbc,nibc->jkni', X2_B, l2)
+        polar1 += contract('jkni,jkni->', tmp2, tmp)
+        tmp = contract('ia,nkab->nkib', X1_A, l2)
+        tmp = contract('jkbc,nkib->jnic', X2_B, tmp)
+        polar1 += contract('jnic,ijnc->', tmp, hbar.Hooov)
+        tmp = contract('ia,nkba->nkbi', X1_A, l2)
+        tmp = contract('jkbc,nkbi->jnci', X2_B, tmp)
+        polar1 += contract('jnci,jinc->', tmp, hbar.Hooov)
+        tmp = contract("klcd,lkdb->cb", X2_A, l2)
+        tmp = contract("jb,cb->jc", X1_B, tmp)
+        polar1 -= contract("jc,jc->", tmp, hbar.Hov)
+        tmp = contract("klcd,ljdc->kj", X2_A, l2)
+        tmp = contract("kj,jb->kb", tmp, X1_B)
+        polar1 -= contract("kb,kb->", tmp, hbar.Hov)
+        # down
+        tmp = contract('lkda,klcd->ac', l2, X2_A)
+        tmp2 = contract('jb,ajcb->ac', X1_B, hbar.Hvovv)
+        polar1 += 2.0 * contract('ac,ac->', tmp, tmp2)
+        tmp = contract('lkda,klcd->ac', l2, X2_A)
+        tmp2 = contract('jb,ajbc->ac', X1_B, hbar.Hvovv)
+        polar1 -= contract('ac,ac->', tmp, tmp2)
+        tmp = contract('jb,ljda->lbda', X1_B, l2)
+        # swapaxes
+        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_A, hbar.Hvovv)
+        tmp2 -= contract('klcd,akcb->ldab', X2_A, hbar.Hvovv)
+        polar1 += contract('lbda,ldab->', tmp, tmp2)
+        tmp = contract('ia,fkba->fkbi', X1_B, hbar.Hvovv)
+        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
+        polar1 -= contract('jkbc,kjbc->', X2_A, tmp)
+        tmp = contract('ia,fjac->fjic', X1_B, hbar.Hvovv)
+        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
+        polar1 -= contract('jkbc,jkbc->', X2_A, tmp)
+        tmp = contract('ia,jkfa->jkfi', X1_B, l2)
+        tmp2 = contract('jkbc,fibc->jkfi', X2_A, hbar.Hvovv)
+        polar1 -= contract('jkfi,jkfi->', tmp2, tmp)
+        tmp = contract('jb,kjib->ki', X1_B, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
+        polar1 -= 2.0 * contract('ki,ki->', tmp, tmp2)
+        tmp = contract('jb,jkib->ki', X1_B, hbar.Hooov)
+        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
+        polar1 += contract('ki,ki->', tmp, tmp2)
+        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_A)
+        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_A)
+        tmp = contract('jild,jb->bild', tmp, X1_B)
+        polar1 -= contract('bild,ilbd->', tmp, l2)
+        tmp = contract('ia,jkna->jkni', X1_B, hbar.Hooov)
+        tmp2 = contract('jkbc,nibc->jkni', X2_A, l2)
+        polar1 += contract('jkni,jkni->', tmp2, tmp)
+        tmp = contract('ia,nkab->nkib', X1_B, l2)
+        tmp = contract('jkbc,nkib->jnic', X2_A, tmp)
+        polar1 += contract('jnic,ijnc->', tmp, hbar.Hooov)
+        tmp = contract('ia,nkba->nkbi', X1_B, l2)
+        tmp = contract('jkbc,nkbi->jnci', X2_A, tmp)
+        polar1 += contract('jnci,jinc->', tmp, hbar.Hooov)
+
+
+        return polar1
 
     def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         """
@@ -340,6 +681,7 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         pertbar_B = self.pertbar[pertkey_b]
         Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
 
+
         # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
         polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
         # Second term (contributions from the L1 and X1)
@@ -400,21 +742,21 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         temp = contract("ijab, ijcb -> ac", l2, X2_A)
         polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
 
-        # <0|(1 + L^(0))[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
-        # <0|[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0> + <0| L^(0) [[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
-        #               Fifth term                          Sixth term
-
-        # Expanding the permutational operator and implementing it explicitly
+        # # <0|(1 + L^(0))[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
+        # # <0|[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0> + <0| L^(0) [[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
+        # #               Fifth term                          Sixth term
+        #
+        # # Expanding the permutational operator and implementing it explicitly
         temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)
         polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
-        # # #
+        # # # #
         temp = contract("je, ja -> ea", X1_A, hbar.Hov)
         temp = contract("ea, ma -> me", temp, X1_B)
         polar -= contract("me, me -> ", temp, l1)
         temp = contract("je, ja -> ea", X1_B, hbar.Hov)
         temp = contract("ea, ma -> me", temp, X1_A)
         polar -= contract("me, me -> ", temp, l1)
-        #
+        # #
         temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
         temp = contract("ieam, ia -> me", temp, X1_B)
         polar -= contract("me, me -> ", temp, l1)
@@ -427,7 +769,7 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
         temp = contract("ajem, je -> ma", temp, X1_A)
         polar += contract("ma, ma -> ", temp, l1)
-        # #
+        # # #
         temp = contract("ejab, jb -> ea", (2.0 * hbar.Hvovv), X1_A)
         temp = contract("ea, ma -> me", temp, X1_B)
         polar += contract("me, me -> ", temp, l1)
@@ -440,7 +782,7 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         temp = contract("ejba, ja -> be", hbar.Hvovv.swapaxes(2, 3), X1_B)
         temp = contract("be, mb -> me", temp, X1_A)
         polar -= contract("me, me -> ", temp, l1)
-
+        #
         temp = contract("jfma, je -> amef", hbar.Hovov, X1_A)
         temp = contract("amef, na -> mnef", temp, X1_B)
         polar -= contract("mnef, mnef -> ", temp, l2)
@@ -472,7 +814,7 @@ def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_B)
         temp = contract("anij, nb -> ijab", temp, X1_A)
         polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
-
+        #
         Goo = contract('mjab,ijab->mi', t2, l2)
         Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
         r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
@@ -723,11 +1065,11 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         polar2 = 0
         pertbar_A = self.pertbar[pertkey_a]
         Avvoo = pertbar_A.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
-        # # <0|Y1(B) * A_bar|0>
-        polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
-        # <0|Y2(B) * A_bar|0>
-        polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
-        polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
+        # # # <0|Y1(B) * A_bar|0>
+        # polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
+        # # <0|Y2(B) * A_bar|0>
+        # polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
+        # polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
         # <0|[A_bar, X(B)]|0>
         polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
         # <0|L1(0) [A_bar, X1(B)]|0>
@@ -735,18 +1077,18 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         polar2 += contract("ac, ac -> ", tmp, pertbar_A.Avv)
         tmp = contract("ia, ka -> ik", l1, X1_B)
         polar2 -= contract("ik, ki -> ", tmp, pertbar_A.Aoo)
-        # <0|L1(0)[a_bar, X2(B)]|0>
+        # # <0|L1(0)[a_bar, X2(B)]|0>
         tmp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
         polar2 += 2.0 * contract("ijab, ijab -> ", tmp, X2_B)
         polar2 += -1.0 * contract("ijab, ijba -> ", tmp, X2_B)
-        # <0|L2(0)[A_bar, X1(B)]|0>
+        # # <0|L2(0)[A_bar, X1(B)]|0>
         tmp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
         polar2 += contract("ia, ia -> ", tmp, X1_B)
         tmp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
         polar2 -= 0.5 * contract("ak, ka -> ", tmp, X1_B)
         tmp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
         polar2 -= 0.5 * contract("bk, kb -> ", tmp, X1_B)
-        # <0|L2(0)[A_bar, X2(B)]|0>
+        # # <0|L2(0)[A_bar, X2(B)]|0>
         tmp = contract("ijab, kjab -> ik", l2, X2_B)
         polar2 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_A.Aoo)
         tmp = contract("ijab, kiba-> jk", l2, X2_B)
@@ -756,7 +1098,7 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         tmp = contract("ijab, ijcb -> ac", l2, X2_B)
         polar2 += 0.5 * contract("ac, ac -> ", tmp, pertbar_A.Avv)
 
-        return -1.0 * (polar1 + polar2)
+        return (polar1 + polar2)
 
 
     def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, max_diis=7, start_diis=1):
diff --git a/pycc/tests/test_036_lr_sym.py b/pycc/tests/test_036_lr_sym.py
index e04da79..516eb56 100644
--- a/pycc/tests/test_036_lr_sym.py
+++ b/pycc/tests/test_036_lr_sym.py
@@ -77,10 +77,10 @@
 # Y_1 = Y(-omega); Y_2 = Y(omega)
 # X_neg = X1(-omega) , X2(-omega)
 # X_pos, Y_neg, Y_pos
-X_1 = {}
-X_2 = {}
-Y_1 = {}
-Y_2 = {}
+# X_1 = {}
+# X_2 = {}
+# Y_1 = {}
+# Y_2 = {}
 X_A = {}
 X_B = {}
 
@@ -96,10 +96,10 @@
     X_A[string] = resp.solve_right(A, omega1)
     X_B[string] = resp.solve_right(A, -omega1)
 
-    X_2[string] = resp.solve_right(A, omega1)
-    Y_2[string] = resp.solve_left(A, omega1)
-    X_1[string] = resp.solve_right(A, -omega1)
-    Y_1[string] = resp.solve_left(A, -omega1)
+    # X_2[string] = resp.solve_right(A, omega1)
+    # Y_2[string] = resp.solve_left(A, omega1)
+    # X_1[string] = resp.solve_right(A, -omega1)
+    # Y_1[string] = resp.solve_left(A, -omega1)
 
 #resp.polar(omega1)
 # Grabbing X, Y and declaring the matrix space for LR
@@ -115,10 +115,10 @@
         # string_a = "MU_" + resp.cart[a]
         # X1_A, X2_A, _ = X_A[string_a]
         string_b = "MU_" + resp.cart[b]
-        Y1_B, Y2_B, _ = Y_2[string_b]
-        X1_B, X2_B, _ = X_2[string_b]
+        # Y1_B, Y2_B, _ = Y_2[string_b]
+        # X1_B, X2_B, _ = X_2[string_b]
         X_1B, X_2B, _ = X_B[string_b]
-        polar_AB_pos[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
+        polar_AB_pos[a,b] = resp.sym_linresp(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
         polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
         # polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
 

From 0fd447a8876b2aeb2e137447e5e63244499acea1 Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Thu, 26 Sep 2024 15:14:30 -0400
Subject: [PATCH 6/7] CCSD linear response symmetric

---
 pycc/ccresponse.py            | 836 ++++++++--------------------------
 pycc/tests/test_036_lr_sym.py | 143 ------
 pycc/tests/test_038_lr_sym.py |  78 ++++
 3 files changed, 266 insertions(+), 791 deletions(-)
 delete mode 100644 pycc/tests/test_036_lr_sym.py
 create mode 100644 pycc/tests/test_038_lr_sym.py

diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py
index 24197fd..2b2a972 100644
--- a/pycc/ccresponse.py
+++ b/pycc/ccresponse.py
@@ -296,24 +296,11 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis
                     X1[X_key], X2[X_key], polar = self.solve_right(self.pertbar[pertkey], -omega, e_conv, r_conv, maxiter, max_diis, start_diis)
                     check.append(polar)
 
-    def sym_linresp(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
-        contract = self.ccwfn.contract
-        o = self.ccwfn.o
-        v = self.ccwfn.v
-
-        l1 = self.cclambda.l1
-        l2 = self.cclambda.l2
-        hbar = self.hbar
-        L = self.ccwfn.H.L
-        t2 = self.ccwfn.t2
-        ERI = self.ccwfn.H.ERI
-
-
     def sym_linresp(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         """
-        	Calculate the CC linear response function for polarizability at field-frequency omega(w1).
+        	Calculate the CC symmetric linear response function for polarizability at field-frequency omega(w1).
 
-        	The linear response function, <<A;B(w1)>> generally reuires the following perturbed wave functions and frequencies:
+        	The linear response function, <<A;B(w1)>> generally requires the following perturbed wave functions and frequencies:
 
         	Parameters
         	----------
@@ -326,712 +313,265 @@ def sym_linresp(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
         	     A value of the chosen linear response function corresponding to compute polariazabiltity in a specified cartesian diresction.
         """
 
+        # Please refer to eqn 94 of [Koch and Jørgensen, “Coupled Cluster Response Functions.”].
+        # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
+        # <<A;B>> = <0|(1 + L^(0)) { [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)] + [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] } |0>
+
         contract = self.ccwfn.contract
+        l2 = self.cclambda.l2
+        t2 = self.ccwfn.t2
         o = self.ccwfn.o
         v = self.ccwfn.v
-
-        # Defining the l1 and l2
-        l1 = self.cclambda.l1
-        l2 = self.cclambda.l2
-        hbar = self.hbar
         L = self.ccwfn.H.L
-        F = self.ccwfn.H.F
-        t2 = self.ccwfn.t2
-        ERI = self.ccwfn.H.ERI
 
-        # Please refer to eqn 45 of [Crawford: 10.1002/wcms.1406].
-        # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
-        # <<A;B>> = <0|(1 + L^(0)) { [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)] + [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] } |0>
+        polar1 = 0.0
 
-        # <0|(1 + L^(0)) [\bar{A}^(0), X^(1)(B)] |0>
-        # <0| [\bar{A}^(0), X^(1)(B)] |0> + <0| L^(0) [\bar{A}^(0), X^(1)(B)] |0>
-        #           First term                      Second term
+        polar1 += self.LCX(pertkey_a, X1_B, X2_B)
+        polar1 += self.LCX(pertkey_b, X1_A, X2_A)
 
-        polar1 = 0.0
+        # <0|(HX1Y1)|0>
+        LHX1Y1 = 2.0 * contract('ijab, ia, jb', L[o, o, v, v], X1_B, X1_A)
+        Goo = contract('mjab,ijab->mi', t2, l2)
+        Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
+        r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
+        r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
+        r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
+        r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
+        LHX1Y1 += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)  # Gvv
+        LHX1Y1 += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)  # Goo
+        polar1 += self.LHX1Y1(X1_B, X1_A)
+        polar1 += self.LHX1Y1(X1_A, X1_B)
+
+        # <0|L2[[H, X2], Y2]|0>
+        temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
+        temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
+        LHX2Y2 = 2 * contract("jlbd, jlbd -> ", temp, l2)
+        polar1 += self.LHX2Y2(X2_A, X2_B)
+        polar1 += self.LHX2Y2(X2_B, X2_A)
+
+        polar1 += self.LHX1Y2(X1_B, X2_A)
+        polar1 += self.LHX1Y2(X1_A, X2_B)
+
+        polar1 += LHX1Y1 + LHX2Y2
+
+        return -1.0 * polar1
+
+    def LCX(self, pertkey_a, X1_B, X2_B):
+        """
+        LCX: Is the first second term contributions to the linear response symmetric function.
+        <0|(1 + L^(0)){ [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)]}|0>
+        """
+
+        contract = self.ccwfn.contract
+        l1 = self.cclambda.l1
+        l2 = self.cclambda.l2
+
+        LCX = 0.0
         pertbar_A = self.pertbar[pertkey_a]
-        pertbar_B = self.pertbar[pertkey_b]
-        Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
 
+        Aov = pertbar_A.Aov
+        Aoo = pertbar_A.Aoo
+        Avv = pertbar_A.Avv
+        Avvvo = pertbar_A.Avvvo
+        Aovoo = pertbar_A.Aovoo
 
-        polar1 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
+        # <0|[A_bar, X1]|0>
+        LCX += 2.0 * contract("ia, ia -> ", Aov, X1_B)
 
-        temp_ov = contract('ab, ib -> ia', pertbar_A.Avv, X1_B)
-        temp_ov -= contract('ji, ja -> ia', pertbar_A.Aoo, X1_B)
-        polar1 += contract('ia, ia', temp_ov, l1)
+        # <0|L1 [A_bar, X1]|0>
+        temp_ov = contract('ab, ib -> ia', Avv, X1_B)
+        temp_ov -= contract('ji, ja -> ia', Aoo, X1_B)
+        LCX += contract('ia, ia', temp_ov, l1)
 
-        temp_ov = -1.0 * contract('ja, ijab -> ib', pertbar_A.Aov, X2_B)
-        temp_ov += 2.0 * contract('ja, jiab -> ib', pertbar_A.Aov, X2_B)
-        polar1 += contract("ib, ib -> ", temp_ov, l1)
+        temp_ov = -1.0 * contract('ja, ijab -> ib', Aov, X2_B)
+        temp_ov += 2.0 * contract('ja, jiab -> ib', Aov, X2_B)
+        LCX += contract("ib, ib -> ", temp_ov, l1)
 
-        # Still need to get the derivation right
-        tmp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
-        polar1 += contract("ia, ia -> ", tmp, X1_B)
-        tmp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
-        polar1 -= 0.5 * contract("ak, ka -> ", tmp, X1_B)
-        tmp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
-        polar1 -= 0.5 * contract("bk, kb -> ", tmp, X1_B)
+        # <0|L2 [A_bar, X1]|0>
+        tmp = contract("ijbc, bcaj -> ia", l2, Avvvo)
+        LCX += contract("ia, ia -> ", tmp, X1_B)
+        tmp = contract("ijab, kbij -> ak", l2, Aovoo)
+        LCX -= 0.5 * contract("ak, ka -> ", tmp, X1_B)
+        tmp = contract("ijab, kaji -> bk", l2, Aovoo)
+        LCX -= 0.5 * contract("bk, kb -> ", tmp, X1_B)
 
         tmp = contract("ijab, kjab -> ik", l2, X2_B)
-        polar1 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_A.Aoo)
+        LCX -= 0.5 * contract("ik, ki -> ", tmp, Aoo)
         tmp = contract("ijab, kiba-> jk", l2, X2_B)
-        polar1 -= 0.5 * contract("jk, kj -> ", tmp, pertbar_A.Aoo)
+        LCX -= 0.5 * contract("jk, kj -> ", tmp, Aoo)
         tmp = contract("ijab, ijac -> bc", l2, X2_B)
-        polar1 += 0.5 * contract("bc, bc -> ", tmp, pertbar_A.Avv)
+        LCX += 0.5 * contract("bc, bc -> ", tmp, Avv)
         tmp = contract("ijab, ijcb -> ac", l2, X2_B)
-        polar1 += 0.5 * contract("ac, ac -> ", tmp, pertbar_A.Avv)
-
-        polar1 += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
+        LCX += 0.5 * contract("ac, ac -> ", tmp, Avv)
 
-        temp_ov = contract('ab, ib -> ia', pertbar_B.Avv, X1_A)
-        temp_ov -= contract('ji, ja -> ia', pertbar_B.Aoo, X1_A)
-        polar1 += contract('ia, ia', temp_ov, l1)
+        return LCX
 
-        temp_ov = -1.0 * contract('ja, ijab -> ib', pertbar_B.Aov, X2_A)
-        temp_ov += 2.0 * contract('ja, jiab -> ib', pertbar_B.Aov, X2_A)
-        polar1 += contract("ib, ib -> ", temp_ov, l1)
-
-        # Still need to get the derivation right
-        tmp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
-        polar1 += contract("ia, ia -> ", tmp, X1_A)
-        tmp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
-        polar1 -= 0.5 * contract("ak, ka -> ", tmp, X1_A)
-        tmp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
-        polar1 -= 0.5 * contract("bk, kb -> ", tmp, X1_A)
+    def LHX1Y1(self, X1_B, X1_A):
+        """
+        LHX1Y1: Is a function for the second term (quadratic term) contribution to the linear response
+        function, where both perturbed amplitudes are first order.
+        # <0|(1 + L^(0)){ [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] }|0>
+        """
 
-        tmp = contract("ijab, kjab -> ik", l2, X2_A)
-        polar1 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_B.Aoo)
-        tmp = contract("ijab, kiba-> jk", l2, X2_A)
-        polar1 -= 0.5 * contract("jk, kj -> ", tmp, pertbar_B.Aoo)
-        tmp = contract("ijab, ijac -> bc", l2, X2_A)
-        polar1 += 0.5 * contract("bc, bc -> ", tmp, pertbar_B.Avv)
-        tmp = contract("ijab, ijcb -> ac", l2, X2_A)
-        polar1 += 0.5 * contract("ac, ac -> ", tmp, pertbar_B.Avv)
+        contract = self.ccwfn.contract
+        o = self.ccwfn.o
+        v = self.ccwfn.v
+        l1 = self.cclambda.l1
+        l2 = self.cclambda.l2
+        t2 = self.ccwfn.t2
+        hbar = self.hbar
+        L = self.ccwfn.H.L
 
-        polar1 += 2.0 * contract('ijab, ia, jb', L[o, o, v, v], X1_B, X1_A)
+        LHX1Y1 = 0.0
 
+        # <0|L1[[H, X1], Y1]|0>
         tmp = contract('ja, ia -> ij', hbar.Hov, X1_A)
         tmp = contract('ij, jb -> ib', tmp, X1_B)
-        polar1 -= contract('ib, ib', tmp, l1)
-        tmp = contract('ja, ia -> ij', hbar.Hov, X1_B)
-        tmp = contract('ij, jb -> ib', tmp, X1_A)
-        polar1 -= contract('ib, ib', tmp, l1)
+        LHX1Y1 -= contract('ib, ib', tmp, l1)
 
         tmp = contract('jika, ia -> jk', 2 * hbar.Hooov, X1_B)
         tmp = contract('jk, jb -> kb', tmp, X1_A)
-        polar1 -= contract('kb, kb', tmp, l1)
+        LHX1Y1 -= contract('kb, kb', tmp, l1)
         tmp = contract('jika, ia -> jk', hbar.Hooov.swapaxes(0, 1), X1_B)
         tmp = contract('jk, jb -> kb', tmp, X1_A)
-        polar1 += contract('kb, kb', tmp, l1)
-        tmp = contract('jika, ia -> jk', 2 * hbar.Hooov, X1_A)
-        tmp = contract('jk, jb -> kb', tmp, X1_B)
-        polar1 -= contract('kb, kb', tmp, l1)
-        tmp = contract('jika, ia -> jk', hbar.Hooov.swapaxes(0, 1), X1_A)
-        tmp = contract('jk, jb -> kb', tmp, X1_B)
-        polar1 += contract('kb, kb', tmp, l1)
+        LHX1Y1 += contract('kb, kb', tmp, l1)
 
         tmp = contract('cjab, jb -> ac', 2 * hbar.Hvovv, X1_A)
         tmp = contract('ac, ia -> ic', tmp, X1_B)
-        polar1 += contract('ic, ic -> ', tmp, l1)
+        LHX1Y1 += contract('ic, ic -> ', tmp, l1)
         tmp = contract('cjab, jb -> ac', hbar.Hvovv.swapaxes(2, 3), X1_A)
         tmp = contract('ac, ia -> ic', tmp, X1_B)
-        polar1 -= contract('ic, ic -> ', tmp, l1)
-        tmp = contract('cjab, jb -> ac', 2 * hbar.Hvovv, X1_B)
-        tmp = contract('ac, ia -> ic', tmp, X1_A)
-        polar1 += contract('ic, ic -> ', tmp, l1)
-        tmp = contract('cjab, jb -> ac', hbar.Hvovv.swapaxes(2, 3), X1_B)
-        tmp = contract('ac, ia -> ic', tmp, X1_A)
-        polar1 -= contract('ic, ic -> ', tmp, l1)
+        LHX1Y1 -= contract('ic, ic -> ', tmp, l1)
 
+        # <0|L2[[H, X1], Y1]|0>
         temp = contract("jcka, ia -> ijkc", hbar.Hovov, X1_A)
         temp = contract("ijkc, jb -> kibc", temp, X1_B)
-        polar1 -= contract("kibc, kibc -> ", temp, l2)
-        temp = contract("jcka, ia -> ijkc", hbar.Hovov, X1_B)
-        temp = contract("ijkc, jb -> kibc", temp, X1_A)
-        polar1 -= contract("kibc, kibc -> ", temp, l2)
+        LHX1Y1 -= contract("kibc, kibc -> ", temp, l2)
         temp = contract("jcak, ia -> ijkc", hbar.Hovvo, X1_A)
         temp = contract("ijkc, jb -> kicb", temp, X1_B)
-        polar1 -= contract("kicb, kicb -> ", temp, l2)
-        temp = contract("jcak, ia -> ijkc", hbar.Hovvo, X1_B)
-        temp = contract("ijkc, jb -> kicb", temp, X1_A)
-        polar1 -= contract("kicb, kicb -> ", temp, l2)
+        LHX1Y1 -= contract("kicb, kicb -> ", temp, l2)
         temp = contract('cdab, ia -> ibcd', hbar.Hvvvv, X1_B)
         temp = contract('ibcd, jb -> ijcd', temp, X1_A)
-        polar1 += 0.5 * contract('ijcd, ijcd', temp, l2)
-        temp = contract('cdab, ia -> ibcd', hbar.Hvvvv, X1_A)
-        temp = contract('ibcd, jb -> ijcd', temp, X1_B)
-        polar1 += 0.5 * contract('ijcd, ijcd', temp, l2)
+        LHX1Y1 += 0.5 * contract('ijcd, ijcd', temp, l2)
         temp = contract('ijkl, ia -> klaj', hbar.Hoooo, X1_B)
         temp = contract('klaj, jb -> klab', temp, X1_A)
-        polar1 += 0.5 * contract('klab, klab', temp, l2)
-        temp = contract('ijkl, ia -> klaj', hbar.Hoooo, X1_A)
-        temp = contract('klaj, jb -> klab', temp, X1_B)
-        polar1 += 0.5 * contract('klab, klab', temp, l2)
-        # I might re-index the three-body term for consistency
-        Goo = contract('mjab,ijab->mi', t2, l2)
-        Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
-        r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
-        r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
-        r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
-        r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
-        polar1 += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)     #Gvv
-        polar1 += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)     #Goo
+        LHX1Y1 += 0.5 * contract('klab, klab', temp, l2)
 
-        temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
-        temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
-        polar1 += 2 * contract("jlbd, jlbd -> ", temp, l2)
-        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
-        temp = contract("bc, klcd -> klbd", temp, X2_B)
-        polar1 -= contract("klbd, klbd -> ", temp, l2)
-        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
-        temp = contract("jk, jlbd -> klbd", temp, X2_B)
-        polar1 -= contract("klbd, klbd -> ", temp, l2)
-        temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A)
-        temp = contract("ikab, ilad -> klbd", temp, X2_B)
-        polar1 -= contract("klbd, klbd -> ", temp, l2)
-        temp = contract("ikad, klcd -> ilac", L[o, o, v, v], X2_B)
-        temp = contract("ilac, ijab -> jlbc", temp, X2_A)
-        polar1 -= contract("jlbc, jlbc -> ", temp, l2)
-        temp = contract("ikad, ikac -> cd", L[o, o, v, v], X2_B)
-        temp = contract("cd, jlbd -> jlbc", temp, X2_A)
-        polar1 -= contract("jlbc, jlbc -> ", temp, l2)
-        temp = contract("ikad, ilad -> kl", L[o, o, v, v], X2_B)
-        temp = contract("kl, jkbc -> jlbc", temp, X2_A)
-        polar1 -= contract("jlbc, jlbc -> ", temp, l2)
-        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A)
-        temp = contract("klij, klcd -> ijcd", temp, X2_B)
-        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_B)
-        temp = contract("klij, klcd -> ijcd", temp, X2_A)
-        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A)
-        temp = contract("ilbc, jlbd -> ijcd", temp, X2_B)
-        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_B)
-        temp = contract("ilbc, jlbd -> ijcd", temp, X2_A)
-        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A)
-        temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
-        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_B)
-        temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
-        polar1 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
-        temp = contract("jc, jkbc -> kb", temp, X2_B)
-        polar1 -= contract("kb, kb", temp, l1)
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
-        temp = contract("jc, jkbc -> kb", temp, X2_A)
-        polar1 -= contract("kb, kb", temp, l1)
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
-        temp = contract("jc, jkcb -> kb", temp, X2_B)
-        polar1 += 2.0 * contract("kb, kb", temp, l1)
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
-        temp = contract("jc, jkcb -> kb", temp, X2_A)
-        polar1 += 2.0 * contract("kb, kb", temp, l1)
-        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_B)
-        temp = contract("jk, jb -> kb", temp, X1_A)
-        polar1 -= contract("kb, kb -> ", temp, l1)
-        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
-        temp = contract("jk, jb -> kb", temp, X1_B)
-        polar1 -= contract("kb, kb -> ", temp, l1)
-        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_B)
-        temp = contract("bc, kc -> kb", temp, X1_A)
-        polar1 -= contract("kb, kb -> ", temp, l1)
-        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
-        temp = contract("bc, kc -> kb", temp, X1_B)
-        polar1 -= contract("kb, kb -> ", temp, l1)
+        return LHX1Y1
 
-        # <O|L2(A)[[Hbar(0),X2(B)],X1(C)]]|0>
-        tmp = contract("klcd,lkdb->cb", X2_B, l2)
-        tmp = contract("jb,cb->jc", X1_A, tmp)
-        polar1 -= contract("jc,jc->", tmp, hbar.Hov)
-        tmp = contract("klcd,ljdc->kj", X2_B, l2)
-        tmp = contract("kj,jb->kb", tmp, X1_A)
-        polar1 -= contract("kb,kb->", tmp, hbar.Hov)
-        # down
-        tmp = contract('lkda,klcd->ac', l2, X2_B)
-        tmp2 = contract('jb,ajcb->ac', X1_A, hbar.Hvovv)
-        polar1 += 2.0 * contract('ac,ac->', tmp, tmp2)
-        tmp = contract('lkda,klcd->ac', l2, X2_B)
-        tmp2 = contract('jb,ajbc->ac', X1_A, hbar.Hvovv)
-        polar1 -= contract('ac,ac->', tmp, tmp2)
-        tmp = contract('jb,ljda->lbda', X1_A, l2)
-        # swapaxes
-        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_B, hbar.Hvovv)
-        tmp2 -= contract('klcd,akcb->ldab', X2_B, hbar.Hvovv)
-        polar1 += contract('lbda,ldab->', tmp, tmp2)
-        tmp = contract('ia,fkba->fkbi', X1_A, hbar.Hvovv)
-        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
-        polar1 -= contract('jkbc,kjbc->', X2_B, tmp)
-        tmp = contract('ia,fjac->fjic', X1_A, hbar.Hvovv)
-        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
-        polar1 -= contract('jkbc,jkbc->', X2_B, tmp)
-        tmp = contract('ia,jkfa->jkfi', X1_A, l2)
-        tmp2 = contract('jkbc,fibc->jkfi', X2_B, hbar.Hvovv)
-        polar1 -= contract('jkfi,jkfi->', tmp2, tmp)
-        tmp = contract('jb,kjib->ki', X1_A, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
-        polar1 -= 2.0 * contract('ki,ki->', tmp, tmp2)
-        tmp = contract('jb,jkib->ki', X1_A, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
-        polar1 += contract('ki,ki->', tmp, tmp2)
-        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_B)
-        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_B)
-        tmp = contract('jild,jb->bild', tmp, X1_A)
-        polar1 -= contract('bild,ilbd->', tmp, l2)
-        tmp = contract('ia,jkna->jkni', X1_A, hbar.Hooov)
-        tmp2 = contract('jkbc,nibc->jkni', X2_B, l2)
-        polar1 += contract('jkni,jkni->', tmp2, tmp)
-        tmp = contract('ia,nkab->nkib', X1_A, l2)
-        tmp = contract('jkbc,nkib->jnic', X2_B, tmp)
-        polar1 += contract('jnic,ijnc->', tmp, hbar.Hooov)
-        tmp = contract('ia,nkba->nkbi', X1_A, l2)
-        tmp = contract('jkbc,nkbi->jnci', X2_B, tmp)
-        polar1 += contract('jnci,jinc->', tmp, hbar.Hooov)
-        tmp = contract("klcd,lkdb->cb", X2_A, l2)
-        tmp = contract("jb,cb->jc", X1_B, tmp)
-        polar1 -= contract("jc,jc->", tmp, hbar.Hov)
-        tmp = contract("klcd,ljdc->kj", X2_A, l2)
-        tmp = contract("kj,jb->kb", tmp, X1_B)
-        polar1 -= contract("kb,kb->", tmp, hbar.Hov)
-        # down
-        tmp = contract('lkda,klcd->ac', l2, X2_A)
-        tmp2 = contract('jb,ajcb->ac', X1_B, hbar.Hvovv)
-        polar1 += 2.0 * contract('ac,ac->', tmp, tmp2)
-        tmp = contract('lkda,klcd->ac', l2, X2_A)
-        tmp2 = contract('jb,ajbc->ac', X1_B, hbar.Hvovv)
-        polar1 -= contract('ac,ac->', tmp, tmp2)
-        tmp = contract('jb,ljda->lbda', X1_B, l2)
-        # swapaxes
-        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_A, hbar.Hvovv)
-        tmp2 -= contract('klcd,akcb->ldab', X2_A, hbar.Hvovv)
-        polar1 += contract('lbda,ldab->', tmp, tmp2)
-        tmp = contract('ia,fkba->fkbi', X1_B, hbar.Hvovv)
-        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
-        polar1 -= contract('jkbc,kjbc->', X2_A, tmp)
-        tmp = contract('ia,fjac->fjic', X1_B, hbar.Hvovv)
-        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
-        polar1 -= contract('jkbc,jkbc->', X2_A, tmp)
-        tmp = contract('ia,jkfa->jkfi', X1_B, l2)
-        tmp2 = contract('jkbc,fibc->jkfi', X2_A, hbar.Hvovv)
-        polar1 -= contract('jkfi,jkfi->', tmp2, tmp)
-        tmp = contract('jb,kjib->ki', X1_B, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
-        polar1 -= 2.0 * contract('ki,ki->', tmp, tmp2)
-        tmp = contract('jb,jkib->ki', X1_B, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
-        polar1 += contract('ki,ki->', tmp, tmp2)
-        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_A)
-        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_A)
-        tmp = contract('jild,jb->bild', tmp, X1_B)
-        polar1 -= contract('bild,ilbd->', tmp, l2)
-        tmp = contract('ia,jkna->jkni', X1_B, hbar.Hooov)
-        tmp2 = contract('jkbc,nibc->jkni', X2_A, l2)
-        polar1 += contract('jkni,jkni->', tmp2, tmp)
-        tmp = contract('ia,nkab->nkib', X1_B, l2)
-        tmp = contract('jkbc,nkib->jnic', X2_A, tmp)
-        polar1 += contract('jnic,ijnc->', tmp, hbar.Hooov)
-        tmp = contract('ia,nkba->nkbi', X1_B, l2)
-        tmp = contract('jkbc,nkbi->jnci', X2_A, tmp)
-        polar1 += contract('jnci,jinc->', tmp, hbar.Hooov)
-
-
-        return polar1
-
-    def linresp_sym(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B):
+    def LHX2Y2(self, X2_A, X2_B):
         """
-        	Calculate the CC linear response function for polarizability at field-frequency omega(w1).
-
-        	The linear response function, <<A;B(w1)>> generally reuires the following perturbed wave functions and frequencies:
-
-        	Parameters
-        	----------
-        	pertkey_a: string
-        		String identifying the one-electron perturbation, A along a cartesian axis
-
-        	Return
-        	------
-        	polar: float
-        	     A value of the chosen linear response function corresponding to compute polariazabiltity in a specified cartesian diresction.
+        LHX2Y2: Is a function for the second term (quadratic term) contribution to the linear response
+        function, where both perturbed amplitudes are second order.
+        # <0|(1 + L^(0)){ [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] }|0>
         """
 
         contract = self.ccwfn.contract
         o = self.ccwfn.o
         v = self.ccwfn.v
-
-        # Defining the l1 and l2
         l1 = self.cclambda.l1
         l2 = self.cclambda.l2
+        t2 = self.ccwfn.t2
         hbar = self.hbar
         L = self.ccwfn.H.L
-        F = self.ccwfn.H.F
-        t2 = self.ccwfn.t2
         ERI = self.ccwfn.H.ERI
 
-        # Please refer to eqn 45 of [Crawford: 10.1002/wcms.1406].
-        # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0)
-        # <<A;B>> = <0|(1 + L^(0)) { [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)] + [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] } |0>
-
-        # <0|(1 + L^(0)) [\bar{A}^(0), X^(1)(B)] |0>
-        # <0| [\bar{A}^(0), X^(1)(B)] |0> + <0| L^(0) [\bar{A}^(0), X^(1)(B)] |0>
-        #           First term                      Second term
-
-        polar = 0.0
-        pertbar_A = self.pertbar[pertkey_a]
-        pertbar_B = self.pertbar[pertkey_b]
-        Aoovv = pertbar_B.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
-
-
-        # First term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
-        polar += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
-        # Second term (contributions from the L1 and X1)
-        temp = contract("ia, ic -> ac", l1, X1_B)
-        polar += contract("ac, ac -> ", temp, pertbar_A.Avv)
-        temp = contract("ia, ka -> ik", l1, X1_B)
-        polar -= contract("ik, ki -> ", temp, pertbar_A.Aoo)
-        # Second term (contributions from the L1 and X2)
-        temp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
-        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_B)
-        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_B)
-        # Second term (contributions from the L2 and X1)
-        temp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
-        polar += contract("ia, ia -> ", temp, X1_B)
-        temp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
-        polar -= 0.5 * contract("ak, ka -> ", temp, X1_B)
-        temp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
-        polar -= 0.5 * contract("bk, kb -> ", temp, X1_B)
-        # Second term (contributions from the L2 and X2)
-        temp = contract("ijab, kjab -> ik", l2, X2_B)
-        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_A.Aoo)
-        temp = contract("ijab, kiba-> jk", l2, X2_B)
-        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_A.Aoo)
-        temp = contract("ijab, ijac -> bc", l2, X2_B)
-        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_A.Avv)
-        temp = contract("ijab, ijcb -> ac", l2, X2_B)
-        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_A.Avv)
-
-        # <0|(1 + L^(0)) [\bar{B}^(1), X^(1)(B)] |0>
-        # <0| [\bar{B}^(1), X^(1)(B)] |0> + <0| L^(0) [\bar{B}^(0), X^(1)(B)] |0>
-        #           Third term                      Fourth term
-
-        # Third term (The contribution is only from the X1 amplitude and no contribution for the X2 amplitudes)
-        polar += 2.0 * contract("ia, ia -> ", pertbar_B.Aov, X1_A)
-        # Fourth term (contributions from the L1 and X1)
-        temp = contract("ia, ic -> ac", l1, X1_A)
-        polar += contract("ac, ac -> ", temp, pertbar_B.Avv)
-        temp = contract("ia, ka -> ik", l1, X1_A)
-        polar -= contract("ik, ki -> ", temp, pertbar_B.Aoo)
-        # Fourth term (contributions from the L1 and X2)
-        temp = contract("ia, jb -> ijab", l1, pertbar_B.Aov)
-        polar += 2.0 * contract("ijab, ijab -> ", temp, X2_A)
-        polar += -1.0 * contract("ijab, ijba -> ", temp, X2_A)
-        # Fourth term (contributions from the L2 and X1)
-        temp = contract("ijbc, bcaj -> ia", l2, pertbar_B.Avvvo)
-        polar += contract("ia, ia -> ", temp, X1_A)
-        temp = contract("ijab, kbij -> ak", l2, pertbar_B.Aovoo)
-        polar -= 0.5 * contract("ak, ka -> ", temp, X1_A)
-        temp = contract("ijab, kaji -> bk", l2, pertbar_B.Aovoo)
-        polar -= 0.5 * contract("bk, kb -> ", temp, X1_A)
-        # Fourth term (contributions from the L2 and X2)
-        temp = contract("ijab, kjab -> ik", l2, X2_A)
-        polar -= 0.5 * contract("ik, ki -> ", temp, pertbar_B.Aoo)
-        temp = contract("ijab, kiba-> jk", l2, X2_A)
-        polar -= 0.5 * contract("jk, kj -> ", temp, pertbar_B.Aoo)
-        temp = contract("ijab, ijac -> bc", l2, X2_A)
-        polar += 0.5 * contract("bc, bc -> ", temp, pertbar_B.Avv)
-        temp = contract("ijab, ijcb -> ac", l2, X2_A)
-        polar += 0.5 * contract("ac, ac -> ", temp, pertbar_B.Avv)
-
-        # # <0|(1 + L^(0))[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
-        # # <0|[[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0> + <0| L^(0) [[\bar{B}^(0), X^(1)(B)], X^(1)(B)]|0>
-        # #               Fifth term                          Sixth term
-        #
-        # # Expanding the permutational operator and implementing it explicitly
-        temp = contract("ijab, ia -> jb", L[o, o, v, v], X1_A)
-        polar += 2.0 * contract("jb, jb -> ", temp, X1_B)
-        # # # #
-        temp = contract("je, ja -> ea", X1_A, hbar.Hov)
-        temp = contract("ea, ma -> me", temp, X1_B)
-        polar -= contract("me, me -> ", temp, l1)
-        temp = contract("je, ja -> ea", X1_B, hbar.Hov)
-        temp = contract("ea, ma -> me", temp, X1_A)
-        polar -= contract("me, me -> ", temp, l1)
-        # #
-        temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_A)
-        temp = contract("ieam, ia -> me", temp, X1_B)
-        polar -= contract("me, me -> ", temp, l1)
-        temp = contract("jima, je -> ieam", (2 * hbar.Hooov), X1_B)
-        temp = contract("ieam, ia -> me", temp, X1_A)
-        polar -= contract("me, me -> ", temp, l1)
-        temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_A)
-        temp = contract("ajem, je -> ma", temp, X1_B)
-        polar += contract("ma, ma -> ", temp, l1)
-        temp = contract("jiem, ia -> ajem", hbar.Hooov.swapaxes(2, 3), X1_B)
-        temp = contract("ajem, je -> ma", temp, X1_A)
-        polar += contract("ma, ma -> ", temp, l1)
-        # # #
-        temp = contract("ejab, jb -> ea", (2.0 * hbar.Hvovv), X1_A)
-        temp = contract("ea, ma -> me", temp, X1_B)
-        polar += contract("me, me -> ", temp, l1)
-        temp = contract("ejab, jb -> ea", hbar.Hvovv.swapaxes(2, 3), X1_A)
-        temp = contract("ea, ma -> me", temp, X1_B)
-        polar -= contract("me, me -> ", temp, l1)
-        temp = contract("ejba, ja -> be", (2.0 * hbar.Hvovv), X1_B)
-        temp = contract("be, mb -> me", temp, X1_A)
-        polar += contract("me, me -> ", temp, l1)
-        temp = contract("ejba, ja -> be", hbar.Hvovv.swapaxes(2, 3), X1_B)
-        temp = contract("be, mb -> me", temp, X1_A)
-        polar -= contract("me, me -> ", temp, l1)
-        #
-        temp = contract("jfma, je -> amef", hbar.Hovov, X1_A)
-        temp = contract("amef, na -> mnef", temp, X1_B)
-        polar -= contract("mnef, mnef -> ", temp, l2)
-
-        temp = contract("jfma, je -> amef", hbar.Hovov, X1_B)
-        temp = contract("amef, na -> mnef", temp, X1_A)
-        polar -= contract("mnef, mnef -> ", temp, l2)
-
-        temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_A)
-        temp = contract("mnej, jc -> mnec", temp, X1_B)
-        polar -= contract("mnec, mnec -> ", temp, l2)
-
-        temp = contract("jeam, na -> mnej", hbar.Hovvo, X1_B)
-        temp = contract("mnej, jc -> mnec", temp, X1_A)
-        polar -= contract("mnec, mnec -> ", temp, l2)
-
-        temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_A)
-        temp = contract("abej, ie -> ijab", temp, X1_B)
-        polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
-
-        temp = contract("abef, jf -> abej", hbar.Hvvvv, X1_B)
-        temp = contract("abej, ie -> ijab", temp, X1_A)
-        polar += 0.5 * contract("ijab, ijab  -> ", temp, l2)
-
-        temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_A)
-        temp = contract("anij, nb -> ijab", temp, X1_B)
-        polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
-
-        temp = contract("mnij, ma -> anij", hbar.Hoooo, X1_B)
-        temp = contract("anij, nb -> ijab", temp, X1_A)
-        polar += 0.5 * contract("ijab, ijab -> ", temp, l2)
-        #
-        Goo = contract('mjab,ijab->mi', t2, l2)
-        Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2)
-        r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v])
-        r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v])
-        r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3)
-        r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3)
-        polar += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B)     #Gvv
-        # print("Gvv\n", polar_Gvv)
-        polar += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B)     #Goo
-        # print("Goo\n", polar_Goo)
-
-        # # Begin HX_2Y_2
-        temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A)
-        temp = contract("kjbc, klcd -> jlbd", temp, X2_B)
-        polar += 2 * contract("jlbd, jlbd -> ", temp, l2)
-
+        LHX2Y2 = 0.0
         temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
         temp = contract("bc, klcd -> klbd", temp, X2_B)
-        polar -= contract("klbd, klbd -> ", temp, l2)
+        LHX2Y2 -= contract("klbd, klbd -> ", temp, l2)
         temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
         temp = contract("jk, jlbd -> klbd", temp, X2_B)
-        polar -= contract("klbd, klbd -> ", temp, l2)
+        LHX2Y2 -= contract("klbd, klbd -> ", temp, l2)
         temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A)
         temp = contract("ikab, ilad -> klbd", temp, X2_B)
-        polar -= contract("klbd, klbd -> ", temp, l2)
-
-        temp = contract("ikad, klcd -> ilac", L[o, o, v, v], X2_B)
-        temp = contract("ilac, ijab -> jlbc", temp, X2_A)
-        polar -= contract("jlbc, jlbc -> ", temp, l2)
-        temp = contract("ikad, ikac -> cd", L[o, o, v, v], X2_B)
-        temp = contract("cd, jlbd -> jlbc", temp, X2_A)
-        polar -= contract("jlbc, jlbc -> ", temp, l2)
-        temp = contract("ikad, ilad -> kl", L[o, o, v, v], X2_B)
-        temp = contract("kl, jkbc -> jlbc", temp, X2_A)
-        polar -= contract("jlbc, jlbc -> ", temp, l2)
-
-
+        LHX2Y2 -= contract("klbd, klbd -> ", temp, l2)
         temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A)
         temp = contract("klij, klcd -> ijcd", temp, X2_B)
-        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-
-        temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_B)
-        temp = contract("klij, klcd -> ijcd", temp, X2_A)
-        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-
+        LHX2Y2 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
         temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A)
         temp = contract("ilbc, jlbd -> ijcd", temp, X2_B)
-        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-
-        temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_B)
-        temp = contract("ilbc, jlbd -> ijcd", temp, X2_A)
-        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-
+        LHX2Y2 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
         temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A)
         temp = contract("ikbd, jkbc -> ijcd", temp, X2_B)
-        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
+        LHX2Y2 += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
 
-        temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_B)
-        temp = contract("ikbd, jkbc -> ijcd", temp, X2_A)
-        polar += 0.5 * contract("ijcd, ijcd -> ", temp, l2)
-        # # End HX_{2}Y_{2}
+        return LHX2Y2
 
-        # Begin L_{1}HX_{1}Y_{2}
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
-        temp = contract("jc, jkbc -> kb", temp, X2_B)
-        polar -= contract("kb, kb", temp, l1)
+    def LHX1Y2(self, X1_B, X2_A):
+        """
+        LHX1Y2: Is a function for the second term (quadratic term) contribution to the linear response
+        function, where perturbed amplitudes is a first order and a second order.
+        # <0|(1 + L^(0)){ [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] }|0>
+        """
+
+        contract = self.ccwfn.contract
+        o = self.ccwfn.o
+        v = self.ccwfn.v
+        l1 = self.cclambda.l1
+        l2 = self.cclambda.l2
+        hbar = self.hbar
+        L = self.ccwfn.H.L
+
+        LHX1Y2 = 0.0
+        # <O|L1(0)[[Hbar(0),X1(B)],X2(A)]]|0>
         temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
         temp = contract("jc, jkbc -> kb", temp, X2_A)
-        polar -= contract("kb, kb", temp, l1)
-        temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_A)
-        temp = contract("jc, jkcb -> kb", temp, X2_B)
-        polar += 2.0 * contract("kb, kb", temp, l1)
+        LHX1Y2 -= contract("kb, kb", temp, l1)
         temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B)
         temp = contract("jc, jkcb -> kb", temp, X2_A)
-        polar += 2.0 * contract("kb, kb", temp, l1)
-
-        temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_B)
-        temp = contract("jk, jb -> kb", temp, X1_A)
-        polar -= contract("kb, kb -> ", temp, l1)
+        LHX1Y2 += 2.0 * contract("kb, kb", temp, l1)
         temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A)
         temp = contract("jk, jb -> kb", temp, X1_B)
-        polar -= contract("kb, kb -> ", temp, l1)
-
-        temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_B)
-        temp = contract("bc, kc -> kb", temp, X1_A)
-        polar -= contract("kb, kb -> ", temp, l1)
+        LHX1Y2 -= contract("kb, kb -> ", temp, l1)
         temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A)
         temp = contract("bc, kc -> kb", temp, X1_B)
-        polar -= contract("kb, kb -> ", temp, l1)
-        # End L_{1}HX_{1}Y_{2}
+        LHX1Y2 -= contract("kb, kb -> ", temp, l1)
 
-        # Begin L_{2}HX_{1}Y_{2}
         # <O|L2(A)[[Hbar(0),X2(B)],X1(C)]]|0>
-        tmp = contract("klcd,lkdb->cb", X2_B, l2)
-        tmp = contract("jb,cb->jc", X1_A, tmp)
-        polar -= contract("jc,jc->", tmp, hbar.Hov)
-
-        tmp = contract("klcd,ljdc->kj", X2_B, l2)
-        tmp = contract("kj,jb->kb", tmp, X1_A)
-        polar -= contract("kb,kb->", tmp, hbar.Hov)
-
-        # down
-        tmp = contract('lkda,klcd->ac', l2, X2_B)
-        tmp2 = contract('jb,ajcb->ac', X1_A, hbar.Hvovv)
-        polar += 2.0 * contract('ac,ac->', tmp, tmp2)
-
-        tmp = contract('lkda,klcd->ac', l2, X2_B)
-        tmp2 = contract('jb,ajbc->ac', X1_A, hbar.Hvovv)
-        polar -= contract('ac,ac->', tmp, tmp2)
-
-        tmp = contract('jb,ljda->lbda', X1_A, l2)
-
-        # swapaxes
-        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_B, hbar.Hvovv)
-        tmp2 -= contract('klcd,akcb->ldab', X2_B, hbar.Hvovv)
-        polar += contract('lbda,ldab->', tmp, tmp2)
-
-        tmp = contract('ia,fkba->fkbi', X1_A, hbar.Hvovv)
-        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
-        polar -= contract('jkbc,kjbc->', X2_B, tmp)
-
-        tmp = contract('ia,fjac->fjic', X1_A, hbar.Hvovv)
-        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
-        polar -= contract('jkbc,jkbc->', X2_B, tmp)
-
-        tmp = contract('ia,jkfa->jkfi', X1_A, l2)
-        tmp2 = contract('jkbc,fibc->jkfi', X2_B, hbar.Hvovv)
-        polar -= contract('jkfi,jkfi->', tmp2, tmp)
-
-        tmp = contract('jb,kjib->ki', X1_A, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
-        polar -= 2.0 * contract('ki,ki->', tmp, tmp2)
-
-        tmp = contract('jb,jkib->ki', X1_A, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_B, l2)
-        polar += contract('ki,ki->', tmp, tmp2)
-
-        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_B)
-        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_B)
-        tmp = contract('jild,jb->bild', tmp, X1_A)
-        polar -= contract('bild,ilbd->', tmp, l2)
-
-        tmp = contract('ia,jkna->jkni', X1_A, hbar.Hooov)
-        tmp2 = contract('jkbc,nibc->jkni', X2_B, l2)
-        polar += contract('jkni,jkni->', tmp2, tmp)
-
-        tmp = contract('ia,nkab->nkib', X1_A, l2)
-        tmp = contract('jkbc,nkib->jnic', X2_B, tmp)
-        polar += contract('jnic,ijnc->', tmp, hbar.Hooov)
-
-        tmp = contract('ia,nkba->nkbi', X1_A, l2)
-        tmp = contract('jkbc,nkbi->jnci', X2_B, tmp)
-        polar += contract('jnci,jinc->', tmp, hbar.Hooov)
-
-        tmp = contract("klcd,lkdb->cb", X2_A, l2)
-        tmp = contract("jb,cb->jc", X1_B, tmp)
-        polar -= contract("jc,jc->", tmp, hbar.Hov)
-
-        tmp = contract("klcd,ljdc->kj", X2_A, l2)
-        tmp = contract("kj,jb->kb", tmp, X1_B)
-        polar -= contract("kb,kb->", tmp, hbar.Hov)
-
-        # down
-        tmp = contract('lkda,klcd->ac', l2, X2_A)
-        tmp2 = contract('jb,ajcb->ac', X1_B, hbar.Hvovv)
-        polar += 2.0 * contract('ac,ac->', tmp, tmp2)
-
-        tmp = contract('lkda,klcd->ac', l2, X2_A)
-        tmp2 = contract('jb,ajbc->ac', X1_B, hbar.Hvovv)
-        polar -= contract('ac,ac->', tmp, tmp2)
-
-        tmp = contract('jb,ljda->lbda', X1_B, l2)
-
-        # swapaxes
-        tmp2 = 2.0 * contract('klcd,akbc->ldab', X2_A, hbar.Hvovv)
-        tmp2 -= contract('klcd,akcb->ldab', X2_A, hbar.Hvovv)
-        polar += contract('lbda,ldab->', tmp, tmp2)
-
-        tmp = contract('ia,fkba->fkbi', X1_B, hbar.Hvovv)
-        tmp = contract('fkbi,jifc->kjbc', tmp, l2)
-        polar -= contract('jkbc,kjbc->', X2_A, tmp)
-
-        tmp = contract('ia,fjac->fjic', X1_B, hbar.Hvovv)
-        tmp = contract('fjic,ikfb->jkbc', tmp, l2)
-        polar -= contract('jkbc,jkbc->', X2_A, tmp)
-
-        tmp = contract('ia,jkfa->jkfi', X1_B, l2)
-        tmp2 = contract('jkbc,fibc->jkfi', X2_A, hbar.Hvovv)
-        polar -= contract('jkfi,jkfi->', tmp2, tmp)
-
-        tmp = contract('jb,kjib->ki', X1_B, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
-        polar -= 2.0 * contract('ki,ki->', tmp, tmp2)
-
-        tmp = contract('jb,jkib->ki', X1_B, hbar.Hooov)
-        tmp2 = contract('klcd,ilcd->ki', X2_A, l2)
-        polar += contract('ki,ki->', tmp, tmp2)
-
-        tmp = 2.0 * contract('jkic,klcd->jild', hbar.Hooov, X2_A)
-        tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_A)
-        tmp = contract('jild,jb->bild', tmp, X1_B)
-        polar -= contract('bild,ilbd->', tmp, l2)
-
-        tmp = contract('ia,jkna->jkni', X1_B, hbar.Hooov)
-        tmp2 = contract('jkbc,nibc->jkni', X2_A, l2)
-        polar += contract('jkni,jkni->', tmp2, tmp)
-
-        tmp = contract('ia,nkab->nkib', X1_B, l2)
-        tmp = contract('jkbc,nkib->jnic', X2_A, tmp)
-        polar += contract('jnic,ijnc->', tmp, hbar.Hooov)
-
-        tmp = contract('ia,nkba->nkbi', X1_B, l2)
-        tmp = contract('jkbc,nkbi->jnci', X2_A, tmp)
-        polar += contract('jnci,jinc->', tmp, hbar.Hooov)
-        # End L_{2}HX_{1}Y_{2}
-
-        return polar#-1.0 * polar
+        temp = contract('ja, ia -> ij', hbar.Hov, X1_B)
+        temp = contract('ij, jkbc -> ikbc', temp, X2_A)
+        LHX1Y2 -= contract('ikbc, ikbc', temp, l2)
+        temp = contract('ja, jb -> ab', hbar.Hov, X1_B)
+        temp = contract('ab, ikac -> ikbc', temp, X2_A)
+        LHX1Y2 -= contract('ikbc, ikbc', temp, l2)
+
+        temp = contract('jima, ia -> jm', (2*hbar.Hooov - hbar.Hooov.swapaxes(0, 1)), X1_B)
+        temp = contract('jm, jkbc -> kmcb', temp, X2_A)
+        LHX1Y2 -= contract('kmcb, kmcb -> ', temp, l2)
+        temp = contract('jima, jb -> imab', (2*hbar.Hooov - hbar.Hooov.swapaxes(0, 1)), X1_B)
+        temp = contract('imab, ikac -> kmcb', temp, X2_A)
+        LHX1Y2 -= contract('kmcb, kmcb', temp, l2)
+
+        temp = contract('djab, ia -> djib', (2*hbar.Hvovv - hbar.Hvovv.swapaxes(2, 3)), X1_B)
+        temp = contract('djib, jkbc -> kicd', temp, X2_A)
+        LHX1Y2 += contract('kicd, kicd', temp, l2)
+        temp = contract('djab, jb -> ad', (2 * hbar.Hvovv - hbar.Hvovv.swapaxes(2, 3)), X1_B)
+        temp = contract('ad, ikac -> kicd', temp, X2_A)
+        LHX1Y2 += contract('kicd, kicd', temp, l2)
+
+        temp = contract('kdab, ia -> ikdb', hbar.Hvovv.swapaxes(0, 1), X1_B)
+        temp = contract('ikdb, jkbc -> ijdc', temp, X2_A)
+        LHX1Y2 -= contract('ijdc, ijdc', temp, l2)
+        temp = contract('kdab, jb -> kjad', hbar.Hvovv.swapaxes(0, 1), X1_B)
+        temp = contract('kjad, ikac -> ijdc', temp, X2_A)
+        LHX1Y2 -= contract('ijdc, ijdc', temp, l2)
+        temp = contract('kdab, kc -> abcd', hbar.Hvovv.swapaxes(0, 1), X1_B)
+        temp = contract('abcd, ijab -> ijdc', temp, X2_A)
+        LHX1Y2 -= contract('ijdc, ijdc', temp, l2)
+
+        temp = contract('jkla, ia -> ijkl', hbar.Hooov, X1_B)
+        temp = contract('ijkl, jkbc -> libc', temp, X2_A)
+        LHX1Y2 += contract('libc, libc', temp, l2)
+        temp = contract('jkla, jb -> klab', hbar.Hooov, X1_B)
+        temp = contract('klab, ikac -> libc', temp, X2_A)
+        LHX1Y2 += contract('libc, libc -> ', temp, l2)
+        temp = contract('jkla, kc -> jlac', hbar.Hooov, X1_B)
+        temp = contract('jlac, ijab -> libc', temp, X2_A)
+        LHX1Y2 += contract('libc, libc -> ', temp, l2)
+
+
+        return LHX1Y2
 
     def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         """
@@ -1065,11 +605,11 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         polar2 = 0
         pertbar_A = self.pertbar[pertkey_a]
         Avvoo = pertbar_A.Avvoo.swapaxes(0, 2).swapaxes(1, 3)
-        # # # <0|Y1(B) * A_bar|0>
-        # polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
-        # # <0|Y2(B) * A_bar|0>
-        # polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
-        # polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
+        # <0|Y1(B) * A_bar|0>
+        polar1 += contract("ai, ia -> ", pertbar_A.Avo, Y1_B)
+        # <0|Y2(B) * A_bar|0>
+        polar1 += 0.5 * contract("abij, ijab -> ", Avvoo, Y2_B)
+        polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B)
         # <0|[A_bar, X(B)]|0>
         polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B)
         # <0|L1(0) [A_bar, X1(B)]|0>
@@ -1077,18 +617,18 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         polar2 += contract("ac, ac -> ", tmp, pertbar_A.Avv)
         tmp = contract("ia, ka -> ik", l1, X1_B)
         polar2 -= contract("ik, ki -> ", tmp, pertbar_A.Aoo)
-        # # <0|L1(0)[a_bar, X2(B)]|0>
+        # <0|L1(0)[a_bar, X2(B)]|0>
         tmp = contract("ia, jb -> ijab", l1, pertbar_A.Aov)
         polar2 += 2.0 * contract("ijab, ijab -> ", tmp, X2_B)
         polar2 += -1.0 * contract("ijab, ijba -> ", tmp, X2_B)
-        # # <0|L2(0)[A_bar, X1(B)]|0>
+        # <0|L2(0)[A_bar, X1(B)]|0>
         tmp = contract("ijbc, bcaj -> ia", l2, pertbar_A.Avvvo)
         polar2 += contract("ia, ia -> ", tmp, X1_B)
         tmp = contract("ijab, kbij -> ak", l2, pertbar_A.Aovoo)
         polar2 -= 0.5 * contract("ak, ka -> ", tmp, X1_B)
         tmp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo)
         polar2 -= 0.5 * contract("bk, kb -> ", tmp, X1_B)
-        # # <0|L2(0)[A_bar, X2(B)]|0>
+        # <0|L2(0)[A_bar, X2(B)]|0>
         tmp = contract("ijab, kjab -> ik", l2, X2_B)
         polar2 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_A.Aoo)
         tmp = contract("ijab, kiba-> jk", l2, X2_B)
@@ -1098,7 +638,7 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B):
         tmp = contract("ijab, ijcb -> ac", l2, X2_B)
         polar2 += 0.5 * contract("ac, ac -> ", tmp, pertbar_A.Avv)
 
-        return (polar1 + polar2)
+        return -1.0 * (polar1 + polar2)
 
 
     def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, max_diis=7, start_diis=1):
diff --git a/pycc/tests/test_036_lr_sym.py b/pycc/tests/test_036_lr_sym.py
deleted file mode 100644
index 516eb56..0000000
--- a/pycc/tests/test_036_lr_sym.py
+++ /dev/null
@@ -1,143 +0,0 @@
-"""
-Test CCSD linear response functions.
-"""
-import sys
-sys.path.append("/Users/jattakumi/pycc/")
-import numpy as np
-# Import package, test suite, and other packages as needed
-import psi4
-from opt_einsum import contract
-
-import pycc
-#from ccwfn import ccwfn
-#from cchbar import cchbar
-#from cclambda import cclambda
-#from ccdensity import ccdensity
-#from ccresponse import ccresponse
-
-geom = """
-O
-H 1 1.8084679
-H 1 1.8084679 2 104.5
-units bohr
-symmetry c1 
-no_reorient
-"""
-
-hf = """
-F  0.000000000000   0.000000000000  0.000000000000
-H  0.000000000000   0.000000000000  -1.732800000000
-units bohr
-no_reorient
-symmetry c1
-"""
-
-hof = """
-          O          -0.947809457408    -0.132934425181     0.000000000000
-          H          -1.513924046286     1.610489987673     0.000000000000
-          F           0.878279174340     0.026485523618     0.000000000000
-unit bohr
-no_reorient
-symmetry c1
-"""
-psi4.core.clean_options()
-psi4.set_memory('2 GiB')
-psi4.core.set_output_file('output.dat', False)
-psi4.set_options({'basis': 'cc-pVDZ',
-                  'scf_type': 'pk',
-                  'mp2_type': 'conv',
-                  'freeze_core': 'true',
-                  'e_convergence': 1e-12,
-                  'd_convergence': 1e-12,
-                  'r_convergence': 1e-12})
-mol = psi4.geometry(hof)
-rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
-
-e_conv = 1e-12
-r_conv = 1e-12
-
-cc = pycc.ccwfn(rhf_wfn)
-ecc = cc.solve_cc(e_conv, r_conv)
-hbar = pycc.cchbar(cc)
-cclambda = pycc.cclambda(cc, hbar)
-lecc = cclambda.solve_lambda(e_conv, r_conv)
-density = pycc.ccdensity(cc, cclambda)
-
-resp = pycc.ccresponse(density)
-
-omega1 = 0.0
-# omega1 = 0.0656
-
-
-#A = resp.pertbar.Aoo
-
-#resp.linresp(omega1)
-# Creating dictionaries
-# X_1 = X(-omega); X_2 = X(omega)
-# Y_1 = Y(-omega); Y_2 = Y(omega)
-# X_neg = X1(-omega) , X2(-omega)
-# X_pos, Y_neg, Y_pos
-# X_1 = {}
-# X_2 = {}
-# Y_1 = {}
-# Y_2 = {}
-X_A = {}
-X_B = {}
-
-for axis in range(0, 3):
-    string = "MU_" + resp.cart[axis]
-
-    A = resp.pertbar[string]
-    # B = resp.pertbar[string]
-    #print("Aoo",A)
-
-    # A -> -omega
-    # B -> +omega
-    X_A[string] = resp.solve_right(A, omega1)
-    X_B[string] = resp.solve_right(A, -omega1)
-
-    # X_2[string] = resp.solve_right(A, omega1)
-    # Y_2[string] = resp.solve_left(A, omega1)
-    # X_1[string] = resp.solve_right(A, -omega1)
-    # Y_1[string] = resp.solve_left(A, -omega1)
-
-#resp.polar(omega1)
-# Grabbing X, Y and declaring the matrix space for LR
-polar_AB_pos = np.zeros((3,3))
-polar_AB_neg = np.zeros((3,3))
-# polar_AB_aveg = np.zeros((3,3))
-
-
-for a in range(0, 3):
-    string_a = "MU_" + resp.cart[a]
-    X1_A, X2_A, _ = X_A[string_a]
-    for b in range(0, 3):
-        # string_a = "MU_" + resp.cart[a]
-        # X1_A, X2_A, _ = X_A[string_a]
-        string_b = "MU_" + resp.cart[b]
-        # Y1_B, Y2_B, _ = Y_2[string_b]
-        # X1_B, X2_B, _ = X_2[string_b]
-        X_1B, X_2B, _ = X_B[string_b]
-        polar_AB_pos[a,b] = resp.sym_linresp(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
-        polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
-        # polar_AB_neg[a, b] = resp.linresp_sym(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
-
-print(polar_AB_pos[0,0])
-print(polar_AB_pos[0,1])
-print(polar_AB_pos[0,2])
-print(polar_AB_pos[1,0])
-print(polar_AB_pos[1,1])
-print(polar_AB_pos[1,2])
-print(polar_AB_pos[2,0])
-print(polar_AB_pos[2,1])
-print(polar_AB_pos[2,2])
-print("Comparing asymmetric versus symmetric")
-print(polar_AB_neg[0,0])
-print(polar_AB_neg[0,1])
-print(polar_AB_neg[0,2])
-print(polar_AB_neg[1,0])
-print(polar_AB_neg[1,1])
-print(polar_AB_neg[1,2])
-print(polar_AB_neg[2,0])
-print(polar_AB_neg[2,1])
-print(polar_AB_neg[2,2])
diff --git a/pycc/tests/test_038_lr_sym.py b/pycc/tests/test_038_lr_sym.py
new file mode 100644
index 0000000..cf82784
--- /dev/null
+++ b/pycc/tests/test_038_lr_sym.py
@@ -0,0 +1,78 @@
+"""
+Test CCSD linear response functions.
+"""
+
+# Import package, test suite, and other packages as needed
+import psi4
+import numpy as np
+import pytest
+import pycc
+from ..data.molecules import *
+
+def test_sym_linresp():
+
+    psi4.core.clean_options()
+    psi4.set_memory('2 GiB')
+    psi4.core.set_output_file('output.dat', False)
+    psi4.set_options({'basis': 'aug-cc-pVDZ',
+                      'scf_type': 'pk',
+                      'mp2_type': 'conv',
+                      'freeze_core': 'true',
+                      'e_convergence': 1e-12,
+                      'd_convergence': 1e-12,
+                      'r_convergence': 1e-12})
+
+    mol = psi4.geometry(moldict["H2O"])
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    e_conv = 1e-12
+    r_conv = 1e-12
+
+    cc = pycc.ccwfn(rhf_wfn)
+    ecc = cc.solve_cc(e_conv, r_conv)
+    hbar = pycc.cchbar(cc)
+    cclambda = pycc.cclambda(cc, hbar)
+    lecc = cclambda.solve_lambda(e_conv, r_conv)
+    density = pycc.ccdensity(cc, cclambda)
+
+    resp = pycc.ccresponse(density)
+
+    omega1 = 0.0656
+
+    # Creating dictionaries
+    X_A = {}
+    X_B = {}
+
+    for axis in range(0, 3):
+        string = "MU_" + resp.cart[axis]
+        A = resp.pertbar[string]
+        X_A[string] = resp.solve_right(A, omega1)
+        X_B[string] = resp.solve_right(A, -omega1)
+
+    # Grabbing X, Y and declaring the matrix space for LR
+    polar_AB = np.zeros((3,3))
+
+    for a in range(0, 3):
+        string_a = "MU_" + resp.cart[a]
+        X1_A, X2_A, _ = X_A[string_a]
+        for b in range(0, 3):
+            string_b = "MU_" + resp.cart[b]
+            X_1B, X_2B, _ = X_B[string_b]
+            polar_AB[a,b] = resp.sym_linresp(string_a, string_b, X1_A, X2_A, X_1B, X_2B)
+
+    print(f"Dynamic Polarizability Tensor @ w = {omega1} a.u.:")
+    print(polar_AB)
+    print("Average Dynamic Polarizability:")
+    polar_AB_avg = np.average([polar_AB[0,0], polar_AB[1,1], polar_AB[2,2]])
+    print(polar_AB_avg)
+
+    # Validating from psi4
+    polar_xx = 9.932240347101651
+    polar_yy = 13.446487681337629
+    polar_zz = 11.344346098120035
+    polar_avg = 11.574358042186
+
+    assert(abs(polar_AB[0,0] - polar_xx) < 1e-7)
+    assert(abs(polar_AB[1,1] - polar_yy) < 1e-7)
+    assert(abs(polar_AB[2,2] - polar_zz) < 1e-7)
+    assert(abs(polar_AB_avg - polar_avg) < 1e-7)

From 601cd3c1d5649977642b3d293fa8e8cbbf76342c Mon Sep 17 00:00:00 2001
From: jattakumi <jattakumi123@vt.edu>
Date: Thu, 26 Sep 2024 15:33:11 -0400
Subject: [PATCH 7/7] Fixing conflicts

---
 pycc/tests/test_036_lr.py | 176 ++++++++++++++------------------------
 1 file changed, 62 insertions(+), 114 deletions(-)

diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py
index a555240..17c61cd 100644
--- a/pycc/tests/test_036_lr.py
+++ b/pycc/tests/test_036_lr.py
@@ -5,56 +5,11 @@
 # Import package, test suite, and other packages as needed
 import psi4
 import numpy as np
-import sys
-sys.path.append("/Users/jattakumi/pycc/pycc/")
 import pytest
 import pycc
-from data.molecules import *
+from ..data.molecules import *
 
 def test_linresp():
-    h2o2 = """
-    O -0.182400 -0.692195 -0.031109
-    O 0.182400 0.692195 -0.031109
-    H 0.533952 -1.077444 0.493728
-    H -0.533952 1.077444 0.493728
-    symmetry c1
-    """
-
-    # H2_4 chain
-    h2_4 = """
-    H 0.000000 0.000000 0.000000
-    H 0.750000 0.000000 0.000000
-    H 0.000000 1.500000 0.000000
-    H 0.375000 1.500000 -0.649520
-    H 0.000000 3.000000 0.000000
-    H -0.375000 3.000000 -0.649520
-    H 0.000000 4.500000 -0.000000
-    H -0.750000 4.500000 -0.000000
-    symmetry c1
-    """
-
-    # H2_7 chain
-    h2_7 = """
-    H 0.000000 0.000000 0.000000
-    H 0.750000 0.000000 0.000000
-    H 0.000000 1.500000 0.000000
-    H 0.375000 1.500000 -0.649520
-    H 0.000000 3.000000 0.000000
-    H -0.375000 3.000000 -0.649520
-    H 0.000000 4.500000 -0.000000
-    H -0.750000 4.500000 -0.000000
-    H 0.000000 6.000000 -0.000000
-    H -0.375000 6.000000 0.649520
-    H 0.000000 7.500000 -0.000000
-    H 0.375000 7.500000 -0.649520
-    H 0.000000 9.000000 -0.000000
-    H 0.750000 9.000000 0.000000
-    symmetry c1
-    """
-
-    threshold = [1e-03, 1e-04, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10, 1e-11, 1e-12, 0]
-    lmo = 'PNO++'
-    geom = "h2_4"
     psi4.core.clean_options()
     psi4.set_memory('2 GiB')
     psi4.set_output_file('output.dat', False)
@@ -67,73 +22,66 @@ def test_linresp():
                       'r_convergence': 1e-12,
     })
 
-    mol = psi4.geometry(h2_4)
+    mol = psi4.geometry(moldict["H2O"])
     rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
 
-    for t in threshold:
-        e_conv = 1e-12
-        r_conv = 1e-12
-
-        cc = pycc.ccwfn(rhf_wfn, local_mos = 'BOYS', local= lmo, local_cutoff = t, filter=True)
-        ecc = cc.solve_cc(e_conv, r_conv)
-        hbar = pycc.cchbar(cc)
-        cclambda = pycc.cclambda(cc, hbar)
-        lecc = cclambda.solve_lambda(e_conv, r_conv)
-        density = pycc.ccdensity(cc, cclambda)
-
-        resp = pycc.ccresponse(density)
-
-        omega1 = 0.0656
-
-        # Creating dictionaries
-        # X_1 = X(-omega); X_2 = X(omega)
-        # Y_1 = Y(-omega); Y_2 = Y(omega)
-        X_1 = {}
-        X_2 = {}
-        Y_1 = {}
-        Y_2 = {}
-
-        for axis in range(0, 3):
-            string = "MU_" + resp.cart[axis]
-
-            A = resp.pertbar[string]
-
-            X_2[string] = resp.solve_right(A, omega1)
-            Y_2[string] = resp.solve_left(A, omega1)
-            X_1[string] = resp.solve_right(A, -omega1)
-            Y_1[string] = resp.solve_left(A, -omega1)
-
-        # Grabbing X, Y and declaring the matrix space for LR
-        polar_AB = np.zeros((3,3))
-
-        for a in range(0, 3):
-            string_a = "MU_" + resp.cart[a]
-            for b in range(0,3):
-                string_b = "MU_" + resp.cart[b]
-                Y1_B, Y2_B, _ = Y_2[string_b]
-                X1_B, X2_B, _ = X_2[string_b]
-                polar_AB[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
-
-        print("Dynamic Polarizability Tensor @ w=0.0656 a.u.:")
-        print(polar_AB)
-        print("Average Dynamic Polarizability:")
-        polar_AB_avg = np.average([polar_AB[0,0], polar_AB[1,1], polar_AB[2,2]])
-        print(polar_AB_avg)
-
-        B_avg = str(geom) + "_Dynpolar_" + str(lmo) + ".txt"
-
-        with open(B_avg, 'a') as f1:
-            f1.write(str(polar_AB_avg) + ' ' + str(t) + '\n')
-
-        del cc, hbar, cclambda, density, resp
-
-    # #validating from psi4
-    # polar_XX = 9.92992070420665
-    # polar_YY = 13.443740151331559
-    # polar_ZZ = 11.342765745046526
-    # polar_avg = 11.572142200333
-    #
-    # assert(abs(polar_AB[0,0] - polar_XX) < 1e-8)
-    # assert(abs(polar_AB[1,1] - polar_YY) < 1e-8)
-    # assert(abs(polar_AB[2,2] - polar_ZZ) < 1e-8)
-    # assert(abs(polar_AB_avg - polar_avg) < 1e-8)
+    # for t in threshold:
+    e_conv = 1e-12
+    r_conv = 1e-12
+
+    cc = pycc.ccwfn(rhf_wfn)
+    ecc = cc.solve_cc(e_conv, r_conv)
+    hbar = pycc.cchbar(cc)
+    cclambda = pycc.cclambda(cc, hbar)
+    lecc = cclambda.solve_lambda(e_conv, r_conv)
+    density = pycc.ccdensity(cc, cclambda)
+
+    resp = pycc.ccresponse(density)
+
+    omega1 = 0.0656
+
+    # Creating dictionaries
+    # X_1 = X(-omega); X_2 = X(omega)
+    # Y_1 = Y(-omega); Y_2 = Y(omega)
+    X_1 = {}
+    X_2 = {}
+    Y_1 = {}
+    Y_2 = {}
+
+    for axis in range(0, 3):
+        string = "MU_" + resp.cart[axis]
+
+        A = resp.pertbar[string]
+
+        X_2[string] = resp.solve_right(A, omega1)
+        Y_2[string] = resp.solve_left(A, omega1)
+        X_1[string] = resp.solve_right(A, -omega1)
+        Y_1[string] = resp.solve_left(A, -omega1)
+
+    # Grabbing X, Y and declaring the matrix space for LR
+    polar_AB = np.zeros((3,3))
+
+    for a in range(0, 3):
+        string_a = "MU_" + resp.cart[a]
+        for b in range(0,3):
+            string_b = "MU_" + resp.cart[b]
+            Y1_B, Y2_B, _ = Y_2[string_b]
+            X1_B, X2_B, _ = X_2[string_b]
+            polar_AB[a,b] = resp.linresp_asym(string_a, X1_B, X2_B, Y1_B, Y2_B)
+
+    print("Dynamic Polarizability Tensor @ w=0.0656 a.u.:")
+    print(polar_AB)
+    print("Average Dynamic Polarizability:")
+    polar_AB_avg = np.average([polar_AB[0,0], polar_AB[1,1], polar_AB[2,2]])
+    print(polar_AB_avg)
+
+    #validating from psi4
+    polar_XX = 9.92992070420665
+    polar_YY = 13.443740151331559
+    polar_ZZ = 11.342765745046526
+    polar_avg = 11.572142200333
+
+    assert(abs(polar_AB[0,0] - polar_XX) < 1e-8)
+    assert(abs(polar_AB[1,1] - polar_YY) < 1e-8)
+    assert(abs(polar_AB[2,2] - polar_ZZ) < 1e-8)
+    assert(abs(polar_AB_avg - polar_avg) < 1e-8)