From 0ef328a9a015f1f244778ab657964e08dbe48aa1 Mon Sep 17 00:00:00 2001 From: jattakumi Date: Fri, 23 Feb 2024 14:02:49 -0500 Subject: [PATCH 01/38] clean draft of asym linear resp --- pycc/ccresponse.py | 417 +++++++++++++++++++++++++++++++++++++- pycc/tests/test_036_lr.py | 78 +++++++ 2 files changed, 491 insertions(+), 4 deletions(-) create mode 100644 pycc/tests/test_036_lr.py diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index 092977c..0ebf0f4 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -1,6 +1,7 @@ """ ccresponse.py: CC Response Functions """ +from opt_einsum import contract if __name__ == "__main__": raise Exception("This file cannot be invoked on its own.") @@ -8,6 +9,7 @@ import numpy as np import time from .utils import helper_diis +from .cclambda import cclambda class ccresponse(object): """ @@ -200,6 +202,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 @@ -283,6 +286,8 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis pertkey = B + "_" + self.cart[axis] X_key = pertkey + "_" + f"{omega:0.6f}" print("Solving right-hand perturbed wave function for %s:" % (X_key)) + #X_2[pertkey] = self.solve_right(self.pertbar[pertkey], omega, e_conv, r_conv, maxiter, max_diis, start_diis) + check.append(polar) 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) if (omega != 0.0): @@ -291,7 +296,58 @@ 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 linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): + + # Defining the l1 and l2 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + + # Please refer to eqn 78 of [Crawford:xxxx]. + # Writing H(1)(omega) = B, T(1)(omega) = X, L(1)(omega) = y + # <> = <0|Y(B) * A_bar|0> + <0| (1 + L(0))[A_bar, X(B)}|0> + # polar1 polar2 + polar1 = 0 + 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|[A_bar, X(B)]|0> + polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B) + # <0|L1(0) [A_bar, X2(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) + polar2 -= contract("ik, ki -> ", tmp, pertbar_A.Aoo) + # <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> + 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, X1(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) + polar2 -= 0.5 * contract("jk, kj -> ", tmp, pertbar_A.Aoo) + tmp = contract("ijab, ijac -> bc", l2, X2_B) + polar2 += 0.5 * contract("bc, bc -> ", tmp, pertbar_A.Avv) + tmp = contract("ijab, ijcb -> ac", l2, X2_B) + polar2 += 0.5 * contract("ac, ac -> ", tmp, pertbar_A.Avv) + + 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): solver_start = time.time() Dia = self.Dia @@ -338,6 +394,63 @@ 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): + 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) + + # need to understand this + 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) + + 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 +504,310 @@ 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 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + cclambda = self.cclambda + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.H.L + + # Inhomogenous terms appearing in Y1 equations + #seems like these imhomogenous terms are computing at the beginning and not involve in the iteration itself + #may require moving to a sperate function + + # 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) + + 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 these next 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) + + #contains regular Gvv as well as Goo, think about just calling it from cclambda instead of generating it + 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() + #homogenous terms appearing in Y1 equations + 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 + + # Inhomogenous terms appearing in Y2 equations + # 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! + # a factor of 0.5 because of the relation/comment just above + # and due to the fact that Y2_ijab = Y2_jiba + 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/tests/test_036_lr.py b/pycc/tests/test_036_lr.py new file mode 100644 index 0000000..3b2bf78 --- /dev/null +++ b/pycc/tests/test_036_lr.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_linresp(): + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + '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_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) From a4318a7e665e17b3c9cdf6533eb483d494b89a17 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Wed, 6 Mar 2024 14:47:27 -0500 Subject: [PATCH 02/38] Working to fix new CI problems. --- .github/workflows/CI.yaml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index dc384ab..b9a4796 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -65,9 +65,10 @@ jobs: - name: Test and generate coverage report run: | - pip install coverage - coverage run -m pytest - coverage xml + # pip install coverage + # coverage run -m pytest + # coverage xml + pytest ls -l - name: Upload coverage reports to Codecov From e605915ab38e78d5465176f228a942d870837cb2 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Wed, 6 Mar 2024 15:19:38 -0500 Subject: [PATCH 03/38] Change pip install. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index b9a4796..52bc132 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -60,7 +60,7 @@ jobs: - name: Install PyCC and Deps run: | - python -m pip install . + pip install -e . pip install pytest - name: Test and generate coverage report From 60ae95f1d7c80d59c996e88454cdb016e1bcbf3c Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Wed, 6 Mar 2024 15:51:30 -0500 Subject: [PATCH 04/38] working --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 52bc132..378b253 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -68,7 +68,7 @@ jobs: # pip install coverage # coverage run -m pytest # coverage xml - pytest + pytest -s ls -l - name: Upload coverage reports to Codecov From 8fa752216fe6608625f2edd8418a56b3036366b0 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Wed, 6 Mar 2024 16:45:35 -0500 Subject: [PATCH 05/38] trying a single test --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 378b253..d5db27c 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -68,7 +68,7 @@ jobs: # pip install coverage # coverage run -m pytest # coverage xml - pytest -s + pytest -s pycc/tests/test_002_ccsd_energy.py ls -l - name: Upload coverage reports to Codecov From fb8316b6847eb17e5d0c47e42336b03d0dd43237 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Wed, 6 Mar 2024 16:53:50 -0500 Subject: [PATCH 06/38] testing libomp --- .github/workflows/CI.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index d5db27c..ffe529a 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -68,6 +68,7 @@ jobs: # pip install coverage # coverage run -m pytest # coverage xml + export KMP_DUPLICATE_LIB_OK=TRUE pytest -s pycc/tests/test_002_ccsd_energy.py ls -l From b3263fdb39f9379141b2bd85e973a0e361738eda Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Wed, 6 Mar 2024 16:59:02 -0500 Subject: [PATCH 07/38] Trying all tests. --- .github/workflows/CI.yaml | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index ffe529a..59c89a1 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -65,11 +65,10 @@ jobs: - name: Test and generate coverage report run: | - # pip install coverage - # coverage run -m pytest - # coverage xml export KMP_DUPLICATE_LIB_OK=TRUE - pytest -s pycc/tests/test_002_ccsd_energy.py + pip install coverage + coverage run -m pytest + coverage xml ls -l - name: Upload coverage reports to Codecov From f06543ee2c4c9493bc7f828e76befcaee1a9e1ed Mon Sep 17 00:00:00 2001 From: Jose Marc Madriaga Date: Mon, 8 Apr 2024 10:35:04 -0400 Subject: [PATCH 08/38] Added CPNO++ to local --- pycc/ccwfn.py | 2 +- pycc/local.py | 77 ++++++++++++++++++++++++++++--- pycc/tests/test_035_cpnoppcc.py | 80 +++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 8 deletions(-) create mode 100644 pycc/tests/test_035_cpnoppcc.py diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index bce9abd..e308ae9 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -87,7 +87,7 @@ 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: diff --git a/pycc/local.py b/pycc/local.py index 64f735e..9e01847 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 @@ -86,6 +87,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 +329,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 @@ -385,7 +388,7 @@ def _build_PNOpp(self): 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 +404,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 +546,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 +559,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 +568,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) diff --git a/pycc/tests/test_035_cpnoppcc.py b/pycc/tests/test_035_cpnoppcc.py new file mode 100644 index 0000000..a265cb7 --- /dev/null +++ b/pycc/tests/test_035_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) From 4ad6b4a29024a1f16b6bc380604a1d8b0ab6057a Mon Sep 17 00:00:00 2001 From: jattakumi Date: Fri, 12 Apr 2024 17:42:25 -0400 Subject: [PATCH 09/38] Pull request comments resolved --- pycc/ccresponse.py | 33 ++++++++++++++++++++++++++------- 1 file changed, 26 insertions(+), 7 deletions(-) diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index 0ebf0f4..a386d8a 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -1,7 +1,6 @@ """ ccresponse.py: CC Response Functions """ -from opt_einsum import contract if __name__ == "__main__": raise Exception("This file cannot be invoked on its own.") @@ -286,7 +285,7 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis pertkey = B + "_" + self.cart[axis] X_key = pertkey + "_" + f"{omega:0.6f}" print("Solving right-hand perturbed wave function for %s:" % (X_key)) - #X_2[pertkey] = self.solve_right(self.pertbar[pertkey], omega, e_conv, r_conv, maxiter, max_diis, start_diis) + X_2[pertkey] = self.solve_right(self.pertbar[pertkey], omega, e_conv, r_conv, maxiter, max_diis, start_diis) check.append(polar) 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) @@ -298,12 +297,28 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis 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). + + The linear response function, <> 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. + """ # Defining the l1 and l2 l1 = self.cclambda.l1 l2 = self.cclambda.l2 - # Please refer to eqn 78 of [Crawford:xxxx]. + # Please refer to eqn 78 of [Crawford: https://crawford.chem.vt.edu/wp-content/uploads/2022/06/cc_response.pdf]. # Writing H(1)(omega) = B, T(1)(omega) = X, L(1)(omega) = y # <> = <0|Y(B) * A_bar|0> + <0| (1 + L(0))[A_bar, X(B)}|0> # polar1 polar2 @@ -395,6 +410,14 @@ def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, m 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 inhomogeneous 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 homogeneous terms. + """ + solver_start = time.time() Dia = self.Dia @@ -514,10 +537,6 @@ def in_Y1(self, pertbar, X1, X2): t2 = self.ccwfn.t2 hbar = self.hbar L = self.H.L - - # Inhomogenous terms appearing in Y1 equations - #seems like these imhomogenous terms are computing at the beginning and not involve in the iteration itself - #may require moving to a sperate function # good r_Y1 = 2.0 * pertbar.Aov.copy() From 9ebf057193c2c203d696c9b714b83de4dd8c2c28 Mon Sep 17 00:00:00 2001 From: jattakumi Date: Wed, 17 Apr 2024 15:08:39 -0400 Subject: [PATCH 10/38] Fixed typos --- pycc/ccresponse.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index a386d8a..9f51250 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -297,7 +297,7 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis 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). The linear response function, <> generally reuires the following perturbed wave functions and frequencies: @@ -312,7 +312,9 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): ------ 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 @@ -410,13 +412,13 @@ def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, m 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 inhomogeneous 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 homogeneous terms. - """ + """ solver_start = time.time() From b5aa312c852930fe206da12f1a077566ae958c8e Mon Sep 17 00:00:00 2001 From: Ethan Fink Date: Mon, 29 Apr 2024 20:21:45 -0400 Subject: [PATCH 11/38] Mini conda version for CI --- .github/workflows/CI.yaml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 59c89a1..75d8afe 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -35,8 +35,9 @@ jobs: ulimit -a - name: Create Psi4 Environment - uses: conda-incubator/setup-miniconda@v2 + uses: conda-incubator/setup-miniconda@v3 with: + miniconda-version: "latest" activate-environment: p4env environment-file: devtools/conda-envs/psi.yaml python-version: ${{ matrix.python-version }} From 754fddff35946c9a53465a9c832cfc8055e0baca Mon Sep 17 00:00:00 2001 From: Ethan Fink Date: Mon, 29 Apr 2024 20:24:16 -0400 Subject: [PATCH 12/38] actions/checkout version update --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 75d8afe..e23434d 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -25,7 +25,7 @@ jobs: python-version: ['3.11'] steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - name: Additional info about the build shell: bash From 6e359fb3f55e71348c18ea7c3358a535932d0750 Mon Sep 17 00:00:00 2001 From: Zhe Wang Date: Wed, 1 May 2024 13:33:38 -0400 Subject: [PATCH 13/38] Missing terms for RT-CC3 with external perturbation --- pycc/ccdensity.py | 45 ++++++++++++++++++++++++++++--------- pycc/cclambda.py | 16 +++++++++++-- pycc/cctriples.py | 57 +++++++++++++++++++++++++++++++++++++++++++++++ pycc/ccwfn.py | 17 +++++++++----- pycc/rt/rtcc.py | 6 ++--- 5 files changed, 121 insertions(+), 20 deletions(-) diff --git a/pycc/ccdensity.py b/pycc/ccdensity.py index f911474..c676a5f 100644 --- a/pycc/ccdensity.py +++ b/pycc/ccdensity.py @@ -8,7 +8,7 @@ import time import numpy as np import torch -from .cctriples import t3c_ijk, t3c_abc, l3_ijk, l3_abc, t3c_bc, l3_bc +from .cctriples import t3c_ijk, t3c_abc, l3_ijk, l3_abc, t3c_bc, l3_bc, t3_pert_ijk, t3_pert_bc class ccdensity(object): """ @@ -152,7 +152,7 @@ def compute_energy(self): return ecc - def compute_onepdm(self, t1, t2, l1, l2): + def compute_onepdm(self, t1, t2, l1, l2, real_time=False): """ Parameters ---------- @@ -196,7 +196,7 @@ def compute_onepdm(self, t1, t2, l1, l2): Wvovv = self.ccwfn.build_cc3_Wamef(o, v, ERI, t1) Wooov = self.ccwfn.build_cc3_Wmnie(o, v, ERI, t1) - opdm[o,v] += self.build_cc3_Dov(o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov) + opdm[o,v] += self.build_cc3_Dov(o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov, real_time=real_time) # Density matrix blocks in contractions with T1-transformed dipole integrals if isinstance(t1, torch.Tensor): @@ -274,7 +274,7 @@ def build_Dov(self, t1, t2, l1, l2): # complete return Dov # CC3 contributions to the one electron densities - def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov): + def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, Wvovv, Wooov, real_time=False): contract = self.contract if isinstance(t1, torch.Tensor): Dov = torch.zeros_like(t1) @@ -290,13 +290,19 @@ def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, W Zlmdi[i,j] += contract('def,ife->di', l3, t2[k]) # Dov_1 t3 = t3c_ijk(o, v, i, j, k, t2, Wvvvo, Wovoo, F, contract) + if real_time is True: + if isinstance(t1, torch.Tensor): + V = F - self.ccwfn.H.F.clone() + else: + V = F - self.ccwfn.H.F.copy() + t3 -= t3_pert_ijk(o, v, i, j, k, t2, V, F, contract) Dov[i] += contract('abc,bc->a', t3 - t3.swapaxes(0,1), l2[j,k]) # Dov_2 Dov -= contract('lmdi, lmda->ia', Zlmdi, t2) return Dov - def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov): + def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False): contract = self.contract if isinstance(l1, torch.Tensor): Doo = torch.zeros_like(l1[:,:no]) @@ -305,12 +311,18 @@ def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv for b in range(nv): for c in range(nv): t3 = t3c_bc(o, v, b, c, t2, Wvvvo, Wovoo, F, contract) + if real_time is True: + if isinstance(t2, torch.Tensor): + V = F - self.ccwfn.H.F.clone() + else: + V = F - self.ccwfn.H.F.copy() + t3 -= t3_pert_bc(o, v, b, c, t2, V, F, contract) l3 = l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract) Doo -= 0.5 * contract('lmia,lmja->ij', t3, l3) return Doo - def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov): + def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv, Wooov, real_time=False): contract = self.contract if isinstance(l1, torch.Tensor): Dvv = torch.zeros_like(l1) @@ -322,6 +334,12 @@ def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv for j in range(no): for k in range(no): t3 = t3c_ijk(o, v, i, j, k, t2, Wvvvo, Wovoo, F, contract) + if real_time is True: + if isinstance(t2, torch.Tensor): + V = F - self.ccwfn.H.F.clone() + else: + V = F - self.ccwfn.H.F.copy() + t3 -= t3_pert_ijk(o, v, i, j, k, t2, V, F, contract) l3 = l3_ijk(i, j, k, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract) Dvv += 0.5 * contract('bdc,adc->ab', t3, l3) @@ -380,6 +398,9 @@ def build_Dooov(self, t1, t2, l1, l2): # complete Dooov += contract('jake,ie->ijka', tmp, t1) tmp = contract('imea,kmef->iakf', t2, l2) Dooov += contract('iakf,jf->ijka', tmp, t1) + + if isinstance(tmp, torch.Tensor): + del tmp, Goo tmp = contract('kmef,jf->kmej', l2, t1) tmp = contract('kmej,ie->kmij', tmp, t1) @@ -391,7 +412,6 @@ def build_Dooov(self, t1, t2, l1, l2): # complete if isinstance(tmp, torch.Tensor): del tmp - del Goo return Dooov @@ -431,6 +451,9 @@ def build_Dvvvo(self, t1, t2, l1, l2): # complete Dvvvo -= contract('iamc,mb->abci', tmp, t1) tmp = contract('mibe,nmce->ibnc', t2, l2) Dvvvo -= contract('ibnc,na->abci', tmp, t1) + + if isinstance(tmp, torch.Tensor): + del tmp, Gvv tmp = contract('nmce,ie->nmci', l2, t1) tmp = contract('nmci,na->amci', tmp, t1) @@ -441,8 +464,7 @@ def build_Dvvvo(self, t1, t2, l1, l2): # complete Dvvvo += self.ccwfn.Gvvvo if isinstance(tmp, torch.Tensor): - del tmp - del Gvv + del tmp return Dvvvo @@ -550,6 +572,9 @@ def build_Doovv(self, t1, t2, l1, l2): tmp = contract('if,mnef->mnei', t1, l2) tmp = contract('mnei,njae->mija', tmp, t2) Doovv += contract('mb,mija->ijab', t1, tmp) + + if isinstance(tmp, torch.Tensor): + del tmp, tmp1, tmp2, Goo, Gvv tmp = contract('jf,mnef->mnej', t1, l2) tmp = contract('ie,mnej->mnij', t1, tmp) @@ -561,7 +586,7 @@ def build_Doovv(self, t1, t2, l1, l2): Doovv += self.ccwfn.Goovv if isinstance(tmp, torch.Tensor): - del tmp, tmp1, tmp2, Goo, Gvv + del tmp return Doovv diff --git a/pycc/cclambda.py b/pycc/cclambda.py index aad49ec..85c62d8 100644 --- a/pycc/cclambda.py +++ b/pycc/cclambda.py @@ -11,7 +11,7 @@ from opt_einsum import contract from .utils import helper_diis import torch -from .cctriples import t3c_ijk, l3_ijk, l3_ijk_alt +from .cctriples import t3c_ijk, l3_ijk, l3_ijk_alt, t3_pert_ijk class cclambda(object): @@ -346,6 +346,12 @@ def residuals(self, F, t1, t2, l1, l2): for n in range(no): for l in range(no): t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True) + if self.ccwfn.real_time is True: + if isinstance(t1, torch.Tensor): + V = F - self.ccwfn.H.F.clone() + else: + V = F - self.ccwfn.H.F.copy() + t3_lmn -= t3_pert_ijk(o, v, l, m, n, t2, V, F, contract) Zmndi[m,n] += contract('def,ief->di', t3_lmn, ERI[o,l,v,v]) Zmndi[m,n] -= contract('fed,ief->di', t3_lmn, L[o,l,v,v]) Zmdfa[m] += contract('def,ea->dfa', t3_lmn, ERI[n,l,v,v]) @@ -386,7 +392,13 @@ def residuals(self, F, t1, t2, l1, l2): for l in range(no): for m in range(no): for n in range(no): - t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True) + t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True) + if self.ccwfn.real_time is True: + if isinstance(t1, torch.Tensor): + V = F - self.ccwfn.H.F.clone() + else: + V = F - self.ccwfn.H.F.copy() + t3_lmn -= t3_pert_ijk(o, v, l, m, n, t2, V, F, contract) Znf[n] += contract('de,def->f', l2[l,m], (t3_lmn - t3_lmn.swapaxes(0,2))) for m in range(no): Y1 += contract('idf,dfa->ia', l2[:,m], Zmdfa[m]) diff --git a/pycc/cctriples.py b/pycc/cctriples.py index 82f9526..8a4e27e 100644 --- a/pycc/cctriples.py +++ b/pycc/cctriples.py @@ -542,3 +542,60 @@ def l3_bc(b, c, o, v, L, l1, l2, Fov, Wvovv, Wooov, F, contract, WithDenom=True) else: return l3 +# Useful for RT-CC3 +# Additional term in T3 equation when an external perturbation is present +def t3_pert_ijk(o, v, i, j, k, t2, V, F, contract, WithDenom=True): + tmp = contract('ld,ad->al', V[o,v], t2[i,j]) + t3 = contract('al,lcb->abc', tmp, t2[k]) + + if WithDenom is True: + if isinstance(t2, torch.Tensor): + Fv = torch.diag(F)[v] + denom = torch.zeros_like(t3) + else: + Fv = np.diag(F)[v] + denom = np.zeros_like(t3) + denom -= Fv.reshape(-1,1,1) + Fv.reshape(-1,1) + Fv + denom += F[i,i] + F[j,j] + F[k,k] + return t3/denom + else: + return t3 + +def t3_pert_abc(o, v, a, b, c, t2, V, F, contract, WithDenom=True): + tmp = contract('ld,ijd->ijl', V[o,v], t2[:,:,a]) + t3 = contract('ijl,kl->ijk', tmp, t2[:,:,c,b]) + + if WithDenom is True: + no = o.stop + if isinstance(t2, torch.Tensor): + Fo = torch.diag(F)[o] + denom = torch.zeros_like(t3) + else: + Fo = np.diag(F)[o] + denom = np.zeros_like(t3) + denom += Fo.reshape(-1,1,1) + Fo.reshape(-1,1) + Fo + denom -= F[a+no,a+no] + F[b+no,b+no] + F[c+no,c+no] + return t3/denom + else: + return t3 + +def t3_pert_bc(o, v, b, c, t2, V, F, contract, WithDenom=True): + tmp = contract('ld,ijad->ijal', V[o,v], t2) + t3 = contract('ijal,kl->ijka', tmp, t2[:,:,c,b]) + + if WithDenom is True: + no = o.stop + if isinstance(t2, torch.Tensor): + Fo = torch.diag(F)[o] + Fv = torch.diag(F)[v] + denom = torch.zeros_like(t3) + else: + Fo = np.diag(F)[o] + Fv = np.diag(F)[v] + denom = np.zeros_like(t3) + denom += Fo.reshape(-1,1,1,1) + Fo.reshape(-1,1,1) + Fo.reshape(-1,1) + denom -= Fv.reshape(1,1,1,-1) + denom -= F[b+no,b+no] + F[c+no,c+no] + return t3/denom + else: + return t3 diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index f0055d0..7ba28cf 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -13,7 +13,7 @@ from .utils import helper_diis, cc_contract from .hamiltonian import Hamiltonian from .local import Local -from .cctriples import t_tjl, t3c_ijk, t3d_ijk, t3c_abc, t3d_abc +from .cctriples import t_tjl, t3c_ijk, t3d_ijk, t3c_abc, t3d_abc, t3_pert_ijk from .lccwfn import lccwfn class ccwfn(object): @@ -87,6 +87,9 @@ def __init__(self, scf_wfn, **kwargs): self.make_t3_density = kwargs.pop('make_t3_density', False) + # RT-CC3 calculations requiring additional terms when an external perturbation is present + self.real_time = kwargs.pop('real_time', False) + valid_local_models = [None, 'PNO', 'PAO','CPNO++','PNO++'] local = kwargs.pop('local', None) # TODO: case-protect this kwarg @@ -325,7 +328,7 @@ def solve_cc(self, e_conv=1e-7, r_conv=1e-7, maxiter=100, max_diis=8, start_diis if niter >= start_diis: self.t1, self.t2 = diis.extrapolate(self.t1, self.t2) - def residuals(self, F, t1, t2): + def residuals(self, F, t1, t2, real_time=False): """ Parameters ---------- @@ -382,9 +385,13 @@ def residuals(self, F, t1, t2): for i in range(no): for j in range(no): for k in range(no): - t3 = t3c_ijk(o, v, i, j, k, t2, Wabei_cc3, Wmbij_cc3, F, -contract, WithDenom=True) - + t3 = t3c_ijk(o, v, i, j, k, t2, Wabei_cc3, Wmbij_cc3, F, contract, WithDenom=True) + if real_time is True: + if isinstance(t1, torch.Tensor): + V = F - self.H.F.clone() + else: + V = F - self.H.F.copy() + t3 -= t3_pert_ijk(o, v, i, j, k, t2, V, F, contract) X1[i] += contract('abc,bc->a', t3 - t3.swapaxes(0,2), L[j,k,v,v]) X2[i,j] += contract('abc,c->ab', t3 - t3.swapaxes(0,2), Fme[k]) X2[i,j] += contract('abc,dbc->ad', 2 * t3 - t3.swapaxes(1,2) - t3.swapaxes(0,2), Wamef_cc3.swapaxes(0,1)[k]) diff --git a/pycc/rt/rtcc.py b/pycc/rt/rtcc.py index 91fae0d..59d1e17 100644 --- a/pycc/rt/rtcc.py +++ b/pycc/rt/rtcc.py @@ -122,7 +122,7 @@ def f(self, t, y): F = self.ccwfn.H.F.copy() + self.mu_tot * self.V(t) # Compute the current residuals - rt1, rt2 = self.ccwfn.residuals(F, t1, t2) + rt1, rt2 = self.ccwfn.residuals(F, t1, t2, real_time=True) rt1 = rt1 * (-1.0j) rt2 = rt2 * (-1.0j) if self.ccwfn.local is not None: @@ -207,7 +207,7 @@ def extract_amps(self, y): return t1, t2, l1, l2, phase - def dipole(self, t1, t2, l1, l2, magnetic = False): + def dipole(self, t1, t2, l1, l2, magnetic = False, real_time=False): """ Parameters ---------- @@ -222,7 +222,7 @@ def dipole(self, t1, t2, l1, l2, magnetic = False): Cartesian components of the correlated dipole moment """ if self.ccwfn.model == 'CC3': - (opdm, opdm_cc3) = self.ccdensity.compute_onepdm(t1, t2, l1, l2) + (opdm, opdm_cc3) = self.ccdensity.compute_onepdm(t1, t2, l1, l2, real_time=real_time) else: opdm = self.ccdensity.compute_onepdm(t1, t2, l1, l2) From 3615e2268381a5cc122b24acd14e654fbcb273c5 Mon Sep 17 00:00:00 2001 From: Zhe Wang Date: Wed, 1 May 2024 18:39:18 -0400 Subject: [PATCH 14/38] Fixed errors after testing --- pycc/ccdensity.py | 12 +++---- pycc/cclambda.py | 4 +-- pycc/ccwfn.py | 2 +- pycc/rt/lasers.py | 31 ++++++++++++++++ pycc/tests/test_037_rtcc3.py | 70 ++++++++++++++++++++++++++++++++++++ 5 files changed, 110 insertions(+), 9 deletions(-) create mode 100644 pycc/tests/test_037_rtcc3.py diff --git a/pycc/ccdensity.py b/pycc/ccdensity.py index c676a5f..cdc42ab 100644 --- a/pycc/ccdensity.py +++ b/pycc/ccdensity.py @@ -290,7 +290,7 @@ def build_cc3_Dov(self, o, v, no, nv, F, L, t1, t2, l1, l2, Wvvvo, Wovoo, Fov, W Zlmdi[i,j] += contract('def,ife->di', l3, t2[k]) # Dov_1 t3 = t3c_ijk(o, v, i, j, k, t2, Wvvvo, Wovoo, F, contract) - if real_time is True: + if real_time is True: if isinstance(t1, torch.Tensor): V = F - self.ccwfn.H.F.clone() else: @@ -311,7 +311,7 @@ def build_cc3_Doo(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv for b in range(nv): for c in range(nv): t3 = t3c_bc(o, v, b, c, t2, Wvvvo, Wovoo, F, contract) - if real_time is True: + if real_time is True: if isinstance(t2, torch.Tensor): V = F - self.ccwfn.H.F.clone() else: @@ -334,7 +334,7 @@ def build_cc3_Dvv(self, o, v, no, nv, F, L, t2, l1, l2, Fov, Wvvvo, Wovoo, Wvovv for j in range(no): for k in range(no): t3 = t3c_ijk(o, v, i, j, k, t2, Wvvvo, Wovoo, F, contract) - if real_time is True: + if real_time is True: if isinstance(t2, torch.Tensor): V = F - self.ccwfn.H.F.clone() else: @@ -399,7 +399,7 @@ def build_Dooov(self, t1, t2, l1, l2): # complete tmp = contract('imea,kmef->iakf', t2, l2) Dooov += contract('iakf,jf->ijka', tmp, t1) - if isinstance(tmp, torch.Tensor): + if isinstance(tmp, torch.Tensor): del tmp, Goo tmp = contract('kmef,jf->kmej', l2, t1) @@ -452,7 +452,7 @@ def build_Dvvvo(self, t1, t2, l1, l2): # complete tmp = contract('mibe,nmce->ibnc', t2, l2) Dvvvo -= contract('ibnc,na->abci', tmp, t1) - if isinstance(tmp, torch.Tensor): + if isinstance(tmp, torch.Tensor): del tmp, Gvv tmp = contract('nmce,ie->nmci', l2, t1) @@ -573,7 +573,7 @@ def build_Doovv(self, t1, t2, l1, l2): tmp = contract('mnei,njae->mija', tmp, t2) Doovv += contract('mb,mija->ijab', t1, tmp) - if isinstance(tmp, torch.Tensor): + if isinstance(tmp, torch.Tensor): del tmp, tmp1, tmp2, Goo, Gvv tmp = contract('jf,mnef->mnej', t1, l2) diff --git a/pycc/cclambda.py b/pycc/cclambda.py index 85c62d8..a958971 100644 --- a/pycc/cclambda.py +++ b/pycc/cclambda.py @@ -346,7 +346,7 @@ def residuals(self, F, t1, t2, l1, l2): for n in range(no): for l in range(no): t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True) - if self.ccwfn.real_time is True: + if self.ccwfn.real_time is True: if isinstance(t1, torch.Tensor): V = F - self.ccwfn.H.F.clone() else: @@ -393,7 +393,7 @@ def residuals(self, F, t1, t2, l1, l2): for m in range(no): for n in range(no): t3_lmn = t3c_ijk(o, v, l, m, n, t2, Wvvvo, Wovoo, F, contract, WithDenom=True) - if self.ccwfn.real_time is True: + if self.ccwfn.real_time is True: if isinstance(t1, torch.Tensor): V = F - self.ccwfn.H.F.clone() else: diff --git a/pycc/ccwfn.py b/pycc/ccwfn.py index 7ba28cf..d6fdc6d 100644 --- a/pycc/ccwfn.py +++ b/pycc/ccwfn.py @@ -88,7 +88,7 @@ def __init__(self, scf_wfn, **kwargs): self.make_t3_density = kwargs.pop('make_t3_density', False) # RT-CC3 calculations requiring additional terms when an external perturbation is present - self.real_time = kwargs.pop('real_time', False) + self.real_time = kwargs.pop('real_time', False) valid_local_models = [None, 'PNO', 'PAO','CPNO++','PNO++'] local = kwargs.pop('local', None) diff --git a/pycc/rt/lasers.py b/pycc/rt/lasers.py index ea1ee7b..82549f2 100644 --- a/pycc/rt/lasers.py +++ b/pycc/rt/lasers.py @@ -58,3 +58,34 @@ def __call__(self, t): pulse = 0 return pulse +# ramped continuous wave (RCW) +# set nr=0 for a regular cosine wave +class lrcw_laser: + def __init__(self, F_str, omega, nr): + self.F_str = F_str + self.omega = omega + self.nr = nr + def __call__(self, t): + tc = 2 * np.pi / self.omega * self.nr + if t <= tc: + pulse = t / tc * self.F_str * np.cos(self.omega * t) + else: + pulse = self.F_str * np.cos(self.omega * t) + return pulse + +class qrcw_laser: + def __init__(self, F_str, omega, nr): + self.F_str = F_str + self.omega = omega + self.nr = nr + def __call__(self, t): + tc = 2 * np.pi / self.omega * self.nr + if t <= 0.5 * tc: + pulse = 2 * t ** 2 / tc ** 2 * self.F_str * np.cos(self.omega * t) + elif t <= tc: + pulse = (1 - 2 * (t - tc) ** 2 / tc ** 2)* self.F_str * np.cos(self.omega * t) + else: + pulse = self.F_str * np.cos(self.omega * t) + return pulse + + diff --git a/pycc/tests/test_037_rtcc3.py b/pycc/tests/test_037_rtcc3.py new file mode 100644 index 0000000..679d10a --- /dev/null +++ b/pycc/tests/test_037_rtcc3.py @@ -0,0 +1,70 @@ +""" +Test RT-CC3 propagation on water molecule. +""" + +# Import package, test suite, and other packages as needed +import psi4 +import pycc +import pytest +from pycc.rt.integrators import rk4 +from pycc.rt.lasers import qrcw_laser +from ..data.molecules import * + +def test_rtcc_he_cc_pvdz(): + """H2O cc-pVDZ""" + 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': 'false', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12, + 'diis': 1}) + mol = psi4.geometry(moldict["H2O_Teach"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + e_conv = 1e-12 + r_conv = 1e-12 + + cc = pycc.ccwfn(rhf_wfn, model='CC3', real_time=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) + + ccdensity = pycc.ccdensity(cc, cclambda) + + # quadratic ramped continuous wave (QRCW) + F_str = 0.002 + omega = 0.078 + nr = 1 + V = qrcw_laser(F_str, omega, nr) + + # RT-CC Setup + phase = 0 + t0 = 0 + tf = 0.05 + h = 0.01 + t = t0 + rtcc = pycc.rtcc(cc, cclambda, ccdensity, V, kick='x') + y0 = rtcc.collect_amps(cc.t1, cc.t2, cclambda.l1, cclambda.l2, phase) + y = y0 + ODE = rk4(h) + t1, t2, l1, l2, phase = rtcc.extract_amps(y0) + mu0_x, mu0_y, mu0_z = rtcc.dipole(t1, t2, l1, l2) + + mu_z_ref = -0.0859645691 #CFOUR + + while t < tf: + y = ODE(rtcc.f, t, y) + t += h + t1, t2, l1, l2, phase = rtcc.extract_amps(y) + mu_x, mu_y, mu_z = rtcc.dipole(t1, t2, l1, l2, real_time=True) + + print(mu_z) + assert (abs(mu_z_ref - mu_z.real) < 1e-10) + From 67a7b5cc8a8e9eb7c0543aa9df832d0299fe386b Mon Sep 17 00:00:00 2001 From: Zhe Wang Date: Wed, 1 May 2024 19:45:01 -0400 Subject: [PATCH 15/38] Resolved warnings for RT-CC3 test --- pycc/rt/rtcc.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pycc/rt/rtcc.py b/pycc/rt/rtcc.py index 59d1e17..e75ee02 100644 --- a/pycc/rt/rtcc.py +++ b/pycc/rt/rtcc.py @@ -240,6 +240,10 @@ def dipole(self, t1, t2, l1, l2, magnetic = False, real_time=False): else: ints_cc3 = np.zeros_like(ints) for i in range(3): + if isinstance(t1, torch.Tensor): + ints_cc3 = ints_cc3.type_as(t1) + else: + ints_cc3 = ints_cc3.astype(t1.dtype) ints_cc3[i][:no,:no] = self.ccdensity.build_Moo(no, nv, ints[i], t1) ints_cc3[i][-nv:,-nv:] = self.ccdensity.build_Mvv(no, nv, ints[i], t1) From d306038ef9eb3f3a22cff2b8974fd1dc8a8e366f Mon Sep 17 00:00:00 2001 From: jattakumi Date: Fri, 23 Feb 2024 14:02:49 -0500 Subject: [PATCH 16/38] clean draft of asym linear resp --- pycc/ccresponse.py | 417 +++++++++++++++++++++++++++++++++++++- pycc/tests/test_036_lr.py | 78 +++++++ 2 files changed, 491 insertions(+), 4 deletions(-) create mode 100644 pycc/tests/test_036_lr.py diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index 092977c..0ebf0f4 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -1,6 +1,7 @@ """ ccresponse.py: CC Response Functions """ +from opt_einsum import contract if __name__ == "__main__": raise Exception("This file cannot be invoked on its own.") @@ -8,6 +9,7 @@ import numpy as np import time from .utils import helper_diis +from .cclambda import cclambda class ccresponse(object): """ @@ -200,6 +202,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 @@ -283,6 +286,8 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis pertkey = B + "_" + self.cart[axis] X_key = pertkey + "_" + f"{omega:0.6f}" print("Solving right-hand perturbed wave function for %s:" % (X_key)) + #X_2[pertkey] = self.solve_right(self.pertbar[pertkey], omega, e_conv, r_conv, maxiter, max_diis, start_diis) + check.append(polar) 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) if (omega != 0.0): @@ -291,7 +296,58 @@ 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 linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): + + # Defining the l1 and l2 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + + # Please refer to eqn 78 of [Crawford:xxxx]. + # Writing H(1)(omega) = B, T(1)(omega) = X, L(1)(omega) = y + # <> = <0|Y(B) * A_bar|0> + <0| (1 + L(0))[A_bar, X(B)}|0> + # polar1 polar2 + polar1 = 0 + 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|[A_bar, X(B)]|0> + polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B) + # <0|L1(0) [A_bar, X2(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) + polar2 -= contract("ik, ki -> ", tmp, pertbar_A.Aoo) + # <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> + 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, X1(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) + polar2 -= 0.5 * contract("jk, kj -> ", tmp, pertbar_A.Aoo) + tmp = contract("ijab, ijac -> bc", l2, X2_B) + polar2 += 0.5 * contract("bc, bc -> ", tmp, pertbar_A.Avv) + tmp = contract("ijab, ijcb -> ac", l2, X2_B) + polar2 += 0.5 * contract("ac, ac -> ", tmp, pertbar_A.Avv) + + 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): solver_start = time.time() Dia = self.Dia @@ -338,6 +394,63 @@ 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): + 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) + + # need to understand this + 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) + + 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 +504,310 @@ 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 + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + cclambda = self.cclambda + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.H.L + + # Inhomogenous terms appearing in Y1 equations + #seems like these imhomogenous terms are computing at the beginning and not involve in the iteration itself + #may require moving to a sperate function + + # 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) + + 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 these next 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) + + #contains regular Gvv as well as Goo, think about just calling it from cclambda instead of generating it + 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() + #homogenous terms appearing in Y1 equations + 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 + + # Inhomogenous terms appearing in Y2 equations + # 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! + # a factor of 0.5 because of the relation/comment just above + # and due to the fact that Y2_ijab = Y2_jiba + 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/tests/test_036_lr.py b/pycc/tests/test_036_lr.py new file mode 100644 index 0000000..3b2bf78 --- /dev/null +++ b/pycc/tests/test_036_lr.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_linresp(): + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pvdz', + 'scf_type': 'pk', + '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_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) From 95a3a020457a6d691b21783ffa669fa88e6ede7f Mon Sep 17 00:00:00 2001 From: jattakumi Date: Fri, 12 Apr 2024 17:42:25 -0400 Subject: [PATCH 17/38] Pull request comments resolved --- pycc/ccresponse.py | 33 ++++++++++++++++++++++++++------- 1 file changed, 26 insertions(+), 7 deletions(-) diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index 0ebf0f4..a386d8a 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -1,7 +1,6 @@ """ ccresponse.py: CC Response Functions """ -from opt_einsum import contract if __name__ == "__main__": raise Exception("This file cannot be invoked on its own.") @@ -286,7 +285,7 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis pertkey = B + "_" + self.cart[axis] X_key = pertkey + "_" + f"{omega:0.6f}" print("Solving right-hand perturbed wave function for %s:" % (X_key)) - #X_2[pertkey] = self.solve_right(self.pertbar[pertkey], omega, e_conv, r_conv, maxiter, max_diis, start_diis) + X_2[pertkey] = self.solve_right(self.pertbar[pertkey], omega, e_conv, r_conv, maxiter, max_diis, start_diis) check.append(polar) 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) @@ -298,12 +297,28 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis 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). + + The linear response function, <> 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. + """ # Defining the l1 and l2 l1 = self.cclambda.l1 l2 = self.cclambda.l2 - # Please refer to eqn 78 of [Crawford:xxxx]. + # Please refer to eqn 78 of [Crawford: https://crawford.chem.vt.edu/wp-content/uploads/2022/06/cc_response.pdf]. # Writing H(1)(omega) = B, T(1)(omega) = X, L(1)(omega) = y # <> = <0|Y(B) * A_bar|0> + <0| (1 + L(0))[A_bar, X(B)}|0> # polar1 polar2 @@ -395,6 +410,14 @@ def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, m 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 inhomogeneous 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 homogeneous terms. + """ + solver_start = time.time() Dia = self.Dia @@ -514,10 +537,6 @@ def in_Y1(self, pertbar, X1, X2): t2 = self.ccwfn.t2 hbar = self.hbar L = self.H.L - - # Inhomogenous terms appearing in Y1 equations - #seems like these imhomogenous terms are computing at the beginning and not involve in the iteration itself - #may require moving to a sperate function # good r_Y1 = 2.0 * pertbar.Aov.copy() From 657e770091b72fbd8e81f7255b0ee4adb3908049 Mon Sep 17 00:00:00 2001 From: jattakumi Date: Wed, 17 Apr 2024 15:08:39 -0400 Subject: [PATCH 18/38] Fixed typos --- pycc/ccresponse.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index a386d8a..9f51250 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -297,7 +297,7 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis 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). The linear response function, <> generally reuires the following perturbed wave functions and frequencies: @@ -312,7 +312,9 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): ------ 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 @@ -410,13 +412,13 @@ def solve_right(self, pertbar, omega, e_conv=1e-12, r_conv=1e-12, maxiter=200, m 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 inhomogeneous 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 homogeneous terms. - """ + """ solver_start = time.time() From d04103d60250f2eec84de15447e9e564c1bb510a Mon Sep 17 00:00:00 2001 From: jattakumi Date: Thu, 2 May 2024 18:46:27 -0400 Subject: [PATCH 19/38] Merged future changes --- pycc/tests/ijk.dat | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 pycc/tests/ijk.dat diff --git a/pycc/tests/ijk.dat b/pycc/tests/ijk.dat new file mode 100644 index 0000000..80ce74f --- /dev/null +++ b/pycc/tests/ijk.dat @@ -0,0 +1,3 @@ +Total number of IJK combinations =: 35 +Num. of IJK with (Gi,Gj,Gk)=(0,0,0) =: 35 + thread 0: first_ijk=0, last_ijk=34 From cd10631dc0988ad907b264d165f894733951b8a0 Mon Sep 17 00:00:00 2001 From: jattakumi Date: Thu, 2 May 2024 19:04:15 -0400 Subject: [PATCH 20/38] Merging changes --- readthedocs.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/readthedocs.yml b/readthedocs.yml index 69d6db5..5f2da0e 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -1,12 +1,12 @@ # readthedocs.yml -version: 2 +version: 3 build: image: latest python: - version: 3.8 + version: 3.11 install: - method: pip path: . From 08c056c80f7012c497bd8a832b5872e1a7255b33 Mon Sep 17 00:00:00 2001 From: jattakumi Date: Thu, 2 May 2024 19:26:24 -0400 Subject: [PATCH 21/38] Merging changes --- .github/workflows/CI.yaml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index e23434d..069071b 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -24,6 +24,7 @@ jobs: os: [macOS-latest, ubuntu-latest] python-version: ['3.11'] + steps: - uses: actions/checkout@v3 @@ -34,6 +35,11 @@ jobs: df -h ulimit -a + - name: Set up Miniconda + uses: conda-incubator/setup-miniconda@v3 + with: + miniconda-version: "latest" + - name: Create Psi4 Environment uses: conda-incubator/setup-miniconda@v3 with: From 108c781d1f348de98f1ab381915169f49d1a9afa Mon Sep 17 00:00:00 2001 From: jattakumi Date: Thu, 2 May 2024 19:36:56 -0400 Subject: [PATCH 22/38] Merging changes --- .github/workflows/CI.yaml | 4 ---- pycc/tests/ijk.dat | 3 --- readthedocs.yml | 4 ++-- 3 files changed, 2 insertions(+), 9 deletions(-) delete mode 100644 pycc/tests/ijk.dat diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 069071b..68cea5c 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -35,10 +35,6 @@ jobs: df -h ulimit -a - - name: Set up Miniconda - uses: conda-incubator/setup-miniconda@v3 - with: - miniconda-version: "latest" - name: Create Psi4 Environment uses: conda-incubator/setup-miniconda@v3 diff --git a/pycc/tests/ijk.dat b/pycc/tests/ijk.dat deleted file mode 100644 index 80ce74f..0000000 --- a/pycc/tests/ijk.dat +++ /dev/null @@ -1,3 +0,0 @@ -Total number of IJK combinations =: 35 -Num. of IJK with (Gi,Gj,Gk)=(0,0,0) =: 35 - thread 0: first_ijk=0, last_ijk=34 diff --git a/readthedocs.yml b/readthedocs.yml index 5f2da0e..69d6db5 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -1,12 +1,12 @@ # readthedocs.yml -version: 3 +version: 2 build: image: latest python: - version: 3.11 + version: 3.8 install: - method: pip path: . From 9b1ffa9eccea72086e9ca875a799cc214539bc11 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 09:08:56 -0400 Subject: [PATCH 23/38] Adding some printing to the test case for additional error checking. --- pycc/tests/test_036_lr.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py index 3b2bf78..9770296 100644 --- a/pycc/tests/test_036_lr.py +++ b/pycc/tests/test_036_lr.py @@ -64,7 +64,11 @@ def test_linresp(): 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 From 0cb7b863803219f62d1b467f7e67ee939c195c5f Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 09:10:02 -0400 Subject: [PATCH 24/38] Limiting CI to test 36 for now. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 68cea5c..f5647bb 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest + coverage run -m pytest pycc/tests/test_036_lr.py coverage xml ls -l From 3619a345a5735eeed6019aa0a512bc7f59c253d1 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 09:16:44 -0400 Subject: [PATCH 25/38] Test 36 is fine, so adding all test back to workflow. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index f5647bb..68cea5c 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_036_lr.py + coverage run -m pytest coverage xml ls -l From 1b6e9aa17e90d2035d054b0edb891e0e6d8bb676 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:18:54 -0400 Subject: [PATCH 26/38] Test 36 dies when running all tests. Trying just test 35 before it. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 68cea5c..d3e0f64 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest + coverage run -m pytest test_035_eomccsd.py test_036_lr.py coverage xml ls -l From 2bd18a14f996c1fd879e20ed1376313c0a5eec37 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:21:13 -0400 Subject: [PATCH 27/38] This time with the paths correct. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index d3e0f64..f078dc1 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest test_035_eomccsd.py test_036_lr.py + coverage run -m pytest pycc/tests/test_035_eomccsd.py pycc/tests/test_036_lr.py coverage xml ls -l From cbb1af4c452ac0d0a8c4b16a0370f6aeaeb33a81 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:24:31 -0400 Subject: [PATCH 28/38] Test 36 dies when test 35 runs before it. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index f078dc1..f5647bb 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_035_eomccsd.py pycc/tests/test_036_lr.py + coverage run -m pytest pycc/tests/test_036_lr.py coverage xml ls -l From c8636b002b6b86762f67506a057a9c15d30538fe Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:27:35 -0400 Subject: [PATCH 29/38] Test 36 still runs correctly alone. Trying 37 after it. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index f5647bb..e17ccf0 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_036_lr.py + coverage run -m pytest pycc/tests/test_036_lr.py pycc/tests/test_037_rtcc3.py coverage xml ls -l From 1a642d73c3cadd405a7f4937660bee78a834483c Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:30:53 -0400 Subject: [PATCH 30/38] Trying test 2 in front of 36. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index e17ccf0..97d13fb 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_036_lr.py pycc/tests/test_037_rtcc3.py + coverage run -m pytest pycc/tests/test_002_ccsd_energy.py pycc/tests/test_036_lr.py coverage xml ls -l From a8806bccdb2439918b2ed38f7b0737886497ae71 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:36:23 -0400 Subject: [PATCH 31/38] Testing to see if 36 runs if I move it to the start of the test queue. --- .github/workflows/CI.yaml | 2 +- pycc/tests/{test_036_lr.py => test_001_lr.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename pycc/tests/{test_036_lr.py => test_001_lr.py} (100%) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 97d13fb..68cea5c 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_002_ccsd_energy.py pycc/tests/test_036_lr.py + coverage run -m pytest coverage xml ls -l diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_001_lr.py similarity index 100% rename from pycc/tests/test_036_lr.py rename to pycc/tests/test_001_lr.py From 0f86646cf0dbcbd182b5aaa2eaf2041b6c8abd94 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:48:00 -0400 Subject: [PATCH 32/38] Returning test to 36 and running 35 ahead of it. Trying a psi4.core.clean call at top of 36. --- .github/workflows/CI.yaml | 2 +- pycc/tests/{test_001_lr.py => test_036_lr.py} | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) rename pycc/tests/{test_001_lr.py => test_036_lr.py} (99%) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 68cea5c..97d13fb 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest + coverage run -m pytest pycc/tests/test_002_ccsd_energy.py pycc/tests/test_036_lr.py coverage xml ls -l diff --git a/pycc/tests/test_001_lr.py b/pycc/tests/test_036_lr.py similarity index 99% rename from pycc/tests/test_001_lr.py rename to pycc/tests/test_036_lr.py index 9770296..da75fdf 100644 --- a/pycc/tests/test_001_lr.py +++ b/pycc/tests/test_036_lr.py @@ -10,6 +10,7 @@ from ..data.molecules import * def test_linresp(): + psi4.core.clean() psi4.set_memory('2 GiB') psi4.core.set_output_file('output.dat', False) psi4.set_options({'basis': 'aug-cc-pvdz', From 221a74899b3e597ee0746547ad6b2a5dd4e2e703 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:50:45 -0400 Subject: [PATCH 33/38] clean made no difference. Trying to remove psi4.core.set_output_file in favor of psi4.set_output_file. --- pycc/tests/test_036_lr.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py index da75fdf..d7ce3bb 100644 --- a/pycc/tests/test_036_lr.py +++ b/pycc/tests/test_036_lr.py @@ -10,9 +10,8 @@ from ..data.molecules import * def test_linresp(): - psi4.core.clean() 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', 'e_convergence': 1e-12, From 159a374ba6d2dede7aff5bacb1a7756898564998 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:54:57 -0400 Subject: [PATCH 34/38] Frozen core gets turned on when another test runs before 36??? --- pycc/tests/test_036_lr.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py index d7ce3bb..4945c7a 100644 --- a/pycc/tests/test_036_lr.py +++ b/pycc/tests/test_036_lr.py @@ -14,6 +14,7 @@ def test_linresp(): psi4.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 From f4b85d55556cd51ef581b8040440e1fb0ace6b72 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 10:57:11 -0400 Subject: [PATCH 35/38] Fixed frozen core input change. Not clear why the option carried over from previous tests. --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 97d13fb..68cea5c 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_002_ccsd_energy.py pycc/tests/test_036_lr.py + coverage run -m pytest coverage xml ls -l From 9bbf3a60b7fd95bfb19f5a73241ce34f67bd20e4 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 11:14:52 -0400 Subject: [PATCH 36/38] Out of curiosity, I'm adding a psi4.core.clean_options to the start of test 36 to see if that fixes the problem generally. --- .github/workflows/CI.yaml | 2 +- pycc/tests/test_036_lr.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index 68cea5c..f078dc1 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest + coverage run -m pytest pycc/tests/test_035_eomccsd.py pycc/tests/test_036_lr.py coverage xml ls -l diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py index 4945c7a..4750880 100644 --- a/pycc/tests/test_036_lr.py +++ b/pycc/tests/test_036_lr.py @@ -10,11 +10,12 @@ from ..data.molecules import * def test_linresp(): + psi4.core.clean_options() psi4.set_memory('2 GiB') psi4.set_output_file('output.dat', False) psi4.set_options({'basis': 'aug-cc-pvdz', 'scf_type': 'pk', - 'freeze_core': 'false', +# 'freeze_core': 'false', 'e_convergence': 1e-12, 'd_convergence': 1e-12, 'r_convergence': 1e-12 From 9474d372de1a505c955041847bd00bee1586ac01 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 11:17:55 -0400 Subject: [PATCH 37/38] clean_options worked --- pycc/tests/test_036_lr.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py index 4750880..1539136 100644 --- a/pycc/tests/test_036_lr.py +++ b/pycc/tests/test_036_lr.py @@ -15,7 +15,6 @@ def test_linresp(): psi4.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 From e5a39e6dcab09c4ebd09b569eccb7902f3425a46 Mon Sep 17 00:00:00 2001 From: "T. Daniel Crawford" Date: Fri, 3 May 2024 11:18:26 -0400 Subject: [PATCH 38/38] running all tests --- .github/workflows/CI.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yaml b/.github/workflows/CI.yaml index f078dc1..68cea5c 100644 --- a/.github/workflows/CI.yaml +++ b/.github/workflows/CI.yaml @@ -70,7 +70,7 @@ jobs: run: | export KMP_DUPLICATE_LIB_OK=TRUE pip install coverage - coverage run -m pytest pycc/tests/test_035_eomccsd.py pycc/tests/test_036_lr.py + coverage run -m pytest coverage xml ls -l