diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 59c89a1..9ef3ec3 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -21,8 +21,8 @@ jobs: strategy: fail-fast: true matrix: - os: [macOS-latest, ubuntu-latest] - python-version: ['3.11'] + os: [macOS-latest] + python-version: ['3.8','3.11'] steps: - uses: actions/checkout@v2 @@ -49,6 +49,7 @@ jobs: hash -r env conda install psi4 python=${{ matrix.python-version }} -c conda-forge + conda install pytorch ls -l $CONDA - name: Test Psi4 Python Loading @@ -67,7 +68,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest + coverage run -m pytest -s pycc/tests/test_038_sim_quadresp.py::test_PNOpp_ccsd_OR coverage xml ls -l diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index 092977c..62ff5f0 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -17,10 +17,18 @@ class ccresponse(object): ------- linresp(): Compute a CC linear response function. + quadresp(): + Compute a CC quadratic response function. + hyperpolar(): + Compute a first electric dipole hyperpolarizability average. solve_right(): Solve the right-hand perturbed wave function equations. + solve_left(): + Solve the left-hand perturbed wave function equations. pertcheck(): Check first-order perturbed wave functions for all available perturbation operators. + pert_quadresp(): + Obtain the solutions of the right- and left-hand perturbed wave function equations for the CC quadritc response function. """ def __init__(self, ccdensity, omega1 = 0, omega2 = 0): @@ -93,6 +101,14 @@ def __init__(self, ccdensity, omega1 = 0, omega2 = 0): self.Dia = eps_occ.reshape(-1,1) - eps_vir self.Dijab = eps_occ.reshape(-1,1,1,1) + eps_occ.reshape(-1,1,1) - eps_vir.reshape(-1,1) - eps_vir + #HBAR-based denominators for simulation code + if self.ccwfn.local is not None: + self.eps_occ = np.diag(self.hbar.Hoo) + self.eps_vir = [] + for ij in range(self.ccwfn.no*self.ccwfn.no): + tmp = self.ccwfn.Local.Q[ij].T @ self.hbar.Hvv @ self.ccwfn.Local.Q[ij] + self.eps_vir.append(np.diag(self.ccwfn.Local.L[ij].T @ tmp @ self.ccwfn.Local.L[ij])) + def pertcheck(self, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis=8, start_diis=1): """ Build first-order perturbed wave functions for all available perturbations and return a dict of their converged pseudoresponse values. Primarily for testing purposes. @@ -200,6 +216,7 @@ def pertcheck(self, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis=8, print("Solving right-hand perturbed wave function for %s:" % (X_key)) X1[X_key], X2[X_key], polar = self.solve_right(self.pertbar[pertkey], -omega, e_conv, r_conv, maxiter, max_diis, start_diis) check[X_key] = polar + return check @@ -291,7 +308,1483 @@ 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 solve_right(self, pertbar, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis=8, start_diis=1): + def pert_quadresp(self, omega1, omega2, e_conv=1e-12, r_conv=1e-12, maxiter=200, max_diis=7, start_diis=1): + """ + Build first-order perturbed wave functions (left- and right-hand) for the electric dipole operator (Mu) + + Parameters + ---------- + omega1: float + First external field frequency. + omega2: float + Second external field frequency. + e_conv : float + convergence condition for the pseudoresponse value (default if 1e-13) + r_conv : float + convergence condition for perturbed wave function rmsd (default if 1e-13) + maxiter : int + maximum allowed number of iterations of the wave function equations (default is 100) + max_diis : int + maximum number of error vectors in the DIIS extrapolation (default is 8; set to 0 to deactivate) + start_diis : int + earliest iteration to start DIIS extrapolations (default is 1) + + To Do + ----- + Organize to only compute the neccesary perturbed wave functions. + """ + + #dictionaries for perturbed waves functions + self.ccpert_om1_X = {} + self.ccpert_om2_X = {} + self.ccpert_om_sum_X = {} + + self.ccpert_om1_2nd_X = {} + self.ccpert_om2_2nd_X = {} + self.ccpert_om_sum_2nd_X = {} + + self.ccpert_om1_Y = {} + self.ccpert_om2_Y = {} + self.ccpert_om_sum_Y = {} + + self.ccpert_om1_2nd_Y = {} + self.ccpert_om2_2nd_Y = {} + self.ccpert_om_sum_2nd_Y = {} + + omega_sum = -(omega1 + omega2) + + for axis in range(0, 3): + + pertkey = "MU_" + self.cart[axis] + + print("Solving right-hand perturbed wave function for omega1 %s:" % (pertkey)) + self.ccpert_om1_X[pertkey] = self.solve_right(self.pertbar[pertkey], omega1, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving left-hand perturbed wave function for %s:" % (pertkey)) + self.ccpert_om1_Y[pertkey] = self.solve_left(self.pertbar[pertkey], omega1, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving right-hand perturbed wave function for omega2 %s:" % (pertkey)) + self.ccpert_om2_X[pertkey] = self.solve_right(self.pertbar[pertkey], omega2, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving left-hand perturbed wave function for %s:" % (pertkey)) + self.ccpert_om2_Y[pertkey] = self.solve_left(self.pertbar[pertkey], omega2, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving right-hand perturbed wave function for omega_sum %s:" % (pertkey)) + self.ccpert_om_sum_X[pertkey] = self.solve_right(self.pertbar[pertkey], omega_sum, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving left-hand perturbed wave function for %s:" % (pertkey)) + self.ccpert_om_sum_Y[pertkey] = self.solve_left(self.pertbar[pertkey], omega_sum, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving right-hand perturbed wave function for -omega1 %s:" % (pertkey)) + self.ccpert_om1_2nd_X[pertkey] = self.solve_right(self.pertbar[pertkey], -omega1, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving left-hand perturbed wave function for %s:" % (pertkey)) + self.ccpert_om1_2nd_Y[pertkey] = self.solve_left(self.pertbar[pertkey], -omega1, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving right-hand perturbed wave function for -omega2 %s:" % (pertkey)) + self.ccpert_om2_2nd_X[pertkey] = self.solve_right(self.pertbar[pertkey], -omega2, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving left-hand perturbed wave function for %s:" % (pertkey)) + self.ccpert_om2_2nd_Y[pertkey] = self.solve_left(self.pertbar[pertkey], -omega2, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving right-hand perturbed wave function for -omega_sum %s:" % (pertkey)) + self.ccpert_om_sum_2nd_X[pertkey] = self.solve_right(self.pertbar[pertkey], -omega_sum, e_conv, r_conv, maxiter, max_diis, start_diis) + + print("Solving left-hand perturbed wave function for %s:" % (pertkey)) + self.ccpert_om_sum_2nd_Y[pertkey] = self.solve_left(self.pertbar[pertkey], -omega_sum, e_conv, r_conv, maxiter, max_diis, start_diis) + + def quadraticresp(self, pertkey_a, pertkey_b, pertkey_c, ccpert_X_A, ccpert_X_B, ccpert_X_C, ccpert_Y_A, ccpert_Y_B, ccpert_Y_C): + """ + Calculate the CC quadratic-response function for one-electron perturbations A,B and C at field-frequency omega1(w1) and omega2(w2). + + The quadratic response function, <>w1, generally requires the following perturbed wave functions and frequencies: + A(-w1-w2), A*(w1+w2), B(w1), B*(-w1), C(w2), C*(w2) + + Parameters + ---------- + pertkey_a: string + String identifying the one-electron perturbation, A, along a cartesian axis + pertkey_b: string + String identifying the one-electron perturbation, B, along a cartesian axis + pertkey_c: string + String identifying the one-electron perturbation, C, along a cartesian axis + ccpert_X_A: + Perturbed right-hand wave functions for A along a cartesian axis + ccpert_X_B: + + Return + ------ + hyper: float + A value of the chosen quadratic response function corresponding to a specified cartesian direction. For example, Beta_xyz. + + Notes + ----- + Only the electric dipole is used for computing second harmonic generation (SHG) where w1 and w2 are identical and optical refractivity (OR) + where w1 = -w2 + + To Do + ----- + - Expand to include all avaiable one-electron perturbations + """ + contract = self.contract + + o = self.ccwfn.o + v = self.ccwfn.v + t1 = self.ccwfn.t1 + t2 = self.ccwfn.t2 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + # Grab X and Y amplitudes corresponding to perturbation A, omega_sum + X1_A = ccpert_X_A[0] + X2_A = ccpert_X_A[1] + Y1_A = ccpert_Y_A[0] + Y2_A = ccpert_Y_A[1] + # Grab X and Y amplitudes corresponding to perturbation B, omega1 + X1_B = ccpert_X_B[0] + X2_B = ccpert_X_B[1] + Y1_B = ccpert_Y_B[0] + Y2_B = ccpert_Y_B[1] + # Grab X and Y amplitudes corresponding to perturbation C, omega2 + X1_C = ccpert_X_C[0] + X2_C = ccpert_X_C[1] + Y1_C = ccpert_Y_C[0] + Y2_C = ccpert_Y_C[1] + # Grab pert integrals + pertbar_A = self.pertbar[pertkey_a] + pertbar_B = self.pertbar[pertkey_b] + pertbar_C = self.pertbar[pertkey_c] + + #Grab H_bar, L and ERI + hbar = self.hbar + L = self.H.L + ERI = self.H. ERI + + self.hyper = 0.0 + self.LAX = 0.0 + self.LAX2 = 0.0 + self.LAX3 = 0.0 + self.LAX4 = 0.0 + self.LAX5 = 0.0 + self.LAX6 = 0.0 + self.LHX1Y1 = 0.0 + self.LHX1Y2 = 0.0 + self.LHX1X2 = 0.0 + self.LHX2Y2 = 0.0 + + # <0|L1(B)[A_bar, X1(C)]|0> + tmp = contract('ia,ic->ac', Y1_B, X1_C) + self.LAX += contract('ac,ac->',tmp, pertbar_A.Avv) + tmp = contract('ia,ka->ik', Y1_B, X1_C) + self.LAX -= contract('ik,ki->', tmp, pertbar_A.Aoo) + + # <0|L1(B)[A_bar, X2(C)]|0> + tmp = contract('ia,jb->ijab', Y1_B, pertbar_A.Aov) + + #swapaxes + self.LAX += 2.0 * contract('ijab,ijab->', tmp, X2_C) + self.LAX -= contract('ijab,ijba->',tmp, X2_C) + + # <0|L2(B)[A_bar, X1(C)]|0> + tmp = contract('ijbc,bcaj->ia', Y2_B, pertbar_A.Avvvo) + self.LAX += contract('ia,ia->', tmp, X1_C) + tmp = contract('ijab,kbij->ak', Y2_B, pertbar_A.Aovoo) + self.LAX -= contract('ak,ka->', tmp, X1_C) + # <0|L2(B)[A_bar, X2(C)]|0> + tmp = contract('ijab,kjab->ik', Y2_B, X2_C) + self.LAX -= contract('ik,ki->', tmp, pertbar_A.Aoo) + tmp = contract('ijab,ijac->bc', Y2_B, X2_C) + self.LAX += contract('bc,bc->', tmp, pertbar_A.Avv) + + self.hyper += self.LAX + + # <0|L1(C)[A_bar, X1(B)]|0> + tmp = contract('ia,ic->ac', Y1_C, X1_B) + self.LAX2 += contract('ac,ac->', tmp, pertbar_A.Avv) + tmp = contract('ia,ka->ik', Y1_C, X1_B) + self.LAX2 -= contract('ik,ki->',tmp, pertbar_A.Aoo) + + # <0|L1(C)[A_bar, X2(B)]|0> + tmp = contract('ia,jb->ijab', Y1_C, pertbar_A.Aov) + + #swapaxes + self.LAX2 += 2.0 * contract('ijab,ijab->', tmp, X2_B) + self.LAX2 -= contract('ijab,ijba->', tmp, X2_B) + + # <0|L2(C)[A_bar, X1(B)]|0> + tmp = contract('ijbc,bcaj->ia', Y2_C, pertbar_A.Avvvo) + self.LAX2 += contract('ia,ia->', tmp, X1_B) + tmp = contract('ijab,kbij->ak', Y2_C, pertbar_A.Aovoo) + self.LAX2 -= contract('ak,ka->', tmp, X1_B) + + # <0|L2(C)[A_bar, X2(B)]|0> + tmp = contract('ijab,kjab->ik', Y2_C, X2_B) + self.LAX2 -= contract('ik,ki->', tmp, pertbar_A.Aoo) + tmp = contract('ijab,ijac->bc', Y2_C, X2_B) + self.LAX2 += contract('bc,bc->', tmp, pertbar_A.Avv) + + self.hyper += self.LAX2 + + # <0|L1(A)[B_bar,X1(C)]|0> + tmp = contract('ia,ic->ac', Y1_A, X1_C) + self.LAX3 += contract('ac,ac->',tmp, pertbar_B.Avv) + tmp = contract('ia,ka->ik', Y1_A, X1_C) + self.LAX3 -= contract('ik,ki->', tmp, pertbar_B.Aoo) + # <0|L1(A)[B_bar, X2(C)]|0> + tmp = contract('ia,jb->ijab', Y1_A, pertbar_B.Aov) + + #swapaxes + self.LAX3 += 2.0 * contract('ijab,ijab->', tmp, X2_C) + self.LAX3 -= contract('ijab,ijba->', tmp, X2_C) + # <0|L2(A)[B_bar, X1(C)]|0> + tmp = contract('ijbc,bcaj->ia', Y2_A, pertbar_B.Avvvo) + self.LAX3 += contract('ia,ia->',tmp, X1_C) + tmp = contract('ijab,kbij->ak', Y2_A, pertbar_B.Aovoo) + self.LAX3 -= contract('ak,ka->', tmp, X1_C) + # <0|L2(A)[B_bar, X2(C)]|0> + tmp = contract('ijab,kjab->ik', Y2_A, X2_C) + self.LAX3 -= contract('ik,ki->', tmp, pertbar_B.Aoo) + tmp = contract('ijab,ijac->bc', Y2_A, X2_C) + self.LAX3 += contract('bc,bc->', tmp, pertbar_B.Avv) + + self.hyper += self.LAX3 + + # <0|L1(C)|[B_bar,X1(A)]|0> + tmp = contract('ia,ic->ac', Y1_C, X1_A) + self.LAX4 += contract('ac,ac->', tmp, pertbar_B.Avv) + tmp = contract('ia,ka->ik', Y1_C, X1_A) + self.LAX4 -= contract('ik,ki->', tmp, pertbar_B.Aoo) + # <0|L1(C)[B_bar, X2(A)]|0> + tmp = contract('ia,jb->ijab',Y1_C, pertbar_B.Aov) + #swapaxes + self.LAX4 += 2.0 * contract('ijab,ijab->',tmp, X2_A) + self.LAX4 -= contract('ijab,ijba->', tmp, X2_A) + # <0|L2(C)[B_bar, X1(A)]|0> + tmp = contract('ijbc,bcaj->ia', Y2_C, pertbar_B.Avvvo) + self.LAX4 += contract('ia,ia->', tmp, X1_A) + tmp = contract('ijab,kbij->ak', Y2_C, pertbar_B.Aovoo) + self.LAX4 -= contract('ak,ka->', tmp, X1_A) + # <0|L2(C)[B_bar, X2(A)]|0> + tmp = contract('ijab,kjab->ik', Y2_C, X2_A) + self.LAX4 -= contract('ik,ki->', tmp, pertbar_B.Aoo) + tmp = contract('ijab,kiba->jk', Y2_C, X2_A) + tmp = contract('ijab,ijac->bc', Y2_C, X2_A) + self.LAX4 += contract('bc,bc->', tmp, pertbar_B.Avv) + + self.hyper += self.LAX4 + + # <0|L1(A)[C_bar,X1(B)]|0> + tmp = contract('ia,ic->ac', Y1_A, X1_B) + self.LAX5 += contract('ac,ac->', tmp, pertbar_C.Avv) + tmp = contract('ia,ka->ik', Y1_A, X1_B) + self.LAX5 -= contract('ik,ki->', tmp, pertbar_C.Aoo) + # <0|L1(A)[C_bar, X2(B)]|0> + tmp = contract('ia,jb->ijab', Y1_A, pertbar_C.Aov) + + #swapaxes + self.LAX5 += 2.0 * contract('ijab,ijab->', tmp, X2_B) + self.LAX5 -= contract('ijab,ijba->', tmp, X2_B) + # <0|L2(A)[C_bar, X1(B)]|0> + tmp = contract('ijbc,bcaj->ia', Y2_A, pertbar_C.Avvvo) + self.LAX5 += contract('ia,ia->', tmp, X1_B) + tmp = contract('ijab,kbij->ak', Y2_A, pertbar_C.Aovoo) + self.LAX5 -= contract('ak,ka->', tmp, X1_B) + # <0|L2(A)[C_bar, X2(B)]|0> + tmp = contract('ijab,kjab->ik', Y2_A, X2_B) + self.LAX5 -= contract('ik,ki->', tmp, pertbar_C.Aoo) + tmp = contract('ijab,ijac->bc', Y2_A, X2_B) + self.LAX5 += contract('bc,bc->', tmp, pertbar_C.Avv) + + self.hyper += self.LAX5 + + # <0|L1(B)|[C_bar,X1(A)]|0> + tmp = contract('ia,ic->ac', Y1_B, X1_A) + self.LAX6 += contract('ac,ac->', tmp, pertbar_C.Avv) + tmp = contract('ia,ka->ik', Y1_B, X1_A) + self.LAX6 -= contract('ik,ki->', tmp, pertbar_C.Aoo) + # <0|L1(B)[C_bar, X2(A)]|0> + tmp = contract('ia,jb->ijab', Y1_B, pertbar_C.Aov) + + #swapaxes + self.LAX6 += 2.0 * contract('ijab,ijab->', tmp, X2_A) + self.LAX6 -= contract('ijab,ijba->', tmp, X2_A) + # <0|L2(B)[C_bar, X1(A)]|0> + tmp = contract('ijbc,bcaj->ia', Y2_B, pertbar_C.Avvvo) + self.LAX6 += contract('ia,ia->', tmp, X1_A) + tmp = contract('ijab,kbij->ak', Y2_B, pertbar_C.Aovoo) + self.LAX6 -= contract('ak,ka->', tmp, X1_A) + # <0|L2(B)[C_bar, X2(A)]|0> + tmp = contract('ijab,kjab->ik', Y2_B, X2_A) + self.LAX6 -= np.einsum('ik,ki->', tmp, pertbar_C.Aoo) + tmp = contract('ijab,ijac->bc', Y2_B, X2_A) + self.LAX6 += contract('bc,bc->', tmp, pertbar_C.Avv) + + self.hyper += self.LAX6 + + self.Fz1 = 0 + self.Fz2 = 0 + self.Fz3 = 0 + + # <0|L1(0)[[A_bar,X1(B)],X1(C)]|0> + tmp = contract('ia,ja->ij', X1_B, pertbar_A.Aov) + tmp2 = contract('ib,jb->ij', l1, X1_C) + self.Fz1 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('jb,ib->ij', X1_C, pertbar_A.Aov) + tmp2 = contract('ia,ja->ij', X1_B, l1) + self.Fz1 -= contract('ij,ij->', tmp2, tmp) + + # <0|L2(0)[[A_bar,X1(B)],X2(C)]|0> + tmp = contract('ia,ja->ij', X1_B, pertbar_A.Aov) + tmp2 = contract('jkbc,ikbc->ij', X2_C, l2) + self.Fz1 -= contract('ij,ij->',tmp2,tmp) + + tmp = contract('ia,jkac->jkic', X1_B, l2) + tmp = contract('jkbc,jkic->ib', X2_C, tmp) + self.Fz1 -= contract('ib,ib->', tmp, pertbar_A.Aov) + + # <0|L2(0)[[A_bar,X2(B)],X1(C)]|0> + tmp = contract('ia,ja->ij', X1_C, pertbar_A.Aov) + tmp2 = contract('jkbc,ikbc->ij', X2_B, l2) + self.Fz1 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('ia,jkac->jkic', X1_C, l2) + tmp = contract('jkbc,jkic->ib', X2_B, tmp) + self.Fz1 -= contract('ib,ib->', tmp, pertbar_A.Aov) + + # <0|L1(0)[B_bar,X1(A)],X1(C)]|0> + tmp = contract('ia,ja->ij', X1_A, pertbar_B.Aov) + tmp2 = contract('ib,jb->ij', l1, X1_C) + self.Fz2 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('jb,ib->ij', X1_C, pertbar_B.Aov) + tmp2 = contract('ia,ja->ij', X1_A, l1) + self.Fz2 -= contract('ij,ij->', tmp2, tmp) + + # <0|L2(0)[[B_bar,X1(A)],X2(C)]|0> + tmp = contract('ia,ja->ij', X1_A, pertbar_B.Aov) + tmp2 = contract('jkbc,ikbc->ij', X2_C, l2) + self.Fz2 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('ia,jkac->jkic', X1_A, l2) + tmp = contract('jkbc,jkic->ib', X2_C, tmp) + self.Fz2 -= contract('ib,ib->', tmp, pertbar_B.Aov) + + # <0|L2(0)[[B_bar,X2(A)],X1(C)]|0> + tmp = contract('ia,ja->ij', X1_C, pertbar_B.Aov) + tmp2 = contract('jkbc,ikbc->ij', X2_A, l2) + self.Fz2 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('ia,jkac->jkic', X1_C, l2) + tmp = contract('jkbc,jkic->ib', X2_A, tmp) + self.Fz2 -= contract('ib,ib->', tmp, pertbar_B.Aov) + + # <0|L1(0)[C_bar,X1(A)],X1(B)]|0> + tmp = contract('ia,ja->ij', X1_A, pertbar_C.Aov) + tmp2 = contract('ib,jb->ij', l1, X1_B) + self.Fz3 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('jb,ib->ij', X1_B, pertbar_C.Aov) + tmp2 = contract('ia,ja->ij', X1_A, l1) + self.Fz3 -= contract('ij,ij->', tmp2, tmp) + + # <0|L2(0)[[C_bar,X1(A)],X2(B)]|0> + tmp = contract('ia,ja->ij', X1_A, pertbar_C.Aov) + tmp2 = contract('jkbc,ikbc->ij', X2_B, l2) + self.Fz3 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('ia,jkac->jkic', X1_A, l2) + tmp = contract('jkbc,jkic->ib', X2_B, tmp) + self.Fz3 -= contract('ib,ib->', tmp, pertbar_C.Aov) + + # <0|L2(0)[[C_bar,X2(A)],X1(B)]|0> + tmp = contract('ia,ja->ij', X1_B, pertbar_C.Aov) + tmp2 = contract('jkbc,ikbc->ij', X2_A, l2) + self.Fz3 -= contract('ij,ij->', tmp2, tmp) + + tmp = contract('ia,jkac->jkic', X1_B, l2) + tmp = contract('jkbc,jkic->ib', X2_A, tmp) + self.Fz3 -= contract('ib,ib->', tmp, pertbar_C.Aov) + + self.hyper += self.Fz1+self.Fz2+self.Fz3 + + self.G = 0 + # + tmp = contract('ia,ijac->jc', X1_A, L[o,o,v,v]) + tmp = contract('kc,jc->jk', X1_C, tmp) + tmp2 = contract('jb,kb->jk', X1_B, l1) + self.G -= contract('jk,jk->', tmp2, tmp) + + tmp = contract('ia,ikab->kb', X1_A, L[o,o,v,v]) + tmp = contract('jb,kb->jk', X1_B, tmp) + tmp2 = contract('jc,kc->jk', l1, X1_C) + self.G -= contract('jk,jk->', tmp2, tmp) + + tmp = contract('jb,jkba->ka', X1_B, L[o,o,v,v]) + tmp = contract('ia,ka->ki', X1_A, tmp) + tmp2 = contract('kc,ic->ki', X1_C, l1) + self.G -= contract('ki,ki->', tmp2, tmp) + + tmp = contract('jb,jibc->ic', X1_B, L[o,o,v,v]) + tmp = contract('kc,ic->ki', X1_C, tmp) + tmp2 = contract('ka,ia->ki', l1, X1_A) + self.G -= contract('ki,ki->', tmp2, tmp) + + tmp = contract('kc,kicb->ib', X1_C, L[o,o,v,v]) + tmp = contract('jb,ib->ji', X1_B, tmp) + tmp2 = contract('ja,ia->ji', l1, X1_A) + self.G -= contract('ji,ji->', tmp2, tmp) + + tmp = contract('kc,kjca->ja', X1_C, L[o,o,v,v]) + tmp = contract('ia,ja->ji', X1_A, tmp) + tmp2 = contract('jb,ib->ji', X1_B, l1) + self.G -= contract('ji,ji->', tmp2, tmp) + + # + tmp = contract('jb,klib->klij', X1_A, hbar.Hooov) + tmp2 = contract('ld,ijcd->ijcl', X1_C, l2) + tmp2 = contract('kc,ijcl->ijkl', X1_B, tmp2) + self.G += contract('ijkl,klij->', tmp2, tmp) + + tmp = contract('jb,lkib->lkij', X1_A, hbar.Hooov) + tmp2 = contract('ld,ijdc->ijlc', X1_C, l2) + tmp2 = contract('kc,ijlc->ijlk', X1_B, tmp2) + self.G += contract('ijlk,lkij->', tmp2, tmp) + + tmp = contract('kc,jlic->jlik', X1_B, hbar.Hooov) + tmp2 = contract('jb,ikbd->ikjd', X1_A, l2) + tmp2 = contract('ld,ikjd->ikjl', X1_C, tmp2) + self.G += contract('ikjl,jlik->', tmp2, tmp) + + tmp = contract('kc,ljic->ljik', X1_B, hbar.Hooov) + tmp2 = contract('jb,ikdb->ikdj', X1_A, l2) + tmp2 = contract('ld,ikdj->iklj', X1_C, tmp2) + self.G += contract('iklj,ljik->', tmp2, tmp) + + tmp = contract('ld,jkid->jkil', X1_C, hbar.Hooov) + tmp2 = contract('jb,ilbc->iljc', X1_A, l2) + tmp2 = contract('kc,iljc->iljk', X1_B, tmp2) + self.G += contract('iljk,jkil->', tmp2, tmp) + + tmp = contract('ld,kjid->kjil', X1_C, hbar.Hooov) + tmp2 = contract('jb,ilcb->ilcj', X1_A, l2) + tmp2 = contract('kc,ilcj->ilkj', X1_B, tmp2) + self.G += contract('ilkj,kjil->', tmp2, tmp) + + tmp = contract('jb,albc->aljc', X1_A, hbar.Hvovv) + tmp = contract('kc,aljc->aljk', X1_B, tmp) + tmp2 = contract('ld,jkad->jkal', X1_C, l2) + self.G -= contract('jkal,aljk->', tmp2, tmp) + + tmp = contract('jb,alcb->alcj', X1_A, hbar.Hvovv) + tmp = contract('kc,alcj->alkj', X1_B, tmp) + tmp2 = contract('ld,jkda->jkla', X1_C, l2) + self.G -= contract('jkla,alkj->', tmp2, tmp) + + tmp = contract('jb,akbd->akjd', X1_A, hbar.Hvovv) + tmp = contract('ld,akjd->akjl', X1_C, tmp) + tmp2 = contract('kc,jlac->jlak', X1_B, l2) + self.G -= contract('jlak,akjl->', tmp2, tmp) + + tmp = contract('jb,akdb->akdj', X1_A, hbar.Hvovv) + tmp = contract('ld,akdj->aklj', X1_C, tmp) + tmp2 = contract('kc,jlca->jlka', X1_B, l2) + self.G -= contract('jlka,aklj->', tmp2, tmp) + + tmp = contract('kc,ajcd->ajkd', X1_B, hbar.Hvovv) + tmp = contract('ld,ajkd->ajkl', X1_C, tmp) + tmp2 = contract('jb,klab->klaj', X1_A, l2) + self.G -= contract('klaj,ajkl->', tmp2, tmp) + + tmp = contract('kc,ajdc->ajdk', X1_B, hbar.Hvovv) + tmp = contract('ld,ajdk->ajlk', X1_C, tmp) + tmp2 = contract('jb,klba->klja', X1_A, l2) + self.G -= contract('klja,ajlk->', tmp2, tmp) + + # + tmp = contract('kc,jlbc->jlbk', X1_B, l2) + tmp2 = contract('ld,ikad->ikal', X1_C, L[o,o,v,v]) + tmp2 = contract('ijab,ikal->jlbk', X2_A, tmp2) + self.G -= contract('jlbk,jlbk->', tmp, tmp2) + + tmp = contract('ld,jkbd->jkbl', X1_C, l2) + tmp2 = contract('kc,ilac->ilak', X1_B, L[o,o,v,v]) + tmp2 = contract('ijab,ilak->jkbl', X2_A, tmp2) + self.G -= contract('jkbl,jkbl->',tmp,tmp2) + + tmp = contract('ijab,jibd->ad', X2_A, l2) + tmp = contract('ld,ad->la', X1_C, tmp) + tmp2 = contract('klca,kc->la', L[o,o,v,v], X1_B) + self.G -= contract('la,la->', tmp, tmp2) + + tmp = contract('ijab,jlba->il', X2_A, l2) + tmp2 = contract('kc,kicd->id', X1_B, L[o,o,v,v]) + tmp2 = contract('ld,id->il', X1_C, tmp2) + self.G -= contract('il,il->', tmp, tmp2) + + tmp = contract('ijab,jkba->ik', X2_A, l2) + tmp2 = contract('ld,lidc->ic', X1_C, L[o,o,v,v]) + tmp2 = contract('kc,ic->ik', X1_B, tmp2) + self.G -= contract('ik,ik->', tmp, tmp2) + + tmp = contract('ijab,jibc->ac', X2_A, l2) + tmp = contract('ac,kc->ka', tmp, X1_B) + tmp2 = contract('ld,lkda->ka', X1_C, L[o,o,v,v]) + self.G -= contract('ka,ka->', tmp, tmp2) + + tmp = contract('ijab,klab->ijkl',X2_A, ERI[o,o,v,v]) + tmp2 = contract('kc,ijcd->ijkd', X1_B, l2) + tmp2 = contract('ld,ijkd->ijkl', X1_C, tmp2) + self.G += contract('ijkl,ijkl->',tmp,tmp2) + + tmp = contract('kc,jlac->jlak', X1_B, ERI[o,o,v,v]) + tmp = contract('ijab,jlak->ilbk', X2_A, tmp) + tmp2 = contract('ikbd,ld->ilbk', l2, X1_C) + self.G += contract('ilbk,ilbk->', tmp, tmp2) + + tmp = contract('kc,ljac->ljak', X1_B, ERI[o,o,v,v]) + tmp = contract('ijab,ljak->ilbk', X2_A, tmp) + tmp2 = contract('ikdb,ld->ilbk', l2, X1_C) + self.G += contract('ilbk,ilbk->', tmp, tmp2) + + tmp = contract('ld,jkad->jkal', X1_C, ERI[o,o,v,v]) + tmp = contract('ijab,jkal->ikbl', X2_A, tmp) + tmp2 = contract('kc,ilbc->ilbk', X1_B, l2) + self.G += contract('ikbl,ilbk->', tmp, tmp2) + + tmp = contract('ld,kjad->kjal', X1_C, ERI[o,o,v,v]) + tmp = contract('ijab,kjal->iklb', X2_A, tmp) + tmp2 = contract('kc,ilcb->ilkb', X1_B, l2) + self.G += contract('iklb,ilkb->', tmp, tmp2) + + tmp = contract('kc,ijcd->ijkd', X1_B, ERI[o,o,v,v]) + tmp = contract('ld,ijkd->ijkl', X1_C, tmp) + tmp2 = contract('ijab,klab->ijkl', X2_A, l2) + self.G += contract('ijkl,ijkl->', tmp, tmp2) + + # + tmp = contract('kc,jlbc->jlbk', X1_A, l2) + tmp2 = contract('ld,ikad->ikal', X1_C, L[o,o,v,v]) + tmp2 = contract('ijab,ikal->jlbk', X2_B, tmp2) + self.G -= contract('jlbk,jlbk->', tmp, tmp2) + + tmp = contract('ld,jkbd->jkbl', X1_C, l2) + tmp2 = contract('kc,ilac->ilak', X1_A, L[o,o,v,v]) + tmp2 = contract('ijab,ilak->jkbl',X2_B, tmp2) + self.G -= contract('jkbl,jkbl->', tmp, tmp2) + + tmp = contract('ijab,jibd->ad', X2_B, l2) + tmp = contract('ld,ad->la', X1_C, tmp) + tmp2 = contract('klca,kc->la', L[o,o,v,v], X1_A) + self.G -= contract('la,la->', tmp, tmp2) + + tmp = contract('ijab,jlba->il', X2_B, l2) + tmp2 = contract('kc,kicd->id', X1_A, L[o,o,v,v]) + tmp2 = contract('ld,id->il', X1_C, tmp2) + self.G -= contract('il,il->', tmp, tmp2) + + tmp = contract('ijab,jkba->ik', X2_B, l2) + tmp2 = contract('ld,lidc->ic', X1_C, L[o,o,v,v]) + tmp2 = contract('kc,ic->ik', X1_A, tmp2) + self.G -= contract('ik,ik->', tmp, tmp2) + + tmp = contract('ijab,jibc->ac', X2_B, l2) + tmp = contract('ac,kc->ka', tmp, X1_A) + tmp2 = contract('ld,lkda->ka', X1_C, L[o,o,v,v]) + self.G -= contract('ka,ka->', tmp, tmp2) + + tmp = contract('ijab,klab->ijkl', X2_B, ERI[o,o,v,v]) + tmp2 = contract('kc,ijcd->ijkd', X1_A, l2) + tmp2 = contract('ld,ijkd->ijkl', X1_C, tmp2) + self.G += contract('ijkl,ijkl->', tmp, tmp2) + + tmp = contract('kc,jlac->jlak', X1_A, ERI[o,o,v,v]) + tmp = contract('ijab,jlak->ilbk', X2_B, tmp) + tmp2 = contract('ikbd,ld->ilbk', l2, X1_C) + self.G += contract('ilbk,ilbk->', tmp, tmp2) + + tmp = contract('kc,ljac->ljak', X1_A, ERI[o,o,v,v]) + tmp = contract('ijab,ljak->ilbk', X2_B, tmp) + tmp2 = contract('ikdb,ld->ilbk', l2, X1_C) + self.G += contract('ilbk,ilbk->', tmp, tmp2) + + tmp = contract('ld,jkad->jkal', X1_C, ERI[o,o,v,v]) + tmp = contract('ijab,jkal->ikbl', X2_B, tmp) + tmp2 = contract('kc,ilbc->ilbk', X1_A, l2) + self.G += contract('ikbl,ilbk->', tmp, tmp2) + + tmp = contract('ld,kjad->kjal', X1_C, ERI[o,o,v,v]) + tmp = contract('ijab,kjal->iklb', X2_B, tmp) + tmp2 = contract('kc,ilcb->ilkb', X1_A, l2) + self.G += contract('iklb,ilkb->', tmp, tmp2) + + tmp = contract('kc,ijcd->ijkd', X1_A, ERI[o,o,v,v]) + tmp = contract('ld,ijkd->ijkl', X1_C, tmp) + tmp2 = contract('ijab,klab->ijkl', X2_B, l2) + self.G += contract('ijkl,ijkl->', tmp, tmp2) + + # + tmp = contract('kc,jlbc->jlbk', X1_A, l2) + tmp2 = contract('ld,ikad->ikal', X1_B, L[o,o,v,v]) + tmp2 = contract('ijab,ikal->jlbk', X2_C, tmp2) + self.G -= contract('jlbk,jlbk->', tmp, tmp2) + + tmp = contract('ld,jkbd->jkbl', X1_B, l2) + tmp2 = contract('kc,ilac->ilak', X1_A, L[o,o,v,v]) + tmp2 = contract('ijab,ilak->jkbl', X2_C, tmp2) + self.G -= contract('jkbl,jkbl->', tmp, tmp2) + + tmp = contract('ijab,jibd->ad', X2_C, l2) + tmp = contract('ld,ad->la', X1_B, tmp) + tmp2 = contract('klca,kc->la', L[o,o,v,v], X1_A) + self.G -= contract('la,la->', tmp, tmp2) + + tmp = contract('ijab,jlba->il', X2_C, l2) + tmp2 = contract('kc,kicd->id', X1_A, L[o,o,v,v]) + tmp2 = contract('ld,id->il', X1_B, tmp2) + self.G -= contract('il,il->', tmp, tmp2) + + tmp = contract('ijab,jkba->ik', X2_C, l2) + tmp2 = contract('ld,lidc->ic', X1_B, L[o,o,v,v]) + tmp2 = contract('kc,ic->ik', X1_A, tmp2) + self.G -= contract('ik,ik->', tmp, tmp2) + + tmp = contract('ijab,jibc->ac', X2_C, l2) + tmp = contract('ac,kc->ka', tmp, X1_A) + tmp2 = contract('ld,lkda->ka', X1_B, L[o,o,v,v]) + self.G -= contract('ka,ka->', tmp, tmp2) + + tmp = contract('ijab,klab->ijkl', X2_C, ERI[o,o,v,v]) + tmp2 = contract('kc,ijcd->ijkd', X1_A, l2) + tmp2 = contract('ld,ijkd->ijkl', X1_B, tmp2) + self.G += contract('ijkl,ijkl->', tmp, tmp2) + + tmp = contract('kc,jlac->jlak', X1_A, ERI[o,o,v,v]) + tmp = contract('ijab,jlak->ilbk', X2_C, tmp) + tmp2 = contract('ikbd,ld->ilbk', l2, X1_B) + self.G += contract('ilbk,ilbk->', tmp, tmp2) + + tmp = contract('kc,ljac->ljak', X1_A, ERI[o,o,v,v]) + tmp = contract('ijab,ljak->ilbk', X2_C, tmp) + tmp2 = contract('ikdb,ld->ilbk', l2, X1_B) + self.G += contract('ilbk,ilbk->', tmp, tmp2) + + tmp = contract('ld,jkad->jkal', X1_B, ERI[o,o,v,v]) + tmp = contract('ijab,jkal->ikbl', X2_C, tmp) + tmp2 = contract('kc,ilbc->ilbk', X1_A, l2) + self.G += contract('ikbl,ilbk->', tmp, tmp2) + + tmp = contract('ld,kjad->kjal', X1_B, ERI[o,o,v,v]) + tmp = contract('ijab,kjal->iklb', X2_C, tmp) + tmp2 = contract('kc,ilcb->ilkb', X1_A, l2) + self.G += contract('iklb,ilkb->', tmp, tmp2) + + tmp = contract('kc,ijcd->ijkd', X1_A, ERI[o,o,v,v]) + tmp = contract('ld,ijkd->ijkl', X1_B, tmp) + tmp2 = contract('ijab,klab->ijkl', X2_C, l2) + self.G += contract('ijkl,ijkl->', tmp, tmp2) + + self.hyper += self.G + + self.Bcon1 = 0 + # + tmp = -1.0* contract('jc,kb->jkcb', hbar.Hov, Y1_A) + tmp -= contract('jc,kb->jkcb', Y1_A, hbar.Hov) + + #swapaxes + tmp -= 2.0* contract('kjib,ic->jkcb', hbar.Hooov, Y1_A) + tmp += contract('jkib,ic->jkcb', hbar.Hooov, Y1_A) + + #swapaxes + tmp -= 2.0* contract('jkic,ib->jkcb', hbar.Hooov, Y1_A) + tmp += np.einsum('kjic,ib->jkcb', hbar.Hooov, Y1_A) + + # swapaxes + tmp += 2.0* contract('ajcb,ka->jkcb', hbar.Hvovv, Y1_A) + tmp -= contract('ajbc,ka->jkcb', hbar.Hvovv, Y1_A) + + # swapaxes + tmp += 2.0* contract('akbc,ja->jkcb', hbar.Hvovv, Y1_A) + tmp -= contract('akcb,ja->jkcb', hbar.Hvovv, Y1_A) + + tmp2 = contract('miae,me->ia', tmp, X1_B) + self.Bcon1 += contract('ia,ia->', tmp2, X1_C) + + # + tmp = -1.0* contract('janc,nkba->jckb', hbar.Hovov, Y2_A) + tmp -= contract('kanb,njca->jckb', hbar.Hovov, Y2_A) + tmp -= contract('jacn,nkab->jckb', hbar.Hovvo, Y2_A) + tmp -= contract('kabn,njac->jckb', hbar.Hovvo, Y2_A) + tmp += 0.5* contract('fabc,jkfa->jckb', hbar.Hvvvv, Y2_A) + tmp += 0.5* contract('facb,kjfa->jckb', hbar.Hvvvv, Y2_A) + tmp += 0.5* contract('kjin,nibc->jckb', hbar.Hoooo, Y2_A) + tmp += 0.5* contract('jkin,nicb->jckb', hbar.Hoooo, Y2_A) + tmp2 = contract('iema,me->ia', tmp, X1_B) + self.Bcon1 += contract('ia,ia->', tmp2, X1_C) + + tmp = contract('ijab,ijdb->ad', t2, Y2_A) + tmp = contract('ld,ad->la', X1_C, tmp) + tmp = contract('la,klca->kc', tmp, L[o,o,v,v]) + self.Bcon1 -= contract('kc,kc->', tmp, X1_B) + + tmp = contract('ijab,jlba->il', t2, Y2_A) + tmp2 = contract('kc,kicd->id', X1_B, L[o,o,v,v]) + tmp2 = contract('id,ld->il', tmp2, X1_C) + self.Bcon1 -= contract('il,il->', tmp2, tmp) + + tmp = contract('ijab,jkba->ik', t2, Y2_A) + tmp2 = contract('ld,lidc->ic', X1_C, L[o,o,v,v]) + tmp2 = contract('ic,kc->ik', tmp2, X1_B) + self.Bcon1 -= contract('ik,ik->', tmp2, tmp) + + tmp = contract('ijab,ijcb->ac', t2, Y2_A) + tmp = contract('kc,ac->ka', X1_B, tmp) + tmp2 = contract('ld,lkda->ka', X1_C, L[o,o,v,v]) + self.Bcon1 -= contract('ka,ka->', tmp2, tmp) + + # + tmp = contract("klcd,ijcd->ijkl", X2_C, Y2_A) + tmp = contract("ijkl,ijab->klab", tmp, X2_B) + self.Bcon1 += 0.5* contract('klab,klab->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,ikbd->jkad", X2_B, Y2_A) + tmp = contract("jkad,klcd->jlac", tmp, X2_C) + self.Bcon1 += contract('jlac,jlac->',tmp, ERI[o,o,v,v]) + + tmp = contract("klcd,ikdb->licb", X2_C, Y2_A) + tmp = contract("licb,ijab->ljca", tmp, X2_B) + self.Bcon1 += contract('ljca,ljac->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,klab->ijkl", X2_B, Y2_A) + tmp = contract("ijkl,klcd->ijcd", tmp, X2_C) + self.Bcon1 += 0.5* contract('ijcd,ijcd->',tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,ijac->bc", X2_B, L[o,o,v,v]) + tmp = contract("bc,klcd->klbd", tmp, X2_C) + self.Bcon1 -= contract("klbd,klbd->", tmp, Y2_A) + tmp = contract("ijab,ikab->jk", X2_B, L[o,o,v,v]) + tmp = contract("jk,klcd->jlcd", tmp, X2_C) + self.Bcon1 -= contract("jlcd,jlcd->", tmp, Y2_A) + tmp = contract("ikbc,klcd->ilbd", L[o,o,v,v], X2_C) + tmp = contract("ilbd,ijab->jlad", tmp, X2_B) + self.Bcon1 -= contract("jlad,jlad->", tmp, Y2_A) + tmp = contract("ijab,jlbc->ilac", X2_B, Y2_A) + tmp = contract("ilac,klcd->ikad", tmp, X2_C) + self.Bcon1 -= contract("ikad,ikad->", tmp, L[o,o,v,v]) + tmp = contract("klca,klcd->ad", L[o,o,v,v], X2_C) + tmp = contract("ad,ijdb->ijab", tmp, Y2_A) + self.Bcon1 -= contract("ijab,ijab->", tmp, X2_B) + tmp = contract("kicd,klcd->il",L[o,o,v,v], X2_C) + tmp = contract("ijab,il->ljab", X2_B, tmp) + self.Bcon1 -= contract("ljab,ljab->", tmp, Y2_A) + + tmp = contract("klcd,ikac->lida", X2_C, Y2_A) + tmp = contract("lida,jlbd->ijab", tmp, L[o,o,v,v]) + self.Bcon1 += 2.0* contract("ijab,ijab->", tmp, X2_B) + + # + tmp = 2.0* contract("jkbc,kc->jb", X2_C, Y1_A) + tmp -= contract("jkcb,kc->jb", X2_C, Y1_A) + tmp = contract('ijab,jb->ia', L[o,o,v,v], tmp) + self.Bcon1 += contract("ia,ia->", tmp, X1_B) + + tmp = contract("jkbc,jkba->ca", X2_C, L[o,o,v,v]) + tmp = contract("ia,ca->ic", X1_B, tmp) + self.Bcon1 -= contract("ic,ic->", tmp, Y1_A) + + tmp = contract("jkbc,jibc->ki", X2_C, L[o,o,v,v]) + tmp = contract("ki,ia->ka", tmp, X1_B) + self.Bcon1 -= contract("ka,ka->", tmp, Y1_A) + + # + tmp = contract("klcd,lkdb->cb", X2_C, Y2_A) + tmp = contract("jb,cb->jc", X1_B, tmp) + self.Bcon1 -= contract("jc,jc->", tmp, hbar.Hov) + + tmp = contract("klcd,ljdc->kj", X2_C, Y2_A) + tmp = contract("kj,jb->kb", tmp, X1_B) + self.Bcon1 -= contract("kb,kb->", tmp, hbar.Hov) + + tmp = contract('lkda,klcd->ac', Y2_A, X2_C) + tmp2 = contract('jb,ajcb->ac', X1_B, hbar.Hvovv) + self.Bcon1 += 2.0* contract('ac,ac->', tmp, tmp2) + + tmp = contract('lkda,klcd->ac', Y2_A, X2_C) + tmp2 = contract('jb,ajbc->ac', X1_B, hbar.Hvovv) + self.Bcon1 -= contract('ac,ac->', tmp, tmp2) + + tmp = contract('jb,ljda->lbda', X1_B, Y2_A) + tmp2 = 2.0* contract('klcd,akbc->ldab', X2_C, hbar.Hvovv) + tmp2 -= contract('klcd,akcb->ldab', X2_C, hbar.Hvovv) + self.Bcon1 += contract('lbda,ldab->', tmp, tmp2) + + tmp = contract('ia,fkba->fkbi', X1_B, hbar.Hvovv) + tmp = contract('fkbi,jifc->kjbc', tmp, Y2_A) + self.Bcon1 -= contract('jkbc,kjbc->', X2_C, tmp) + + tmp = contract('ia,fjac->fjic', X1_B, hbar.Hvovv) + tmp = contract('fjic,ikfb->jkbc', tmp, Y2_A) + self.Bcon1 -= contract('jkbc,jkbc->', X2_C, tmp) + + tmp = contract('ia,jkfa->jkfi', X1_B, Y2_A) + tmp2 = contract('jkbc,fibc->jkfi', X2_C, hbar.Hvovv) + self.Bcon1 -= contract('jkfi,jkfi->', tmp2, tmp) + + tmp = contract('jb,kjib->ki', X1_B, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_C, Y2_A) + self.Bcon1 -= 2.0*contract('ki,ki->', tmp, tmp2) + + tmp = contract('jb,jkib->ki', X1_B, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_C, Y2_A) + self.Bcon1 += contract('ki,ki->', tmp, tmp2) + + tmp = 2.0* contract('jkic,klcd->jild', hbar.Hooov, X2_C) + tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_C) + tmp = contract('jild,jb->bild', tmp, X1_B) + self.Bcon1 -= contract('bild,ilbd->', tmp, Y2_A) + + tmp = contract('ia,jkna->jkni', X1_B, hbar.Hooov) + tmp2 = contract('jkbc,nibc->jkni', X2_C, Y2_A) + self.Bcon1 += contract('jkni,jkni->', tmp2, tmp) + + tmp = contract('ia,nkab->nkib', X1_B, Y2_A) + tmp = contract('jkbc,nkib->jnic', X2_C, tmp) + self.Bcon1 += contract('jnic,ijnc->', tmp, hbar.Hooov) + + tmp = contract('ia,nkba->nkbi', X1_B, Y2_A) + tmp = contract('jkbc,nkbi->jnci', X2_C, tmp) + self.Bcon1 += contract('jnci,jinc->', tmp, hbar.Hooov) + + # + #swapaxes + tmp = 2.0* contract("jkbc,kc->jb", X2_B, Y1_A) + tmp -= contract("jkcb,kc->jb", X2_B, Y1_A) + tmp = contract('ijab,jb->ia', L[o,o,v,v], tmp) + self.Bcon1 += contract("ia,ia->", tmp, X1_C) + + tmp = contract("jkbc,jkba->ca", X2_B, L[o,o,v,v]) + tmp = contract("ia,ca->ic", X1_C, tmp) + self.Bcon1 -= contract("ic,ic->", tmp, Y1_A) + + tmp = contract("jkbc,jibc->ki", X2_B, L[o,o,v,v]) + tmp = contract("ki,ia->ka", tmp, X1_C) + self.Bcon1 -= contract("ka,ka->", tmp, Y1_A) + + # + tmp = contract("klcd,lkdb->cb", X2_B, Y2_A) + tmp = contract("jb,cb->jc", X1_C, tmp) + self.Bcon1 -= contract("jc,jc->", tmp, hbar.Hov) + + tmp = contract("klcd,ljdc->kj", X2_B, Y2_A) + tmp = contract("kj,jb->kb", tmp, X1_C) + self.Bcon1 -= contract("kb,kb->", tmp, hbar.Hov) + + tmp = contract('lkda,klcd->ac', Y2_A, X2_B) + tmp2 = contract('jb,ajcb->ac', X1_C, hbar.Hvovv) + self.Bcon1 += 2.0* contract('ac,ac->', tmp, tmp2) + + tmp = contract('lkda,klcd->ac', Y2_A, X2_B) + tmp2 = contract('jb,ajbc->ac', X1_C, hbar.Hvovv) + self.Bcon1 -= contract('ac,ac->', tmp, tmp2) + + tmp = contract('jb,ljda->lbda', X1_C, Y2_A) + + #swapaxes + tmp2 = 2.0* contract('klcd,akbc->ldab', X2_B, hbar.Hvovv) + tmp2 -= contract('klcd,akcb->ldab', X2_B, hbar.Hvovv) + self.Bcon1 += contract('lbda,ldab->', tmp, tmp2) + + tmp = contract('ia,fkba->fkbi', X1_C, hbar.Hvovv) + tmp = contract('fkbi,jifc->kjbc', tmp, Y2_A) + self.Bcon1 -= contract('jkbc,kjbc->', X2_B, tmp) + + tmp = contract('ia,fjac->fjic', X1_C, hbar.Hvovv) + tmp = contract('fjic,ikfb->jkbc', tmp, Y2_A) + self.Bcon1 -= contract('jkbc,jkbc->', X2_B, tmp) + + tmp = contract('ia,jkfa->jkfi', X1_C, Y2_A) + tmp2 = contract('jkbc,fibc->jkfi', X2_B, hbar.Hvovv) + self.Bcon1 -= contract('jkfi,jkfi->', tmp2, tmp) + + tmp = contract('jb,kjib->ki', X1_C, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_B, Y2_A) + self.Bcon1 -= 2.0* contract('ki,ki->', tmp, tmp2) + + tmp = contract('jb,jkib->ki', X1_C, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_B, Y2_A) + self.Bcon1 += 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_C) + self.Bcon1 -= contract('bild,ilbd->', tmp, Y2_A) + + tmp = contract('ia,jkna->jkni', X1_C, hbar.Hooov) + tmp2 = contract('jkbc,nibc->jkni', X2_B, Y2_A) + self.Bcon1 += contract('jkni,jkni->', tmp2, tmp) + + tmp = contract('ia,nkab->nkib', X1_C, Y2_A) + tmp = contract('jkbc,nkib->jnic', X2_B, tmp) + self.Bcon1 += contract('jnic,ijnc->', tmp, hbar.Hooov) + + tmp = contract('ia,nkba->nkbi', X1_C, Y2_A) + tmp = contract('jkbc,nkbi->jnci', X2_B, tmp) + self.Bcon1 += contract('jnci,jinc->', tmp, hbar.Hooov) + + self.Bcon2 = 0 + # + tmp = -1.0* contract('jc,kb->jkcb', hbar.Hov, Y1_B) + tmp -= contract('jc,kb->jkcb', Y1_B, hbar.Hov) + + #swapaxes + tmp -= 2.0* contract('kjib,ic->jkcb',hbar.Hooov, Y1_B) + tmp += contract('jkib,ic->jkcb', hbar.Hooov, Y1_B) + + #swapaxes + tmp -= 2.0* contract('jkic,ib->jkcb', hbar.Hooov, Y1_B) + tmp += contract('kjic,ib->jkcb', hbar.Hooov, Y1_B) + + tmp += 2.0* contract('ajcb,ka->jkcb', hbar.Hvovv, Y1_B) + tmp -= contract('ajbc,ka->jkcb', hbar.Hvovv, Y1_B) + tmp += 2.0* contract('akbc,ja->jkcb', hbar.Hvovv, Y1_B) + tmp -= contract('akcb,ja->jkcb', hbar.Hvovv, Y1_B) + + tmp2 = contract('miae,me->ia', tmp, X1_A) + self.Bcon2 += contract('ia,ia->', tmp2, X1_C) + + # + tmp = -1.0* contract('janc,nkba->jckb', hbar.Hovov, Y2_B) + tmp -= contract('kanb,njca->jckb', hbar.Hovov, Y2_B) + tmp -= contract('jacn,nkab->jckb', hbar.Hovvo, Y2_B) + tmp -= contract('kabn,njac->jckb', hbar.Hovvo, Y2_B) + + # swapaxes? + tmp += 0.5*contract('fabc,jkfa->jckb', hbar.Hvvvv, Y2_B) + tmp += 0.5*contract('facb,kjfa->jckb', hbar.Hvvvv, Y2_B) + + # swapaxes? + tmp += 0.5* contract('kjin,nibc->jckb', hbar.Hoooo, Y2_B) + tmp += 0.5* contract('jkin,nicb->jckb', hbar.Hoooo, Y2_B) + tmp2 = contract('iema,me->ia', tmp, X1_A) + self.Bcon2 += contract('ia,ia->', tmp2, X1_C) + + tmp = contract('ijab,ijdb->ad', t2, Y2_B) + tmp = contract('ld,ad->la', X1_C, tmp) + tmp = contract('la,klca->kc', tmp, L[o,o,v,v]) + self.Bcon2 -= contract('kc,kc->', tmp, X1_A) + + tmp = contract('ijab,jlba->il', t2, Y2_B) + tmp2 = contract('kc,kicd->id', X1_A, L[o,o,v,v]) + tmp2 = contract('id,ld->il', tmp2, X1_C) + self.Bcon2 -= contract('il,il->', tmp2, tmp) + + tmp = contract('ijab,jkba->ik', t2, Y2_B) + tmp2 = contract('ld,lidc->ic', X1_C, L[o,o,v,v]) + tmp2 = contract('ic,kc->ik', tmp2, X1_A) + self.Bcon2 -= contract('ik,ik->', tmp2, tmp) + + tmp = contract('ijab,ijcb->ac', t2, Y2_B) + tmp = contract('kc,ac->ka', X1_A, tmp) + tmp2 = contract('ld,lkda->ka', X1_C, L[o,o,v,v]) + self.Bcon2 -= contract('ka,ka->', tmp2, tmp) + + # + tmp = contract("klcd,ijcd->ijkl", X2_C, Y2_B) + tmp = contract("ijkl,ijab->klab", tmp, X2_A) + self.Bcon2 += 0.5* contract('klab,klab->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,ikbd->jkad", X2_A, Y2_B) + tmp = contract("jkad,klcd->jlac", tmp, X2_C) + self.Bcon2 += contract('jlac,jlac->', tmp, ERI[o,o,v,v]) + + tmp = contract("klcd,ikdb->licb", X2_C, Y2_B) + tmp = contract("licb,ijab->ljca", tmp, X2_A) + self.Bcon2 += contract('ljca,ljac->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,klab->ijkl", X2_A, Y2_B) + tmp = contract("ijkl,klcd->ijcd", tmp, X2_C) + self.Bcon2 += 0.5* contract('ijcd,ijcd->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,ijac->bc", X2_A, L[o,o,v,v]) + tmp = contract("bc,klcd->klbd", tmp, X2_C) + self.Bcon2 -= contract("klbd,klbd->", tmp, Y2_B) + tmp = contract("ijab,ikab->jk", X2_A, L[o,o,v,v]) + tmp = contract("jk,klcd->jlcd", tmp, X2_C) + self.Bcon2 -= contract("jlcd,jlcd->", tmp, Y2_B) + tmp = contract("ikbc,klcd->ilbd", L[o,o,v,v], X2_C) + tmp = contract("ilbd,ijab->jlad", tmp, X2_A) + self.Bcon2 -= contract("jlad,jlad->", tmp, Y2_B) + tmp = contract("ijab,jlbc->ilac", X2_A, Y2_B) + tmp = contract("ilac,klcd->ikad", tmp, X2_C) + self.Bcon2 -= contract("ikad,ikad->", tmp, L[o,o,v,v]) + tmp = contract("klca,klcd->ad", L[o,o,v,v], X2_C) + tmp = contract("ad,ijdb->ijab", tmp, Y2_B) + self.Bcon2 -= contract("ijab,ijab->", tmp, X2_A) + tmp = contract("kicd,klcd->il", L[o,o,v,v], X2_C) + tmp = contract("ijab,il->ljab", X2_A, tmp) + self.Bcon2 -= contract("ljab,ljab->", tmp, Y2_B) + + tmp = contract("klcd,ikac->lida", X2_C, Y2_B) + tmp = contract("lida,jlbd->ijab", tmp, L[o,o,v,v]) + self.Bcon2 += 2.0* contract("ijab,ijab->", tmp, X2_A) + + # + #swapaxes + tmp = 2.0* contract("jkbc,kc->jb", X2_C, Y1_B) + tmp -= contract("jkcb,kc->jb", X2_C, Y1_B) + tmp = contract('ijab,jb->ia', L[o,o,v,v], tmp) + self.Bcon2 += contract("ia,ia->", tmp, X1_A) + + tmp = contract("jkbc,jkba->ca", X2_C, L[o,o,v,v]) + tmp = contract("ia,ca->ic", X1_A, tmp) + self.Bcon2 -= contract("ic,ic->", tmp, Y1_B) + + tmp = contract("jkbc,jibc->ki", X2_C, L[o,o,v,v]) + tmp = contract("ki,ia->ka", tmp, X1_A) + self.Bcon2 -= contract("ka,ka->", tmp, Y1_B) + + # + tmp = contract("klcd,lkdb->cb", X2_C, Y2_B) + tmp = contract("jb,cb->jc", X1_A, tmp) + self.Bcon2 -= contract("jc,jc->", tmp, hbar.Hov) + + tmp = contract("klcd,ljdc->kj", X2_C, Y2_B) + tmp = contract("kj,jb->kb", tmp, X1_A) + self.Bcon2 -= contract("kb,kb->", tmp, hbar.Hov) + + tmp = contract('lkda,klcd->ac', Y2_B, X2_C) + tmp2 = contract('jb,ajcb->ac', X1_A, hbar.Hvovv) + self.Bcon2 += 2.0* contract('ac,ac->', tmp, tmp2) + + tmp = contract('lkda,klcd->ac', Y2_B, X2_C) + tmp2 = contract('jb,ajbc->ac', X1_A, hbar.Hvovv) + self.Bcon2 -= contract('ac,ac->', tmp, tmp2) + + tmp = contract('jb,ljda->lbda', X1_A, Y2_B) + + #swapaxes + tmp2 = 2.0* contract('klcd,akbc->ldab', X2_C, hbar.Hvovv) + tmp2 -= contract('klcd,akcb->ldab', X2_C, hbar.Hvovv) + self.Bcon2 += contract('lbda,ldab->', tmp, tmp2) + + tmp = contract('ia,fkba->fkbi', X1_A, hbar.Hvovv) + tmp = contract('fkbi,jifc->kjbc', tmp, Y2_B) + self.Bcon2 -= contract('jkbc,kjbc->', X2_C, tmp) + + tmp = contract('ia,fjac->fjic', X1_A, hbar.Hvovv) + tmp = contract('fjic,ikfb->jkbc', tmp, Y2_B) + self.Bcon2 -= contract('jkbc,jkbc->', X2_C, tmp) + + tmp = contract('ia,jkfa->jkfi', X1_A, Y2_B) + tmp2 = contract('jkbc,fibc->jkfi', X2_C, hbar.Hvovv) + self.Bcon2 -= contract('jkfi,jkfi->', tmp2, tmp) + + tmp = contract('jb,kjib->ki', X1_A, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_C, Y2_B) + self.Bcon2 -= 2.0* contract('ki,ki->', tmp, tmp2) + + tmp = contract('jb,jkib->ki', X1_A, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_C, Y2_B) + self.Bcon2 += contract('ki,ki->', tmp, tmp2) + + #swapaxes + tmp = 2.0* contract('jkic,klcd->jild', hbar.Hooov, X2_C) + tmp -= contract('kjic,klcd->jild', hbar.Hooov, X2_C) + tmp = contract('jild,jb->bild', tmp, X1_A) + self.Bcon2 -= contract('bild,ilbd->', tmp, Y2_B) + + tmp = contract('ia,jkna->jkni', X1_A, hbar.Hooov) + tmp2 = contract('jkbc,nibc->jkni', X2_C, Y2_B) + self.Bcon2 += contract('jkni,jkni->', tmp2, tmp) + + tmp = contract('ia,nkab->nkib', X1_A, Y2_B) + tmp = contract('jkbc,nkib->jnic', X2_C, tmp) + self.Bcon2 += contract('jnic,ijnc->', tmp, hbar.Hooov) + + tmp = contract('ia,nkba->nkbi', X1_A, Y2_B) + tmp = contract('jkbc,nkbi->jnci', X2_C,tmp) + self.Bcon2 += contract('jnci,jinc->', tmp, hbar.Hooov) + + # + # swapaxes + tmp = 2.0* contract("jkbc,kc->jb", X2_A, Y1_B) + tmp -= contract("jkcb,kc->jb", X2_A, Y1_B) + tmp = contract('ijab,jb->ia', L[o,o,v,v], tmp) + self.Bcon2 += contract("ia,ia->", tmp, X1_C) + + tmp = contract("jkbc,jkba->ca", X2_A, L[o,o,v,v]) + tmp = contract("ia,ca->ic", X1_C, tmp) + self.Bcon2 -= contract("ic,ic->", tmp, Y1_B) + + tmp = contract("jkbc,jibc->ki", X2_A, L[o,o,v,v]) + tmp = contract("ki,ia->ka", tmp, X1_C) + self.Bcon2 -= contract("ka,ka->", tmp, Y1_B) + + # + tmp = contract("klcd,lkdb->cb", X2_A, Y2_B) + tmp = contract("jb,cb->jc", X1_C, tmp) + self.Bcon2 -= contract("jc,jc->", tmp, hbar.Hov) + + tmp = contract("klcd,ljdc->kj", X2_A, Y2_B) + tmp = contract("kj,jb->kb", tmp, X1_C) + self.Bcon2 -= contract("kb,kb->",tmp, hbar.Hov) + + tmp = contract('lkda,klcd->ac', Y2_B, X2_A) + tmp2 = contract('jb,ajcb->ac', X1_C, hbar.Hvovv) + self.Bcon2 += 2.0* contract('ac,ac->', tmp, tmp2) + + tmp = contract('lkda,klcd->ac', Y2_B, X2_A) + tmp2 = contract('jb,ajbc->ac', X1_C, hbar.Hvovv) + self.Bcon2 -= contract('ac,ac->', tmp, tmp2) + + tmp = contract('jb,ljda->lbda', X1_C, Y2_B) + + #swapaxes + tmp2 = 2.0* contract('klcd,akbc->ldab', X2_A, hbar.Hvovv) + tmp2 -= contract('klcd,akcb->ldab', X2_A, hbar.Hvovv) + self.Bcon2 += contract('lbda,ldab->', tmp, tmp2) + + tmp = contract('ia,fkba->fkbi', X1_C, hbar.Hvovv) + tmp = contract('fkbi,jifc->kjbc', tmp, Y2_B) + self.Bcon2 -= contract('jkbc,kjbc->', X2_A, tmp) + + tmp = contract('ia,fjac->fjic', X1_C, hbar.Hvovv) + tmp = contract('fjic,ikfb->jkbc', tmp, Y2_B) + self.Bcon2 -= contract('jkbc,jkbc->', X2_A, tmp) + + tmp = contract('ia,jkfa->jkfi', X1_C, Y2_B) + tmp2 = contract('jkbc,fibc->jkfi', X2_A, hbar.Hvovv) + self.Bcon2 -= contract('jkfi,jkfi->', tmp2, tmp) + + tmp = contract('jb,kjib->ki', X1_C, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_A, Y2_B) + self.Bcon2 -= 2.0* contract('ki,ki->', tmp, tmp2) + + tmp = contract('jb,jkib->ki', X1_C, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_A, Y2_B) + self.Bcon2 += 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_C) + self.Bcon2 -= contract('bild,ilbd->', tmp, Y2_B) + + tmp = contract('ia,jkna->jkni', X1_C, hbar.Hooov) + tmp2 = contract('jkbc,nibc->jkni', X2_A, Y2_B) + self.Bcon2 += contract('jkni,jkni->', tmp2, tmp) + + tmp = contract('ia,nkab->nkib', X1_C, Y2_B) + tmp = contract('jkbc,nkib->jnic', X2_A, tmp) + self.Bcon2 += contract('jnic,ijnc->', tmp, hbar.Hooov) + + tmp = contract('ia,nkba->nkbi', X1_C, Y2_B) + tmp = contract('jkbc,nkbi->jnci', X2_A, tmp) + self.Bcon2 += contract('jnci,jinc->', tmp, hbar.Hooov) + + self.Bcon3 = 0 + # <0|L1(C)[[Hbar(0),X1(A),X1(B)]]|0> + tmp = -1.0* contract('jc,kb->jkcb', hbar.Hov, Y1_C) + tmp -= contract('jc,kb->jkcb', Y1_C, hbar.Hov) + + #swapaxes + tmp -= 2.0* contract('kjib,ic->jkcb', hbar.Hooov, Y1_C) + tmp += contract('jkib,ic->jkcb', hbar.Hooov, Y1_C) + + #swapaxes + tmp -= 2.0* contract('jkic,ib->jkcb', hbar.Hooov, Y1_C) + tmp += contract('kjic,ib->jkcb', hbar.Hooov, Y1_C) + + #swapaxes + tmp += 2.0* contract('ajcb,ka->jkcb', hbar.Hvovv, Y1_C) + tmp -= contract('ajbc,ka->jkcb', hbar.Hvovv, Y1_C) + tmp += 2.0* contract('akbc,ja->jkcb', hbar.Hvovv, Y1_C) + tmp -= contract('akcb,ja->jkcb', hbar.Hvovv, Y1_C) + + tmp2 = contract('miae,me->ia', tmp, X1_A) + self.Bcon3 += contract('ia,ia->', tmp2, X1_B) + + # <0|L2(C)|[[Hbar(0),X1(A)],X1(B)]|0> + tmp = -1.0* contract('janc,nkba->jckb', hbar.Hovov, Y2_C) + tmp -= contract('kanb,njca->jckb', hbar.Hovov, Y2_C) + tmp -= contract('jacn,nkab->jckb', hbar.Hovvo, Y2_C) + tmp -= contract('kabn,njac->jckb', hbar.Hovvo, Y2_C) + + #swapaxes? + tmp += 0.5* contract('fabc,jkfa->jckb', hbar.Hvvvv, Y2_C) + tmp += 0.5* contract('facb,kjfa->jckb', hbar.Hvvvv, Y2_C) + + #swapaxes? + tmp += 0.5* contract('kjin,nibc->jckb', hbar.Hoooo, Y2_C) + tmp += 0.5* contract('jkin,nicb->jckb', hbar.Hoooo, Y2_C) + tmp2 = contract('iema,me->ia', tmp, X1_A) + self.Bcon3 += contract('ia,ia->', tmp2, X1_B) + + tmp = contract('ijab,ijdb->ad', t2, Y2_C) + tmp = contract('ld,ad->la', X1_B, tmp) + tmp = contract('la,klca->kc', tmp, L[o,o,v,v]) + self.Bcon3 -= contract('kc,kc->',tmp, X1_A) + + tmp = contract('ijab,jlba->il', t2, Y2_C) + tmp2 = contract('kc,kicd->id', X1_A, L[o,o,v,v]) + tmp2 = contract('id,ld->il', tmp2, X1_B) + self.Bcon3 -= contract('il,il->', tmp2, tmp) + + tmp = contract('ijab,jkba->ik', t2, Y2_C) + tmp2 = contract('ld,lidc->ic', X1_B, L[o,o,v,v]) + tmp2 = contract('ic,kc->ik', tmp2, X1_A) + self.Bcon3 -= contract('ik,ik->', tmp2, tmp) + + tmp = contract('ijab,ijcb->ac', t2, Y2_C) + tmp = contract('kc,ac->ka', X1_A, tmp) + tmp2 = contract('ld,lkda->ka', X1_B, L[o,o,v,v]) + self.Bcon3 -= contract('ka,ka->', tmp2, tmp) + + # <0|L2(C)[[Hbar(0),X2(A)],X2(B)]|0> + tmp = contract("klcd,ijcd->ijkl", X2_B, Y2_C) + tmp = contract("ijkl,ijab->klab", tmp, X2_A) + self.Bcon3 += 0.5* contract('klab,klab->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,ikbd->jkad", X2_A, Y2_C) + tmp = contract("jkad,klcd->jlac", tmp, X2_B) + self.Bcon3 += contract('jlac,jlac->', tmp, ERI[o,o,v,v]) + + tmp = contract("klcd,ikdb->licb", X2_B, Y2_C) + tmp = contract("licb,ijab->ljca", tmp, X2_A) + self.Bcon3 += contract('ljca,ljac->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,klab->ijkl", X2_A, Y2_C) + tmp = contract("ijkl,klcd->ijcd", tmp, X2_B) + self.Bcon3 += 0.5* contract('ijcd,ijcd->', tmp, ERI[o,o,v,v]) + + tmp = contract("ijab,ijac->bc", X2_A, L[o,o,v,v]) + tmp = contract("bc,klcd->klbd", tmp, X2_B) + self.Bcon3 -= contract("klbd,klbd->", tmp, Y2_C) + tmp = contract("ijab,ikab->jk", X2_A, L[o,o,v,v]) + tmp = contract("jk,klcd->jlcd", tmp, X2_B) + self.Bcon3 -= contract("jlcd,jlcd->", tmp, Y2_C) + tmp = contract("ikbc,klcd->ilbd", L[o,o,v,v], X2_B) + tmp = contract("ilbd,ijab->jlad", tmp, X2_A) + self.Bcon3 -= contract("jlad,jlad->", tmp, Y2_C) + tmp = contract("ijab,jlbc->ilac", X2_A, Y2_C) + tmp = contract("ilac,klcd->ikad", tmp, X2_B) + self.Bcon3 -= contract("ikad,ikad->", tmp, L[o,o,v,v]) + tmp = contract("klca,klcd->ad", L[o,o,v,v], X2_B) + tmp = contract("ad,ijdb->ijab", tmp, Y2_C) + self.Bcon3 -= contract("ijab,ijab->", tmp, X2_A) + tmp = contract("kicd,klcd->il", L[o,o,v,v], X2_B) + tmp = contract("ijab,il->ljab", X2_A, tmp) + self.Bcon3 -= contract("ljab,ljab->", tmp, Y2_C) + + tmp = contract("klcd,ikac->lida", X2_B, Y2_C) + tmp = contract("lida,jlbd->ijab", tmp, L[o,o,v,v]) + self.Bcon3 += 2.0* contract("ijab,ijab->", tmp, X2_A) + + # <0|L1(C)[[Hbar(0),X1(A)],X2(B)]]|0> + #swapaxes + tmp = 2.0 * contract("jkbc,kc->jb", X2_B, Y1_C) + tmp -= contract("jkcb,kc->jb", X2_B, Y1_C) + tmp = contract('ijab,jb->ia', L[o,o,v,v], tmp) + self.Bcon3 += contract("ia,ia->", tmp, X1_A) + + tmp = contract("jkbc,jkba->ca", X2_B, L[o,o,v,v]) + tmp = contract("ia,ca->ic", X1_A, tmp) + self.Bcon3 -= contract("ic,ic->", tmp, Y1_C) + + tmp = contract("jkbc,jibc->ki", X2_B, L[o,o,v,v]) + tmp = contract("ki,ia->ka", tmp, X1_A) + self.Bcon3 -= contract("ka,ka->", tmp, Y1_C) + + # <0|L2(C)[[Hbar(0),X1(A)],X2(B)]]|0> + tmp = contract("klcd,lkdb->cb", X2_B, Y2_C) + tmp = contract("jb,cb->jc", X1_A, tmp) + self.Bcon3 -= contract("jc,jc->", tmp, hbar.Hov) + + tmp = contract("klcd,ljdc->kj", X2_B, Y2_C) + tmp = contract("kj,jb->kb", tmp, X1_A) + self.Bcon3 -= contract("kb,kb->", tmp, hbar.Hov) + + tmp = contract('lkda,klcd->ac', Y2_C, X2_B) + tmp2 = contract('jb,ajcb->ac', X1_A, hbar.Hvovv) + self.Bcon3 += 2.0* contract('ac,ac->', tmp, tmp2) + + tmp = contract('lkda,klcd->ac', Y2_C, X2_B) + tmp2 = contract('jb,ajbc->ac', X1_A, hbar.Hvovv) + self.Bcon3 -= contract('ac,ac->', tmp, tmp2) + + tmp = contract('jb,ljda->lbda', X1_A, Y2_C) + + #swapaxes + tmp2 = 2.0* contract('klcd,akbc->ldab', X2_B, hbar.Hvovv) + tmp2 -= contract('klcd,akcb->ldab', X2_B, hbar.Hvovv) + self.Bcon3 += contract('lbda,ldab->', tmp, tmp2) + + tmp = contract('ia,fkba->fkbi', X1_A, hbar.Hvovv) + tmp = contract('fkbi,jifc->kjbc', tmp, Y2_C) + self.Bcon3 -= contract('jkbc,kjbc->', X2_B, tmp) + + tmp = contract('ia,fjac->fjic', X1_A, hbar.Hvovv) + tmp = contract('fjic,ikfb->jkbc', tmp, Y2_C) + self.Bcon3 -= contract('jkbc,jkbc->', X2_B, tmp) + + tmp = contract('ia,jkfa->jkfi', X1_A, Y2_C) + tmp2 = contract('jkbc,fibc->jkfi', X2_B, hbar.Hvovv) + self.Bcon3 -= contract('jkfi,jkfi->', tmp2, tmp) + + tmp = contract('jb,kjib->ki', X1_A, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_B, Y2_C) + self.Bcon3 -= 2.0* contract('ki,ki->', tmp, tmp2) + + tmp = contract('jb,jkib->ki', X1_A, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_B, Y2_C) + self.Bcon3 += contract('ki,ki->', tmp, tmp2) + + #swapaxes + 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) + self.Bcon3 -= contract('bild,ilbd->', tmp, Y2_C) + + tmp = contract('ia,jkna->jkni', X1_A, hbar.Hooov) + tmp2 = contract('jkbc,nibc->jkni', X2_B, Y2_C) + self.Bcon3 += contract('jkni,jkni->', tmp2, tmp) + + tmp = contract('ia,nkab->nkib', X1_A, Y2_C) + tmp = contract('jkbc,nkib->jnic', X2_B, tmp) + self.Bcon3 += contract('jnic,ijnc->', tmp, hbar.Hooov) + + tmp = contract('ia,nkba->nkbi', X1_A, Y2_C) + tmp = contract('jkbc,nkbi->jnci', X2_B, tmp) + self.Bcon3 += contract('jnci,jinc->', tmp, hbar.Hooov) + + # <0|L1(C)[[Hbar(0),X2(A)],X1(B)]]|0> + tmp = 2.0* contract("jkbc,kc->jb", X2_A, Y1_C) + tmp -= contract("jkcb,kc->jb", X2_A, Y1_C) + tmp = contract('ijab,jb->ia', L[o,o,v,v], tmp) + self.Bcon3 += contract("ia,ia->", tmp, X1_B) + + tmp = contract("jkbc,jkba->ca", X2_A, L[o,o,v,v]) + tmp = contract("ia,ca->ic", X1_B, tmp) + self.Bcon3 -= contract("ic,ic->", tmp, Y1_C) + + tmp = contract("jkbc,jibc->ki", X2_A, L[o,o,v,v]) + tmp = contract("ki,ia->ka", tmp, X1_B) + self.Bcon3 -= contract("ka,ka->", tmp, Y1_C) + + # <0|L1(C)[[Hbar(0),X2(A)],X1(B)]]|0> + tmp = contract("klcd,lkdb->cb", X2_A, Y2_C) + tmp = contract("jb,cb->jc", X1_B, tmp) + self.Bcon3 -= contract("jc,jc->", tmp, hbar.Hov) + + tmp = contract("klcd,ljdc->kj", X2_A, Y2_C) + tmp = contract("kj,jb->kb", tmp, X1_B) + self.Bcon3 -= contract("kb,kb->",tmp, hbar.Hov) + + tmp = contract('lkda,klcd->ac', Y2_C, X2_A) + tmp2 = contract('jb,ajcb->ac', X1_B, hbar.Hvovv) + self.Bcon3 += 2.0* contract('ac,ac->', tmp, tmp2) + + tmp = contract('lkda,klcd->ac', Y2_C, X2_A) + tmp2 = contract('jb,ajbc->ac', X1_B, hbar.Hvovv) + self.Bcon3 -= contract('ac,ac->', tmp, tmp2) + + tmp = contract('jb,ljda->lbda', X1_B, Y2_C) + + #swapaxes + tmp2 = 2.0* contract('klcd,akbc->ldab', X2_A, hbar.Hvovv) + tmp2 -= contract('klcd,akcb->ldab', X2_A, hbar.Hvovv) + self.Bcon3 += contract('lbda,ldab->', tmp, tmp2) + + tmp = contract('ia,fkba->fkbi', X1_B, hbar.Hvovv) + tmp = contract('fkbi,jifc->kjbc', tmp, Y2_C) + self.Bcon3 -= contract('jkbc,kjbc->', X2_A, tmp) + + tmp = contract('ia,fjac->fjic', X1_B, hbar.Hvovv) + tmp = contract('fjic,ikfb->jkbc', tmp, Y2_C) + self.Bcon3 -= contract('jkbc,jkbc->', X2_A, tmp) + + tmp = contract('ia,jkfa->jkfi', X1_B, Y2_C) + tmp2 = contract('jkbc,fibc->jkfi', X2_A, hbar.Hvovv) + self.Bcon3 -= contract('jkfi,jkfi->', tmp2, tmp) + + tmp = contract('jb,kjib->ki', X1_B, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_A, Y2_C) + self.Bcon3 -= 2.0*contract('ki,ki->', tmp, tmp2) + + tmp = contract('jb,jkib->ki', X1_B, hbar.Hooov) + tmp2 = contract('klcd,ilcd->ki', X2_A, Y2_C) + self.Bcon3 += contract('ki,ki->', tmp, tmp2) + + #swapaxes + 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) + self.Bcon3 -= contract('bild,ilbd->', tmp, Y2_C) + + tmp = contract('ia,jkna->jkni', X1_B, hbar.Hooov) + tmp2 = contract('jkbc,nibc->jkni', X2_A, Y2_C) + self.Bcon3 += contract('jkni,jkni->', tmp2, tmp) + + tmp = contract('ia,nkab->nkib', X1_B, Y2_C) + tmp = contract('jkbc,nkib->jnic', X2_A, tmp) + self.Bcon3 += contract('jnic,ijnc->', tmp, hbar.Hooov) + + tmp = contract('ia,nkba->nkbi', X1_B, Y2_C) + tmp = contract('jkbc,nkbi->jnci', X2_A, tmp) + self.Bcon3 += contract('jnci,jinc->', tmp, hbar.Hooov) + + self.hyper += self.Bcon1 + self.Bcon2 + self.Bcon3 + + return self.hyper + + def hyperpolar(self): + """ + Return + ------ + Beta_avg: float + Hyperpolarizability average + """ + solver_start = time.time() + + ccpert_om1_X = self.ccpert_om1_X + ccpert_om2_X = self.ccpert_om2_X + ccpert_om_sum_X = self.ccpert_om_sum_X + + ccpert_om1_2nd_X = self.ccpert_om1_2nd_X + ccpert_om2_2nd_X = self.ccpert_om2_2nd_X + ccpert_om_sum_2nd_X = self.ccpert_om_sum_2nd_X + + ccpert_om1_Y = self.ccpert_om1_Y + ccpert_om2_Y = self.ccpert_om2_Y + ccpert_om_sum_Y = self.ccpert_om_sum_Y + + ccpert_om1_2nd_Y = self.ccpert_om1_2nd_Y + ccpert_om2_2nd_Y = self.ccpert_om2_2nd_Y + ccpert_om_sum_2nd_Y = self.ccpert_om_sum_2nd_Y + + hyper_AB_1st = np.zeros((3,3,3)) + hyper_AB_2nd = np.zeros((3,3,3)) + hyper_AB = np.zeros((3,3,3)) + + for a in range(0, 3): + pertkey_a = "MU_" + self.cart[a] + for b in range(0, 3): + pertkey_b = "MU_" + self.cart[b] + for c in range(0, 3): + pertkey_c = "MU_" + self.cart[c] + + hyper_AB_1st[a,b,c] = self.quadraticresp(pertkey_a, pertkey_b, pertkey_c, ccpert_om_sum_X[pertkey_a], ccpert_om1_X[pertkey_b], ccpert_om2_X[pertkey_c], ccpert_om_sum_Y[pertkey_a], ccpert_om1_Y[pertkey_b], ccpert_om2_Y[pertkey_c] ) + hyper_AB_2nd[a,b,c] = self.quadraticresp(pertkey_a, pertkey_b, pertkey_c, ccpert_om_sum_2nd_X[pertkey_a], ccpert_om1_2nd_X[pertkey_b], ccpert_om2_2nd_X[pertkey_c], ccpert_om_sum_2nd_Y[pertkey_a], ccpert_om1_2nd_Y[pertkey_b], ccpert_om2_2nd_Y[pertkey_c]) + hyper_AB[a,b,c] = (hyper_AB_1st[a,b,c] + hyper_AB_2nd[a,b,c] )/2 + + Beta_avg = 0 + for i in range(0,3): + Beta_avg += (hyper_AB[2,i,i] + hyper_AB[i,2,i] + hyper_AB[i,i,2])/5 + + print("Beta_avg = %10.12lf" %(Beta_avg)) + print("\n First Dipole Hyperpolarizability computed in %.3f seconds.\n" % (time.time() - solver_start)) + + return Beta_avg + + def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, max_diis=7, start_diis=1): solver_start = time.time() Dia = self.Dia @@ -319,12 +1812,21 @@ def solve_right(self, pertbar, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, m r1 = self.r_X1(pertbar, omega) r2 = self.r_X2(pertbar, omega) - self.X1 += r1/(Dia + omega) - self.X2 += r2/(Dijab + omega) + if self.ccwfn.local is not None: + inc1, inc2 = self.ccwfn.Local.filter_pertamps(r1, r2, self.eps_occ, self.eps_vir, omega) + self.X1 += inc1 + self.X2 += inc2 + + rms = contract('ia,ia->', np.conj(inc1/(Dia+omega)), inc1/(Dia+omega)) + rms += contract('ijab,ijab->', np.conj(inc2/(Dijab+omega)), inc2/(Dijab+omega)) + rms = np.sqrt(rms) + else: + self.X1 += r1/(Dia + omega) + self.X2 += r2/(Dijab + omega) - rms = contract('ia,ia->', np.conj(r1/(Dia+omega)), r1/(Dia+omega)) - rms += contract('ijab,ijab->', np.conj(r2/(Dijab+omega)), r2/(Dijab+omega)) - rms = np.sqrt(rms) + rms = contract('ia,ia->', np.conj(r1/(Dia+omega)), r1/(Dia+omega)) + rms += contract('ijab,ijab->', np.conj(r2/(Dijab+omega)), r2/(Dijab+omega)) + rms = np.sqrt(rms) pseudo = self.pseudoresponse(pertbar, self.X1, self.X2) pseudodiff = np.abs(pseudo - pseudo_last) @@ -338,6 +1840,78 @@ def solve_right(self, pertbar, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, m if niter >= start_diis: self.X1, self.X2 = diis.extrapolate(self.X1, self.X2) + def solve_left(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, max_diis=7, start_diis=1): + ''' + Notes + ----- + The first-order lambda equations are partition into two expressions: inhomogeneous (in_Y1 and in_Y2) and homogeneous terms (r_Y1 and r_Y2), + the inhomogenous terms contains only terms that are not changing over the iterative process of obtaining the solutions for these equations. Therefore, it is + computed only once and is called when solving for the homogenous terms. + ''' + solver_start = time.time() + + Dia = self.Dia + Dijab = self.Dijab + + # initial guess + X1_guess = pertbar.Avo.T/(Dia + omega) + X2_guess = pertbar.Avvoo/(Dijab + omega) + + # initial guess + Y1 = 2.0 * X1_guess.copy() + Y2 = 4.0 * X2_guess.copy() + Y2 -= 2.0 * X2_guess.copy().swapaxes(2,3) + + pseudo = self.pseudoresponse(pertbar, Y1, Y2) + print(f"Iter {0:3d}: CC Pseudoresponse = {pseudo.real:.15f} dP = {pseudo.real:.5E}") + + diis = helper_diis(Y1, Y2, max_diis) + contract = self.ccwfn.contract + + self.Y1 = Y1 + self.Y2 = Y2 + + # uses updated X1 and X2 + self.im_Y1 = self.in_Y1(pertbar, self.X1, self.X2) + self.im_Y2 = self.in_Y2(pertbar, self.X1, self.X2) + + for niter in range(1, maxiter+1): + pseudo_last = pseudo + + Y1 = self.Y1 + Y2 = self.Y2 + + r1 = self.r_Y1(pertbar, omega) + r2 = self.r_Y2(pertbar, omega) + + if self.ccwfn.local is not None: + inc1, inc2 = self.ccwfn.Local.filter_pertamps(r1, r2, self.eps_occ, self.eps_vir, omega) + self.Y1 += inc1 + self.Y2 += inc2 + + rms = contract('ia,ia->', np.conj(inc1/(Dia+omega)), inc1/(Dia+omega)) + rms += contract('ijab,ijab->', np.conj(inc2/(Dijab+omega)), inc2/(Dijab+omega)) + rms = np.sqrt(rms) + else: + self.Y1 += r1/(Dia + omega) + self.Y2 += r2/(Dijab + omega) + + rms = contract('ia,ia->', np.conj(r1/(Dia+omega)), r1/(Dia+omega)) + rms += contract('ijab,ijab->', np.conj(r2/(Dijab+omega)), r2/(Dijab+omega)) + rms = np.sqrt(rms) + + pseudo = self.pseudoresponse(pertbar, self.Y1, self.Y2) + pseudodiff = np.abs(pseudo - pseudo_last) + print(f"Iter {niter:3d}: CC Pseudoresponse = {pseudo.real:.15f} dP = {pseudodiff:.5E} rms = {rms.real:.5E}") + + if ((abs(pseudodiff) < e_conv) and abs(rms) < r_conv): + print("\nPerturbed wave function converged in %.3f seconds.\n" % (time.time() - solver_start)) + return self.Y1, self.Y2 , pseudo + + diis.add_error_vector(self.Y1, self.Y2) + if niter >= start_diis: + self.Y1, self.Y2 = diis.extrapolate(self.Y1, self.Y2) + def r_X1(self, pertbar, omega): contract = self.contract o = self.ccwfn.o @@ -391,14 +1965,304 @@ def r_X2(self, pertbar, omega): return r_X2 + def in_Y1(self, pertbar, X1, X2): + contract = self.contract + o = self.ccwfn.o + v = self.ccwfn.v + Y1 = self.Y1 + + Y2 = self.Y2 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + cclambda = self.cclambda + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.H.L + + # good + r_Y1 = 2.0 * pertbar.Aov.copy() + # good + r_Y1 -= contract('im,ma->ia', pertbar.Aoo, l1) + r_Y1 += contract('ie,ea->ia', l1, pertbar.Avv) + # + r_Y1 += contract('imfe,feam->ia', l2, pertbar.Avvvo) + + #can combine the next two to swapaxes type contraction + r_Y1 -= 0.5 * contract('ienm,mnea->ia', pertbar.Aovoo, l2) + r_Y1 -= 0.5 * contract('iemn,mnae->ia', pertbar.Aovoo, l2) + + # good + r_Y1 += 2.0 * contract('imae,me->ia', L[o,o,v,v], X1) + + # + tmp = -1.0 * contract('ma,ie->miae', hbar.Hov, l1) + tmp -= contract('ma,ie->miae', l1, hbar.Hov) + tmp -= 2.0 * contract('mina,ne->miae', hbar.Hooov, l1) + + #double check this one + tmp += contract('imna,ne->miae', hbar.Hooov, l1) + + #can combine the next two to swapaxes type contraction + tmp -= 2.0 * contract('imne,na->miae', hbar.Hooov, l1) + tmp += contract('mine,na->miae', hbar.Hooov, l1) + + #can combine the next two to swapaxes type contraction + tmp += 2.0 * contract('fmae,if->miae', hbar.Hvovv, l1) + tmp -= contract('fmea,if->miae', hbar.Hvovv, l1) + + #can combine the next two to swapaxes type contraction + tmp += 2.0 * contract('fiea,mf->miae', hbar.Hvovv, l1) + tmp -= contract('fiae,mf->miae', hbar.Hvovv, l1) + r_Y1 += contract('miae,me->ia', tmp, X1) + + # good + + #can combine the next two to swapaxes type contraction + tmp = 2.0 * contract('mnef,nf->me', X2, l1) + tmp -= contract('mnfe,nf->me', X2, l1) + r_Y1 += contract('imae,me->ia', L[o,o,v,v], tmp) + r_Y1 -= contract('ni,na->ia', cclambda.build_Goo(X2, L[o,o,v,v]), l1) + r_Y1 += contract('ie,ea->ia', l1, cclambda.build_Gvv(L[o,o,v,v], X2)) + + # good + + # can reorganize thesenext four to two swapaxes type contraction + tmp = -1.0 * contract('nief,mfna->iema', l2, hbar.Hovov) + tmp -= contract('ifne,nmaf->iema', hbar.Hovov, l2) + tmp -= contract('inef,mfan->iema', l2, hbar.Hovvo) + tmp -= contract('ifen,nmfa->iema', hbar.Hovvo, l2) + + #can combine the next two to swapaxes type contraction + tmp += 0.5 * contract('imfg,fgae->iema', l2, hbar.Hvvvv) + tmp += 0.5 * contract('imgf,fgea->iema', l2, hbar.Hvvvv) + + #can combine the next two to swapaxes type contraction + tmp += 0.5 * contract('imno,onea->iema', hbar.Hoooo, l2) + tmp += 0.5 * contract('mino,noea->iema', hbar.Hoooo, l2) + r_Y1 += contract('iema,me->ia', tmp, X1) + + tmp = contract('nb,fb->nf', X1, cclambda.build_Gvv(l2, t2)) + r_Y1 += contract('inaf,nf->ia', L[o,o,v,v], tmp) + tmp = contract('me,fa->mefa', X1, cclambda.build_Gvv(l2, t2)) + r_Y1 += contract('mief,mefa->ia', L[o,o,v,v], tmp) + tmp = contract('me,ni->meni', X1, cclambda.build_Goo(t2, l2)) + r_Y1 -= contract('meni,mnea->ia', tmp, L[o,o,v,v]) + tmp = contract('jf,nj->fn', X1, cclambda.build_Goo(t2, l2)) + r_Y1 -= contract('inaf,fn->ia', L[o,o,v,v], tmp) + + # + r_Y1 -= contract('mi,ma->ia', cclambda.build_Goo(X2, l2), hbar.Hov) + r_Y1 += contract('ie,ea->ia', hbar.Hov, cclambda.build_Gvv(l2, X2)) + tmp = contract('imfg,mnef->igne', l2, X2) + r_Y1 -= contract('igne,gnea->ia', tmp, hbar.Hvovv) + tmp = contract('mifg,mnef->igne', l2, X2) + r_Y1 -= contract('igne,gnae->ia', tmp, hbar.Hvovv) + tmp = contract('mnga,mnef->gaef', l2, X2) + r_Y1 -= contract('gief,gaef->ia', hbar.Hvovv, tmp) + + #can combine the next two to swapaxes type contraction + tmp = 2.0 * contract('gmae,mnef->ganf', hbar.Hvovv, X2) + tmp -= contract('gmea,mnef->ganf', hbar.Hvovv, X2) + r_Y1 += contract('nifg,ganf->ia', l2, tmp) + + #can combine the next two to swapaxes type contraction + r_Y1 -= 2.0 * contract('giea,ge->ia', hbar.Hvovv, cclambda.build_Gvv(X2, l2)) + r_Y1 += contract('giae,ge->ia', hbar.Hvovv, cclambda.build_Gvv(X2, l2)) + tmp = contract('oief,mnef->oimn', l2, X2) + r_Y1 += contract('oimn,mnoa->ia', tmp, hbar.Hooov) + tmp = contract('mofa,mnef->oane', l2, X2) + r_Y1 += contract('inoe,oane->ia', hbar.Hooov, tmp) + tmp = contract('onea,mnef->oamf', l2, X2) + r_Y1 += contract('miof,oamf->ia', hbar.Hooov, tmp) + + #can combine the next two to swapaxes type contraction + r_Y1 -= 2.0 * contract('mioa,mo->ia', hbar.Hooov, cclambda.build_Goo(X2, l2)) + r_Y1 += contract('imoa,mo->ia', hbar.Hooov, cclambda.build_Goo(X2, l2)) + + #can combine the next two to swapaxes type contraction + tmp = -2.0 * contract('imoe,mnef->ionf', hbar.Hooov, X2) + tmp += contract('mioe,mnef->ionf', hbar.Hooov, X2) + r_Y1 += contract('ionf,nofa->ia', tmp, l2) + + return r_Y1 + + def r_Y1(self, pertbar, omega): + contract = self.contract + o = self.ccwfn.o + v = self.ccwfn.v + Y1 = self.Y1 + Y2 = self.Y2 + l2 = self.cclambda.l2 + cclambda = self.cclambda + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.H.L + + #imhomogenous terms + r_Y1 = self.im_Y1.copy() + r_Y1 += omega * Y1 + r_Y1 += contract('ie,ea->ia', Y1, hbar.Hvv) + r_Y1 -= contract('im,ma->ia', hbar.Hoo, Y1) + r_Y1 += 2.0 * contract('ieam,me->ia', hbar.Hovvo, Y1) + r_Y1 -= contract('iema,me->ia', hbar.Hovov, Y1) + r_Y1 += contract('imef,efam->ia', Y2, hbar.Hvvvo) + r_Y1 -= contract('iemn,mnae->ia', hbar.Hovoo, Y2) + + #can combine the next two to swapaxes type contraction + r_Y1 -= 2.0 * contract('eifa,ef->ia', hbar.Hvovv, cclambda.build_Gvv(t2, Y2)) + r_Y1 += contract('eiaf,ef->ia', hbar.Hvovv, cclambda.build_Gvv(t2, Y2)) + + #can combine the next two to swapaxes type contraction + r_Y1 -= 2.0 * contract('mina,mn->ia', hbar.Hooov, cclambda.build_Goo(t2, Y2)) + r_Y1 += contract('imna,mn->ia', hbar.Hooov, cclambda.build_Goo(t2, Y2)) + + return r_Y1 + + def in_Y2(self, pertbar, X1, X2): + contract = self.contract + o = self.ccwfn.o + v = self.ccwfn.v + #X1 = self.X1 + #X2 = self.X2 + Y1 = self.Y1 + Y2 = self.Y2 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + cclambda = self.cclambda + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.H.L + ERI = self.H.ERI + + # good + + #next two turn to swapaxes contraction + r_Y2 = 2.0 * contract('ia,jb->ijab', l1, pertbar.Aov.copy()) + r_Y2 -= contract('ja,ib->ijab', l1, pertbar.Aov) + + # good + r_Y2 += contract('ijeb,ea->ijab', l2, pertbar.Avv) + r_Y2 -= contract('im,mjab->ijab', pertbar.Aoo, l2) + + # good + tmp = contract('me,ja->meja', X1, l1) + r_Y2 -= contract('mieb,meja->ijab', L[o,o,v,v], tmp) + tmp = contract('me,mb->eb', X1, l1) + r_Y2 -= contract('ijae,eb->ijab', L[o,o,v,v], tmp) + tmp = contract('me,ie->mi', X1, l1) + r_Y2 -= contract('mi,jmba->ijab', tmp, L[o,o,v,v]) + tmp = 2.0 *contract('me,jb->mejb', X1, l1) + r_Y2 += contract('imae,mejb->ijab', L[o,o,v,v], tmp) + + # + tmp = contract('me,ma->ea', X1, hbar.Hov) + r_Y2 -= contract('ijeb,ea->ijab', l2, tmp) + tmp = contract('me,ie->mi', X1, hbar.Hov) + r_Y2 -= contract('mi,jmba->ijab', tmp, l2) + tmp = contract('me,ijef->mijf', X1, l2) + r_Y2 -= contract('mijf,fmba->ijab', tmp, hbar.Hvovv) + tmp = contract('me,imbf->eibf', X1, l2) + r_Y2 -= contract('eibf,fjea->ijab', tmp, hbar.Hvovv) + tmp = contract('me,jmfa->ejfa', X1, l2) + r_Y2 -= contract('fibe,ejfa->ijab', hbar.Hvovv, tmp) + + #swapaxes contraction + tmp = 2.0 * contract('me,fmae->fa', X1, hbar.Hvovv) + tmp -= contract('me,fmea->fa', X1, hbar.Hvovv) + r_Y2 += contract('ijfb,fa->ijab', l2, tmp) + + #swapaxes contraction + tmp = 2.0 * contract('me,fiea->mfia', X1, hbar.Hvovv) + tmp -= contract('me,fiae->mfia', X1, hbar.Hvovv) + r_Y2 += contract('mfia,jmbf->ijab', tmp, l2) + tmp = contract('me,jmna->ejna', X1, hbar.Hooov) + r_Y2 += contract('ineb,ejna->ijab', l2, tmp) + + tmp = contract('me,mjna->ejna', X1, hbar.Hooov) + r_Y2 += contract('nieb,ejna->ijab', l2, tmp) + tmp = contract('me,nmba->enba', X1, l2) + r_Y2 += contract('jine,enba->ijab', hbar.Hooov, tmp) + + #swapaxes + tmp = 2.0 * contract('me,mina->eina', X1, hbar.Hooov) + tmp -= contract('me,imna->eina', X1, hbar.Hooov) + r_Y2 -= contract('eina,njeb->ijab', tmp, l2) + + #swapaxes + tmp = 2.0 * contract('me,imne->in', X1, hbar.Hooov) + tmp -= contract('me,mine->in', X1, hbar.Hooov) + r_Y2 -= contract('in,jnba->ijab', tmp, l2) + + # + tmp = 0.5 * contract('ijef,mnef->ijmn', l2, X2) + r_Y2 += contract('ijmn,mnab->ijab', tmp, ERI[o,o,v,v]) + tmp = 0.5 * contract('ijfe,mnef->ijmn', ERI[o,o,v,v], X2) + r_Y2 += contract('ijmn,mnba->ijab', tmp, l2) + tmp = contract('mifb,mnef->ibne', l2, X2) + r_Y2 += contract('ibne,jnae->ijab', tmp, ERI[o,o,v,v]) + tmp = contract('imfb,mnef->ibne', l2, X2) + r_Y2 += contract('ibne,njae->ijab', tmp, ERI[o,o,v,v]) + tmp = contract('mjfb,mnef->jbne', l2, X2) + r_Y2 -= contract('jbne,inae->ijab', tmp, L[o,o,v,v]) + + #temp intermediate? + r_Y2 -= contract('in,jnba->ijab', cclambda.build_Goo(L[o,o,v,v], X2), l2) + r_Y2 += contract('ijfb,af->ijab', l2, cclambda.build_Gvv(X2, L[o,o,v,v])) + r_Y2 += contract('ijae,be->ijab', L[o,o,v,v], cclambda.build_Gvv(X2, l2)) + r_Y2 -= contract('imab,jm->ijab', L[o,o,v,v], cclambda.build_Goo(l2, X2)) + tmp = contract('nifb,mnef->ibme', l2, X2) + r_Y2 -= contract('ibme,mjea->ijab', tmp, L[o,o,v,v]) + tmp = 2.0 * contract('njfb,mnef->jbme', l2, X2) + r_Y2 += contract('imae,jbme->ijab', L[o,o,v,v], tmp) + + return r_Y2 + + def r_Y2(self, pertbar, omega): + contract = self.contract + o = self.ccwfn.o + v = self.ccwfn.v + Y1 = self.Y1 + Y2 = self.Y2 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + cclambda = self.cclambda + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.H.L + ERI = self.H.ERI + + #inhomogenous terms + r_Y2 = self.im_Y2.copy() + # Homogenous terms now! + r_Y2 += 0.5 * omega * self.Y2.copy() + r_Y2 += 2.0 * contract('ia,jb->ijab', Y1, hbar.Hov) + r_Y2 -= contract('ja,ib->ijab', Y1, hbar.Hov) + r_Y2 += contract('ijeb,ea->ijab', Y2, hbar.Hvv) + r_Y2 -= contract('im,mjab->ijab', hbar.Hoo, Y2) + r_Y2 += 0.5 * contract('ijmn,mnab->ijab', hbar.Hoooo, Y2) + r_Y2 += 0.5 * contract('ijef,efab->ijab', Y2, hbar.Hvvvv) + r_Y2 += 2.0 * contract('ie,ejab->ijab', Y1, hbar.Hvovv) + r_Y2 -= contract('ie,ejba->ijab', Y1, hbar.Hvovv) + r_Y2 -= 2.0 * contract('mb,jima->ijab', Y1, hbar.Hooov) + r_Y2 += contract('mb,ijma->ijab', Y1, hbar.Hooov) + r_Y2 += 2.0 * contract('ieam,mjeb->ijab', hbar.Hovvo, Y2) + r_Y2 -= contract('iema,mjeb->ijab', hbar.Hovov, Y2) + r_Y2 -= contract('mibe,jema->ijab', Y2, hbar.Hovov) + r_Y2 -= contract('mieb,jeam->ijab', Y2, hbar.Hovvo) + r_Y2 += contract('ijeb,ae->ijab', L[o,o,v,v], cclambda.build_Gvv(t2, Y2)) + r_Y2 -= contract('mi,mjab->ijab', cclambda.build_Goo(t2, Y2), L[o,o,v,v]) + + r_Y2 = r_Y2 + r_Y2.swapaxes(0,1).swapaxes(2,3) + + return r_Y2 + def pseudoresponse(self, pertbar, X1, X2): contract = self.ccwfn.contract polar1 = 2.0 * contract('ai,ia->', np.conj(pertbar.Avo), X1) polar2 = 2.0 * contract('ijab,ijab->', np.conj(pertbar.Avvoo), (2.0*X2 - X2.swapaxes(2,3))) - return -2.0*(polar1 + polar2) - - + return -2.0*(polar1 + polar2) + class pertbar(object): def __init__(self, pert, ccwfn): o = ccwfn.o diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index c90e710..5ab03ae 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -87,13 +87,14 @@ def __init__(self, scf_wfn, **kwargs): self.make_t3_density = kwargs.pop('make_t3_density', False) - valid_local_models = [None, 'PNO', 'PAO','PNO++'] + valid_local_models = [None, 'PNO', 'PAO','CPNO++','PNO++'] local = kwargs.pop('local', None) # TODO: case-protect this kwarg if local not in valid_local_models: raise Exception("%s is not an allowed local-CC model." % (local)) self.local = local self.local_cutoff = kwargs.pop('local_cutoff', 1e-5) + self.ed_omega = kwargs.pop('omega', 0) valid_local_MOs = ['PIPEK_MEZEY', 'BOYS'] local_MOs = kwargs.pop('local_mos', 'PIPEK_MEZEY') @@ -151,7 +152,7 @@ def __init__(self, scf_wfn, **kwargs): self.H = Hamiltonian(self.ref, self.C, self.C, self.C, self.C) if local is not None: - self.Local = Local(local, self.C, self.nfzc, self.no, self.nv, self.H, self.local_cutoff,self.it2_opt) + self.Local = Local(local, self.C, self.nfzc, self.no, self.nv, self.H, self.local_cutoff,self.it2_opt, self.ed_omega) if filter is not True: self.Local.trans_integrals(self.o, self.v) self.Local.overlaps(self.Local.QL) diff --git a/pycc/data/molecules.py b/pycc/data/molecules.py index 0ea22a1..8676250 100644 --- a/pycc/data/molecules.py +++ b/pycc/data/molecules.py @@ -26,6 +26,16 @@ symmetry c1 """ +# CO +co = """ +C 0.000000 0.000000 0.000000 +O 0.000000 0.000000 -2.132000 +noreorient +nocom +units au +symmetry c1 +""" + # Håkon's H2O test case h2o_hek=""" units au @@ -43,8 +53,8 @@ symmetry c1 units bohr """ -### Water molecule +### Water molecule h2o = """ O H 1 1.1 @@ -52,6 +62,15 @@ symmetry c1 """ +### Water molecule from Dalton +h2o_dalton = """ +O 0.000000000000 0.000000000000 0.000000000000 +H 0.000000000000 -0.756689920000 0.585891940000 +H 0.000000000000 0.756689920000 0.585891940000 +units Angstrom +symmetry c1 +""" + ### Water cluster ## Number of water molecules ## 2 @@ -269,10 +288,12 @@ moldict["He"] = he moldict["Be"] = be moldict["LiH"] = lih +moldict["CO"] = co moldict["H2"] = h2 moldict["H2O_HEK"] = h2o_hek moldict["H2O_Teach"] = h2o_tutorial moldict["H2O"] = h2o +moldict["H2O_D"] = h2o_dalton moldict["(H2O)_2"] = h2o_2 moldict["(H2O)_3"] = h2o_3 moldict["(H2O)_4"] = h2o_4 diff --git a/pycc/local.py b/pycc/local.py index 64f735e..73c8905 100644 --- a/pycc/local.py +++ b/pycc/local.py @@ -49,6 +49,7 @@ class Local(object): _build_PAO(): build PAO orbital rotation tensors _build_PNO(): build PNO orbital rotation tensors _build_PNOpp(): build PNO++ orbital rotation tensors + _build_cPNOpp(): build cPNO++ orbital rotation tensors trans_integrals(): transform Fock matrix and ERI from the MO basis to a local basis Notes @@ -58,7 +59,7 @@ class Local(object): to run local MP2, uncomment the necessary lines within the _build_"local" functions which are at the end """ - def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt, + def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt, omega, core_cut=5E-2, lindep_cut=1E-6, e_conv=1e-12, @@ -72,6 +73,7 @@ def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt, self.C = C.to_array() self.local = local self.it2_opt = it2_opt + self.omega = omega self.core_cut = core_cut self.lindep_cut = lindep_cut self.e_conv = e_conv @@ -86,6 +88,8 @@ def _build(self): self._build_PAO() elif self.local.upper() in ["PNO++"]: self._build_PNOpp() + elif self.local.upper() in ["CPNO++"]: + self._build_cPNOpp() else: raise Exception("Not a valid local type!") @@ -326,14 +330,14 @@ def _build_PNO(self): # MP2 loop (optional) if self.it2_opt: self._MP2_loop(t2,self.H.F,self.H.ERI,self.H.L,Dijab) - + print("Computing PNOs. Canonical VMO dim: %d" % (self.nv)) #Constructing the PNO density D = self._pairdensity(t2) - + # Now obtain the Q and L - Q, L, eps, dim = self.QL_tensors(v,self.local,t2,D) + Q, L, eps, dim = self.QL_tensors(v,t2,D,local ='PNO') self.Q = Q # transform between canonical VMO and local spaces self.L = L # transform between local and semicanonical local spaces @@ -375,17 +379,18 @@ def _build_PNOpp(self): Dijab = eps_occ.reshape(-1,1,1,1) + eps_occ.reshape(-1,1,1) - eps_vir.reshape(-1,1) - eps_vir # initial guess amplitudes - t2 = self.H.ERI[o,o,v,v]/Dijab + t2 = self.H.ERI[o,o,v,v]/(Dijab + self.omega) + # MP2 loop (optional) if self.it2_opt: - self._MP2_loop(t2,self.H.F,self.H.ERI,self.H.L,Dijab) + self._MP2_loop(t2,self.H.F,self.H.ERI,self.H.L,Dijab + self.omega) # Construct the perturbed pair density, Eqn. 10 D = self._pert_pairdensity(t2) # Now obtain Q and L - Q, L, eps, dim = self.QL_tensors(v,self.local,t2,D) + Q, L, eps, dim = self.QL_tensors(v,t2,D,local ='PNO++') self.Q = Q # transform between canonical VMO and local spaces self.L = L # transform between local and semicanonical local spaces @@ -401,6 +406,66 @@ def _build_PNOpp(self): self.Q[ji] = self.Q[ij] self.L[ji] = self.L[ij] + def _build_cPNOpp(self): + """ + Perform MP2 loop in non-canonical MO basis, then construct pair density based on t2 amplitudes + then build MP2-level PNO + + Attributes + ---------- + Q: transform between canonical VMO and PNO++ spaces + L: transform between LPNO and semicanonical PNO++ spaces + dim: dimension of cPNO++ space + eps: semicananonical cPNO++ energies + + Notes + ----- + Equations from D'Cunha & Crawford 2021 [10.1021/acs.jctc.0c01086] + """ + v = slice(self.no, self.no+self.nv) + + self._build_PNO() + Q_PNO = self.Q + + self._build_PNOpp() + Q_PNOpp = self.Q + + self.Q = [] # truncated PNO list + self.dim = np.zeros((self.no*self.no), dtype=int) # dimension of local space for each pair + self.L = [] # semicanonical PNO list + self.eps = [] # approximated virtual orbital energies + T2_local = 0 + for ij in range(self.no*self.no): + i = ij // self.no + j = ij % self.no + + Q_comb = np.hstack((Q_PNO[ij], Q_PNOpp[ij])) + Q_ortho, trash = np.linalg.qr(Q_comb) + self.Q.append(Q_ortho) + # Compute semicanonical virtual space + F = Q_ortho.T @ self.H.F[v,v] @ Q_ortho # Fock matrix in local basis + eval, evec = np.linalg.eigh(F) + self.eps.append(eval) + self.L.append(evec) + self.dim[ij] = Q_ortho.shape[1] + T2_local += self.dim[ij] * self.dim[ij] + print(self.local + " dimension of pair %d = %d" % (ij, self.dim[ij])) + + print("Average " + self.local + " dimension: %2.3f" % (np.average(self.dim))) + print("T2 " + self.local + ": %d" % (T2_local)) + T2_full = (self.no*self.no)*(self.nv*self.nv) + print("T2 full: %d" % (T2_full)) + print("T2 Ratio: %3.12f" % (T2_local/T2_full)) + + #temporary way to generate make sure the phase factor of Q_ij and L_ij matches with Q_ji and L_ji + for i in range(self.no): + for j in range(0,i): + ij = i*self.no + j + ji = j*self.no + i + + self.Q[ji] = self.Q[ij] + self.L[ji] = self.L[ij] + def _pert_pairdensity(self,t2): ''' Constructing the approximated perturbed pair density @@ -483,7 +548,7 @@ def _pairdensity(self, t_ijab): D[ij] *= 0.5 return D - def QL_tensors(self,v,local,t2,D): + def QL_tensors(self,v,t2,D,local): # Create list for Q, L and eps Q_full = np.zeros_like(t2.copy().reshape((self.no*self.no,self.nv,self.nv))) Q = [] # truncated PNO list @@ -496,7 +561,7 @@ def QL_tensors(self,v,local,t2,D): for ij in range(self.no*self.no): i = ij // self.no j = ij % self.no - + # Compute local and truncate occ[ij], Q_full[ij] = np.linalg.eigh(D[ij]) if (occ[ij] < 0).any(): # Check for negative occupation numbers @@ -505,7 +570,7 @@ def QL_tensors(self,v,local,t2,D): Using absolute values - please check if your input is correct.".format(neg)) dim[ij] = (np.abs(occ[ij]) > self.cutoff).sum() Q.append(Q_full[ij, :, (self.nv-dim[ij]):]) - + # Compute semicanonical virtual space F = Q[ij].T @ self.H.F[v,v] @ Q[ij] # Fock matrix in local basis eval, evec = np.linalg.eigh(F) @@ -780,6 +845,41 @@ def filter_amps(self, r1, r2): return t1, t2 + def filter_pertamps(self, r1, r2, eps_occ, eps_vir, omega): + no = self.no + nv = self.nv + dim = self.dim + + t1 = np.zeros((no,nv)) + for i in range(no): + ii = i * no + i + + X = self.Q[ii].T @ r1[i] + Y = self.L[ii].T @ X + + for a in range(dim[ii]): + Y[a] = Y[a]/(eps_occ[i] - eps_vir[ii][a] + omega) + + X = self.L[ii] @ Y + t1[i] = self.Q[ii] @ X + + t2 = np.zeros((no,no,nv,nv)) + for ij in range(no*no): + i = ij // no + j = ij % no + + X = self.Q[ij].T @ r2[i,j] @ self.Q[ij] + Y = self.L[ij].T @ X @ self.L[ij] + + for a in range(dim[ij]): + for b in range(dim[ij]): + Y[a,b] = Y[a,b]/(eps_occ[i] + eps_occ[j] - eps_vir[ij][a] - eps_vir[ij][b] + omega) + + X = self.L[ij] @ Y @ self.L[ij].T + t2[i,j] = self.Q[ij] @ X @ self.Q[ij].T + + return t1, t2 + def filter_res(self, r1, r2): no = self.no nv = self.nv diff --git a/pycc/tests/test_016_chk/ref_wfn.npy b/pycc/tests/test_016_chk/ref_wfn.npy deleted file mode 100644 index fee1227..0000000 Binary files a/pycc/tests/test_016_chk/ref_wfn.npy and /dev/null differ diff --git a/pycc/tests/test_036_cpnoppcc.py b/pycc/tests/test_036_cpnoppcc.py new file mode 100644 index 0000000..a265cb7 --- /dev/null +++ b/pycc/tests/test_036_cpnoppcc.py @@ -0,0 +1,80 @@ +""" +Test basic CPNO++-CCSD energy +""" + +# Import package, test suite, and other packages as needed +import psi4 +import pycc +import pytest +from ..data.molecules import * + + +def test_cpnopp_ccsd(): + """H2O CPNO++-CCSD Test""" + # Psi4 Setup + psi4.set_memory('2 GB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'cc-pVDZ', + 'scf_type': 'pk', + 'mp2_type': 'conv', + 'freeze_core': 'false', + 'e_convergence': 1e-13, + 'd_convergence': 1e-13, + 'r_convergence': 1e-13, + 'diis': 8}) + mol = psi4.geometry(moldict["H2O"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + maxiter = 75 + e_conv = 1e-12 + r_conv = 1e-12 + max_diis = 8 + + ccsd = pycc.ccwfn(rhf_wfn, local='CPNO++', local_cutoff=1e-7, it2_opt=False) + eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter) + + hbar = pycc.cchbar(ccsd) + cclambda = pycc.cclambda(ccsd, hbar) + lccsd = cclambda.solve_lambda(e_conv, r_conv, maxiter, max_diis) + + #Ruhee's ccsd_lpno code + esim = -0.22303320613504354 + lsim = -0.21890326836263854 + + assert (abs(esim - eccsd) < 1e-7) + assert (abs(lsim - lccsd) < 1e-7) + +def test_pnopp_ccsd_opt(): + """H2O CPNO++-CCSD with Optimized Initial T2 Amplitudes""" + # Psi4 Setup + psi4.set_memory('2 GB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'cc-pVDZ', + 'scf_type': 'pk', + 'mp2_type': 'conv', + 'freeze_core': 'false', + 'e_convergence': 1e-13, + 'd_convergence': 1e-13, + 'r_convergence': 1e-13, + 'diis': 8}) + mol = psi4.geometry(moldict["H2O"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + maxiter = 75 + e_conv = 1e-12 + r_conv = 1e-12 + max_diis = 8 + + ccsd = pycc.ccwfn(rhf_wfn, local='CPNO++', local_cutoff=1e-7) + eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter) + + hbar = pycc.cchbar(ccsd) + cclambda = pycc.cclambda(ccsd, hbar) + lccsd = cclambda.solve_lambda(e_conv, r_conv, maxiter, max_diis) + + #Comparing against simulation code + esim = -0.223866820104919 + lsim = -0.21966259490352782 + + assert (abs(esim - eccsd) < 1e-7) + assert (abs(lsim - lccsd) < 1e-7) diff --git a/pycc/tests/test_037_quadresp.py b/pycc/tests/test_037_quadresp.py new file mode 100644 index 0000000..be198b2 --- /dev/null +++ b/pycc/tests/test_037_quadresp.py @@ -0,0 +1,93 @@ +""" +Test CCSD quadratic response function. +""" + +# Import package, test suite, and other packages as needed +import psi4 +import pycc +import pytest +from ..data.molecules import * + +def test_ccsd_SHG(): + """ H2O Second Harmonic Generation CCSD Test """ + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + 'freeze_core': 'false', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12}) + mol = psi4.geometry(moldict["H2O_D"]) + 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) + + # SHG frequencies + omega1 = 0.0656 + omega2 = 0.0656 + + resp.pert_quadresp(omega1, omega2) + SHG = resp.hyperpolar() + + #Dalton + H2O_SHG = -19.7591180824 + + # Other system's Dalton SHG value, all validated #PyCC SHG values + #CO_SHG = -32.4314012752 #CO = -32.431401275449 + #HF_SHG = 10.23985062319 #HF = 10.239850624779 + #HCl_SHG = -21.348258444480 #HCl = -21.348258543702 + + assert(abs(SHG - H2O_SHG) < 1e-7) + +def test_ccsd_OR(): + """ CO Optical Refractivity CCSD Test """ + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + 'freeze_core': 'false', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12}) + mol = psi4.geometry(moldict["CO"]) + 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) + + #OR frequencies + omega1 = -0.0656 + omega2 = 0.0656 + + resp.pert_quadresp(omega1, omega2) + OR = resp.hyperpolar() + + #Dalton + CO_OR = -29.3364710172 + + # Other system's Dalton OR value, all validated #PyCC OR values + #H2O_OR = -17.257079066 #H2O = 17.257079081719 + #HF_OR = 9.6113353644 #HF = 9.611335367514 + #HCl_OR = -19.34209144399 #HCl = -19.342091493471 + + assert(abs(OR - CO_OR) < 1e-7) diff --git a/pycc/tests/test_038_sim_quadresp.py b/pycc/tests/test_038_sim_quadresp.py new file mode 100644 index 0000000..a6e2765 --- /dev/null +++ b/pycc/tests/test_038_sim_quadresp.py @@ -0,0 +1,134 @@ +""" +Test simulation code of CCSD quadratic response function. +""" + +# Import package, test suite, and other packages as needed +import psi4 +import pycc +import pytest +from ..data.molecules import * + +def test_PNO_ccsd_SHG(): + """ H2O Second Harmonic Generation PNO-CCSD Test using local filter""" + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + 'freeze_core': 'false', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12}) + mol = psi4.geometry(moldict["H2O_D"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + e_conv = 1e-12 + r_conv = 1e-12 + + cc = pycc.ccwfn(rhf_wfn, local="PNO", local_cutoff=0, 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) + + # SHG frequencies + omega1 = 0.0656 + omega2 = 0.0656 + + resp.pert_quadresp(omega1, omega2) + SHG = resp.hyperpolar() + + #Dalton + H2O_SHG = -19.7591180824 + + ### Untruncated ### + # Other system's Dalton SHG value, all validated #PyCC SHG values + #CO_SHG = -32.4314012752 #CO = -32.431401275449 + #HF_SHG = 10.23985062319 #HF = 10.239850624779 + #HCl_SHG = -21.348258444480 #HCl = -21.348258543702 + + assert(abs(SHG - H2O_SHG) < 1e-7) + +def test_PNOpp_ccsd_OR(): + """ CO Optical Refractivity PNO++-CCSD Test using local filter""" + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + 'freeze_core': 'false', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12}) + mol = psi4.geometry(moldict["CO"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + e_conv = 1e-12 + r_conv = 1e-12 + + cc = pycc.ccwfn(rhf_wfn, local="PNO++", local_cutoff = 1e-7, 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) + + #OR frequencies + omega1 = -0.0656 + omega2 = 0.0656 + + resp.pert_quadresp(omega1, omega2) + OR = resp.hyperpolar() + + #Dalton + CO_OR = -29.3364710172 + + ### Untruncated ### + # Other system's Dalton OR value, all validated #PyCC OR values + #H2O_OR = -17.257079066 #H2O = 17.257079081719 + #HF_OR = 9.6113353644 #HF = 9.611335367514 + #HCl_OR = -19.34209144399 #HCl = -19.342091493471 + + #can only compare to itself + OR_sim = -14.631284641281 + assert(abs(OR - OR_sim) < 1e-10) + +def test_CPNOpp_ccsd_SHG(): + """ H2O Second Harmonic Generation CPNO++-CCSD Test using local filter""" + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + 'freeze_core': 'false', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12}) + mol = psi4.geometry(moldict["H2O_D"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + e_conv = 1e-12 + r_conv = 1e-12 + + #Testing different localization scheme for MOS -> Foster-Boys + cc = pycc.ccwfn(rhf_wfn, local_mos = 'BOYS', local="CPNO++", local_cutoff = 1e-7, 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) + + #SHG frequencies + omega1 = 0.0656 + omega2 = 0.0656 + + resp.pert_quadresp(omega1, omega2) + SHG = resp.hyperpolar() + + #can only compare to itself + SHG_sim = -17.874897845906 + assert(abs(SHG - SHG_sim) < 1e-10) diff --git a/setup.py b/setup.py index a3e9f99..71fd45e 100644 --- a/setup.py +++ b/setup.py @@ -35,9 +35,9 @@ install_requires=[ "numpy", "opt_einsum", - "scipy", - "torch" - ], + "scipy"], + # "torch" + #], # Which Python importable modules should be included when your package is installed # Handled automatically by setuptools. Use 'exclude' to prevent some specific @@ -59,7 +59,7 @@ 'Mac OS-X', 'Unix'], # 'Windows'], # Valid platforms your code works on, adjust to your flavor - python_requires=">=3.11" # Python version restrictions + python_requires=">=3.7" # Python version restrictions # Manual control if final package is compressible or not, set False to prevent the .egg from being made # zip_safe=False,