Merge pull request #77 from zw0930/cc3-correct

Adding the missing term in T3 equation for RT-CC3

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
 

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])

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

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])

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
+
+

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)
 
@@ -240,6 +240,10 @@ def dipole(self, t1, t2, l1, l2, magnetic = 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)