diff --git a/pycc/ccresponse.py b/pycc/ccresponse.py index 9f51250..2b2a972 100644 --- a/pycc/ccresponse.py +++ b/pycc/ccresponse.py @@ -9,6 +9,7 @@ import time from .utils import helper_diis from .cclambda import cclambda +from .hamiltonian import Hamiltonian class ccresponse(object): """ @@ -295,6 +296,282 @@ def linresp(self, A, B, omega, e_conv=1e-13, r_conv=1e-13, maxiter=200, max_diis X1[X_key], X2[X_key], polar = self.solve_right(self.pertbar[pertkey], -omega, e_conv, r_conv, maxiter, max_diis, start_diis) check.append(polar) + def sym_linresp(self, pertkey_a, pertkey_b, X1_A, X2_A, X1_B, X2_B): + """ + Calculate the CC symmetric linear response function for polarizability at field-frequency omega(w1). + + The linear response function, <> generally requires 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. + """ + + # Please refer to eqn 94 of [Koch and Jørgensen, “Coupled Cluster Response Functions.”]. + # Writing H(1)(omega) = B, T^(1)(omega) = X, Lambda = L^(0) + # <> = <0|(1 + L^(0)) { [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)] + [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] } |0> + + contract = self.ccwfn.contract + l2 = self.cclambda.l2 + t2 = self.ccwfn.t2 + o = self.ccwfn.o + v = self.ccwfn.v + L = self.ccwfn.H.L + + polar1 = 0.0 + + polar1 += self.LCX(pertkey_a, X1_B, X2_B) + polar1 += self.LCX(pertkey_b, X1_A, X2_A) + + # <0|(HX1Y1)|0> + LHX1Y1 = 2.0 * contract('ijab, ia, jb', L[o, o, v, v], X1_B, X1_A) + Goo = contract('mjab,ijab->mi', t2, l2) + Gvv = -1.0 * contract('ijeb,ijab->ae', t2, l2) + r2_Gvv = contract('ae,ijeb->ijab', Gvv, L[o, o, v, v]) + r2_Goo = -1.0 * contract('mi,mjab->ijab', Goo, L[o, o, v, v]) + r2_Gvv = r2_Gvv + r2_Gvv.swapaxes(0, 1).swapaxes(2, 3) + r2_Goo = r2_Goo + r2_Goo.swapaxes(0, 1).swapaxes(2, 3) + LHX1Y1 += contract('ijab,ia,jb', r2_Gvv, X1_A, X1_B) # Gvv + LHX1Y1 += contract('ijab,ia,jb', r2_Goo, X1_A, X1_B) # Goo + polar1 += self.LHX1Y1(X1_B, X1_A) + polar1 += self.LHX1Y1(X1_A, X1_B) + + # <0|L2[[H, X2], Y2]|0> + temp = contract("ikac, ijab -> kjbc", L[o, o, v, v], X2_A) + temp = contract("kjbc, klcd -> jlbd", temp, X2_B) + LHX2Y2 = 2 * contract("jlbd, jlbd -> ", temp, l2) + polar1 += self.LHX2Y2(X2_A, X2_B) + polar1 += self.LHX2Y2(X2_B, X2_A) + + polar1 += self.LHX1Y2(X1_B, X2_A) + polar1 += self.LHX1Y2(X1_A, X2_B) + + polar1 += LHX1Y1 + LHX2Y2 + + return -1.0 * polar1 + + def LCX(self, pertkey_a, X1_B, X2_B): + """ + LCX: Is the first second term contributions to the linear response symmetric function. + <0|(1 + L^(0)){ [\bar{A}^(0), X^(1)(B)] + [\bar{B}^(1), X^(1)(B)]}|0> + """ + + contract = self.ccwfn.contract + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + + LCX = 0.0 + pertbar_A = self.pertbar[pertkey_a] + + Aov = pertbar_A.Aov + Aoo = pertbar_A.Aoo + Avv = pertbar_A.Avv + Avvvo = pertbar_A.Avvvo + Aovoo = pertbar_A.Aovoo + + # <0|[A_bar, X1]|0> + LCX += 2.0 * contract("ia, ia -> ", Aov, X1_B) + + # <0|L1 [A_bar, X1]|0> + temp_ov = contract('ab, ib -> ia', Avv, X1_B) + temp_ov -= contract('ji, ja -> ia', Aoo, X1_B) + LCX += contract('ia, ia', temp_ov, l1) + + temp_ov = -1.0 * contract('ja, ijab -> ib', Aov, X2_B) + temp_ov += 2.0 * contract('ja, jiab -> ib', Aov, X2_B) + LCX += contract("ib, ib -> ", temp_ov, l1) + + # <0|L2 [A_bar, X1]|0> + tmp = contract("ijbc, bcaj -> ia", l2, Avvvo) + LCX += contract("ia, ia -> ", tmp, X1_B) + tmp = contract("ijab, kbij -> ak", l2, Aovoo) + LCX -= 0.5 * contract("ak, ka -> ", tmp, X1_B) + tmp = contract("ijab, kaji -> bk", l2, Aovoo) + LCX -= 0.5 * contract("bk, kb -> ", tmp, X1_B) + + tmp = contract("ijab, kjab -> ik", l2, X2_B) + LCX -= 0.5 * contract("ik, ki -> ", tmp, Aoo) + tmp = contract("ijab, kiba-> jk", l2, X2_B) + LCX -= 0.5 * contract("jk, kj -> ", tmp, Aoo) + tmp = contract("ijab, ijac -> bc", l2, X2_B) + LCX += 0.5 * contract("bc, bc -> ", tmp, Avv) + tmp = contract("ijab, ijcb -> ac", l2, X2_B) + LCX += 0.5 * contract("ac, ac -> ", tmp, Avv) + + return LCX + + def LHX1Y1(self, X1_B, X1_A): + """ + LHX1Y1: Is a function for the second term (quadratic term) contribution to the linear response + function, where both perturbed amplitudes are first order. + # <0|(1 + L^(0)){ [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] }|0> + """ + + contract = self.ccwfn.contract + o = self.ccwfn.o + v = self.ccwfn.v + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.ccwfn.H.L + + LHX1Y1 = 0.0 + + # <0|L1[[H, X1], Y1]|0> + tmp = contract('ja, ia -> ij', hbar.Hov, X1_A) + tmp = contract('ij, jb -> ib', tmp, X1_B) + LHX1Y1 -= contract('ib, ib', tmp, l1) + + tmp = contract('jika, ia -> jk', 2 * hbar.Hooov, X1_B) + tmp = contract('jk, jb -> kb', tmp, X1_A) + LHX1Y1 -= contract('kb, kb', tmp, l1) + tmp = contract('jika, ia -> jk', hbar.Hooov.swapaxes(0, 1), X1_B) + tmp = contract('jk, jb -> kb', tmp, X1_A) + LHX1Y1 += contract('kb, kb', tmp, l1) + + tmp = contract('cjab, jb -> ac', 2 * hbar.Hvovv, X1_A) + tmp = contract('ac, ia -> ic', tmp, X1_B) + LHX1Y1 += contract('ic, ic -> ', tmp, l1) + tmp = contract('cjab, jb -> ac', hbar.Hvovv.swapaxes(2, 3), X1_A) + tmp = contract('ac, ia -> ic', tmp, X1_B) + LHX1Y1 -= contract('ic, ic -> ', tmp, l1) + + # <0|L2[[H, X1], Y1]|0> + temp = contract("jcka, ia -> ijkc", hbar.Hovov, X1_A) + temp = contract("ijkc, jb -> kibc", temp, X1_B) + LHX1Y1 -= contract("kibc, kibc -> ", temp, l2) + temp = contract("jcak, ia -> ijkc", hbar.Hovvo, X1_A) + temp = contract("ijkc, jb -> kicb", temp, X1_B) + LHX1Y1 -= contract("kicb, kicb -> ", temp, l2) + temp = contract('cdab, ia -> ibcd', hbar.Hvvvv, X1_B) + temp = contract('ibcd, jb -> ijcd', temp, X1_A) + LHX1Y1 += 0.5 * contract('ijcd, ijcd', temp, l2) + temp = contract('ijkl, ia -> klaj', hbar.Hoooo, X1_B) + temp = contract('klaj, jb -> klab', temp, X1_A) + LHX1Y1 += 0.5 * contract('klab, klab', temp, l2) + + return LHX1Y1 + + def LHX2Y2(self, X2_A, X2_B): + """ + LHX2Y2: Is a function for the second term (quadratic term) contribution to the linear response + function, where both perturbed amplitudes are second order. + # <0|(1 + L^(0)){ [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] }|0> + """ + + contract = self.ccwfn.contract + o = self.ccwfn.o + v = self.ccwfn.v + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + t2 = self.ccwfn.t2 + hbar = self.hbar + L = self.ccwfn.H.L + ERI = self.ccwfn.H.ERI + + LHX2Y2 = 0.0 + temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A) + temp = contract("bc, klcd -> klbd", temp, X2_B) + LHX2Y2 -= contract("klbd, klbd -> ", temp, l2) + temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A) + temp = contract("jk, jlbd -> klbd", temp, X2_B) + LHX2Y2 -= contract("klbd, klbd -> ", temp, l2) + temp = contract("ijac, jkbc -> ikab", L[o, o, v, v], X2_A) + temp = contract("ikab, ilad -> klbd", temp, X2_B) + LHX2Y2 -= contract("klbd, klbd -> ", temp, l2) + temp = contract("klab, ijab -> klij", ERI[o, o, v, v].copy(), X2_A) + temp = contract("klij, klcd -> ijcd", temp, X2_B) + LHX2Y2 += 0.5 * contract("ijcd, ijcd -> ", temp, l2) + temp = contract("klab, ikac -> ilbc", ERI[o, o, v, v].copy(), X2_A) + temp = contract("ilbc, jlbd -> ijcd", temp, X2_B) + LHX2Y2 += 0.5 * contract("ijcd, ijcd -> ", temp, l2) + temp = contract("klab, ilad -> ikbd", ERI[o, o, v, v].copy(), X2_A) + temp = contract("ikbd, jkbc -> ijcd", temp, X2_B) + LHX2Y2 += 0.5 * contract("ijcd, ijcd -> ", temp, l2) + + return LHX2Y2 + + def LHX1Y2(self, X1_B, X2_A): + """ + LHX1Y2: Is a function for the second term (quadratic term) contribution to the linear response + function, where perturbed amplitudes is a first order and a second order. + # <0|(1 + L^(0)){ [[\bar{B}^(0), X^(1)(B)], X^(1)(B)] }|0> + """ + + contract = self.ccwfn.contract + o = self.ccwfn.o + v = self.ccwfn.v + l1 = self.cclambda.l1 + l2 = self.cclambda.l2 + hbar = self.hbar + L = self.ccwfn.H.L + + LHX1Y2 = 0.0 + # + temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B) + temp = contract("jc, jkbc -> kb", temp, X2_A) + LHX1Y2 -= contract("kb, kb", temp, l1) + temp = contract("ijac, ia -> jc", L[o, o, v, v], X1_B) + temp = contract("jc, jkcb -> kb", temp, X2_A) + LHX1Y2 += 2.0 * contract("kb, kb", temp, l1) + temp = contract("ijac, ikac -> jk", L[o, o, v, v], X2_A) + temp = contract("jk, jb -> kb", temp, X1_B) + LHX1Y2 -= contract("kb, kb -> ", temp, l1) + temp = contract("ijac, ijab -> bc", L[o, o, v, v], X2_A) + temp = contract("bc, kc -> kb", temp, X1_B) + LHX1Y2 -= contract("kb, kb -> ", temp, l1) + + # + temp = contract('ja, ia -> ij', hbar.Hov, X1_B) + temp = contract('ij, jkbc -> ikbc', temp, X2_A) + LHX1Y2 -= contract('ikbc, ikbc', temp, l2) + temp = contract('ja, jb -> ab', hbar.Hov, X1_B) + temp = contract('ab, ikac -> ikbc', temp, X2_A) + LHX1Y2 -= contract('ikbc, ikbc', temp, l2) + + temp = contract('jima, ia -> jm', (2*hbar.Hooov - hbar.Hooov.swapaxes(0, 1)), X1_B) + temp = contract('jm, jkbc -> kmcb', temp, X2_A) + LHX1Y2 -= contract('kmcb, kmcb -> ', temp, l2) + temp = contract('jima, jb -> imab', (2*hbar.Hooov - hbar.Hooov.swapaxes(0, 1)), X1_B) + temp = contract('imab, ikac -> kmcb', temp, X2_A) + LHX1Y2 -= contract('kmcb, kmcb', temp, l2) + + temp = contract('djab, ia -> djib', (2*hbar.Hvovv - hbar.Hvovv.swapaxes(2, 3)), X1_B) + temp = contract('djib, jkbc -> kicd', temp, X2_A) + LHX1Y2 += contract('kicd, kicd', temp, l2) + temp = contract('djab, jb -> ad', (2 * hbar.Hvovv - hbar.Hvovv.swapaxes(2, 3)), X1_B) + temp = contract('ad, ikac -> kicd', temp, X2_A) + LHX1Y2 += contract('kicd, kicd', temp, l2) + + temp = contract('kdab, ia -> ikdb', hbar.Hvovv.swapaxes(0, 1), X1_B) + temp = contract('ikdb, jkbc -> ijdc', temp, X2_A) + LHX1Y2 -= contract('ijdc, ijdc', temp, l2) + temp = contract('kdab, jb -> kjad', hbar.Hvovv.swapaxes(0, 1), X1_B) + temp = contract('kjad, ikac -> ijdc', temp, X2_A) + LHX1Y2 -= contract('ijdc, ijdc', temp, l2) + temp = contract('kdab, kc -> abcd', hbar.Hvovv.swapaxes(0, 1), X1_B) + temp = contract('abcd, ijab -> ijdc', temp, X2_A) + LHX1Y2 -= contract('ijdc, ijdc', temp, l2) + + temp = contract('jkla, ia -> ijkl', hbar.Hooov, X1_B) + temp = contract('ijkl, jkbc -> libc', temp, X2_A) + LHX1Y2 += contract('libc, libc', temp, l2) + temp = contract('jkla, jb -> klab', hbar.Hooov, X1_B) + temp = contract('klab, ikac -> libc', temp, X2_A) + LHX1Y2 += contract('libc, libc -> ', temp, l2) + temp = contract('jkla, kc -> jlac', hbar.Hooov, X1_B) + temp = contract('jlac, ijab -> libc', temp, X2_A) + LHX1Y2 += contract('libc, libc -> ', temp, l2) + + + return LHX1Y2 def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): """ @@ -335,7 +612,7 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): polar1 += 0.5 * contract("baji, ijab -> ", Avvoo, Y2_B) # <0|[A_bar, X(B)]|0> polar2 += 2.0 * contract("ia, ia -> ", pertbar_A.Aov, X1_B) - # <0|L1(0) [A_bar, X2(B)]|0> + # <0|L1(0) [A_bar, X1(B)]|0> tmp = contract("ia, ic -> ac", l1, X1_B) polar2 += contract("ac, ac -> ", tmp, pertbar_A.Avv) tmp = contract("ia, ka -> ik", l1, X1_B) @@ -351,7 +628,7 @@ def linresp_asym(self, pertkey_a, X1_B, X2_B, Y1_B, Y2_B): polar2 -= 0.5 * contract("ak, ka -> ", tmp, X1_B) tmp = contract("ijab, kaji -> bk", l2, pertbar_A.Aovoo) polar2 -= 0.5 * contract("bk, kb -> ", tmp, X1_B) - # <0|L2(0)[A_bar, X1(B)]|0> + # <0|L2(0)[A_bar, X2(B)]|0> tmp = contract("ijab, kjab -> ik", l2, X2_B) polar2 -= 0.5 * contract("ik, ki -> ", tmp, pertbar_A.Aoo) tmp = contract("ijab, kiba-> jk", l2, X2_B) diff --git a/pycc/data/molecules.py b/pycc/data/molecules.py index 0ea22a1..d8a563c 100644 --- a/pycc/data/molecules.py +++ b/pycc/data/molecules.py @@ -265,6 +265,14 @@ symmetry c1 """ +h2o2 = """ + O -0.182400 -0.692195 -0.031109 + O 0.182400 0.692195 -0.031109 + H 0.533952 -1.077444 0.493728 + H -0.533952 1.077444 0.493728 + symmetry c1 + """ + moldict = {} moldict["He"] = he moldict["Be"] = be @@ -285,3 +293,4 @@ moldict["(S)-dimethylallene"] = sdma moldict["(S)-2-chloropropionitrile"] = s2cpn moldict["(R)-methylthiirane"] = rmt +moldict["hydrogenperoxide"] = h2o2 diff --git a/pycc/tests/test_036_lr.py b/pycc/tests/test_036_lr.py index 3b2bf78..17c61cd 100644 --- a/pycc/tests/test_036_lr.py +++ b/pycc/tests/test_036_lr.py @@ -10,17 +10,22 @@ from ..data.molecules import * def test_linresp(): + psi4.core.clean_options() psi4.set_memory('2 GiB') - psi4.core.set_output_file('output.dat', False) + psi4.set_output_file('output.dat', False) psi4.set_options({'basis': 'aug-cc-pvdz', 'scf_type': 'pk', + 'mp2_type': 'conv', + 'freeze_core': 'true', 'e_convergence': 1e-12, 'd_convergence': 1e-12, - 'r_convergence': 1e-12 + 'r_convergence': 1e-12, }) + mol = psi4.geometry(moldict["H2O"]) rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + # for t in threshold: e_conv = 1e-12 r_conv = 1e-12 @@ -64,7 +69,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 diff --git a/pycc/tests/test_038_lr_sym.py b/pycc/tests/test_038_lr_sym.py new file mode 100644 index 0000000..cf82784 --- /dev/null +++ b/pycc/tests/test_038_lr_sym.py @@ -0,0 +1,78 @@ +""" +Test CCSD linear response functions. +""" + +# Import package, test suite, and other packages as needed +import psi4 +import numpy as np +import pytest +import pycc +from ..data.molecules import * + +def test_sym_linresp(): + + psi4.core.clean_options() + psi4.set_memory('2 GiB') + psi4.core.set_output_file('output.dat', False) + psi4.set_options({'basis': 'aug-cc-pVDZ', + 'scf_type': 'pk', + 'mp2_type': 'conv', + 'freeze_core': 'true', + 'e_convergence': 1e-12, + 'd_convergence': 1e-12, + 'r_convergence': 1e-12}) + + mol = psi4.geometry(moldict["H2O"]) + rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True) + + e_conv = 1e-12 + r_conv = 1e-12 + + cc = pycc.ccwfn(rhf_wfn) + ecc = cc.solve_cc(e_conv, r_conv) + hbar = pycc.cchbar(cc) + cclambda = pycc.cclambda(cc, hbar) + lecc = cclambda.solve_lambda(e_conv, r_conv) + density = pycc.ccdensity(cc, cclambda) + + resp = pycc.ccresponse(density) + + omega1 = 0.0656 + + # Creating dictionaries + X_A = {} + X_B = {} + + for axis in range(0, 3): + string = "MU_" + resp.cart[axis] + A = resp.pertbar[string] + X_A[string] = resp.solve_right(A, omega1) + X_B[string] = resp.solve_right(A, -omega1) + + # Grabbing X, Y and declaring the matrix space for LR + polar_AB = np.zeros((3,3)) + + for a in range(0, 3): + string_a = "MU_" + resp.cart[a] + X1_A, X2_A, _ = X_A[string_a] + for b in range(0, 3): + string_b = "MU_" + resp.cart[b] + X_1B, X_2B, _ = X_B[string_b] + polar_AB[a,b] = resp.sym_linresp(string_a, string_b, X1_A, X2_A, X_1B, X_2B) + + print(f"Dynamic Polarizability Tensor @ w = {omega1} a.u.:") + print(polar_AB) + print("Average Dynamic Polarizability:") + polar_AB_avg = np.average([polar_AB[0,0], polar_AB[1,1], polar_AB[2,2]]) + print(polar_AB_avg) + + # Validating from psi4 + polar_xx = 9.932240347101651 + polar_yy = 13.446487681337629 + polar_zz = 11.344346098120035 + polar_avg = 11.574358042186 + + assert(abs(polar_AB[0,0] - polar_xx) < 1e-7) + assert(abs(polar_AB[1,1] - polar_yy) < 1e-7) + assert(abs(polar_AB[2,2] - polar_zz) < 1e-7) + assert(abs(polar_AB_avg - polar_avg) < 1e-7)