Third-order ISR properties

adcc/ElectronicTransition.py

@@ -78,11 +78,11 @@ def __init__(self, data, method=None, property_method=None):
         if not isinstance(self.method, AdcMethod):
             self.method = AdcMethod(self.method)
         if property_method is None:
-            if self.method.level < 3:
-                property_method = self.method
-            else:
+            if self.method.level == 3:
                 # Auto-select ADC(2) properties for ADC(3) calc
                 property_method = self.method.at_level(2)
+            else:
+                property_method = self.method
         elif not isinstance(property_method, AdcMethod):
             property_method = AdcMethod(property_method)
         self._property_method = property_method
@@ -160,13 +160,12 @@ def transition_dipole_moment(self):
             warnings.warn("ADC(0) transition dipole moments are known to be "
                           "faulty in some cases.")
         dipole_integrals = self.operators.electric_dipole
-        return np.array([
+        ret = np.array([
             [product_trace(comp, tdm) for comp in dipole_integrals]
             for tdm in self.transition_dm
         ])
+        return ret
 
-    @cached_property
-    @mark_excitation_property()
     @timed_member_call(timer="_property_timer")
     def transition_dipole_moment_velocity(self):
         """List of transition dipole moments in the

adcc/LazyMp.py

@@ -93,6 +93,93 @@ def td2(self, space):
             - 0.5 * self.t2eri(b.oovv, b.oo)
         ) / denom
 
+    @cached_member_function
+    def tt2(self, space):
+        """Second order triples amplitudes"""
+        hf = self.reference_state
+        t2 = self.t2(b.oovv).evaluate()
+        denom = - direct_sum('ia,jkbc->ijkabc', self.df(b.ov),
+                             direct_sum('jb,kc->jkbc', self.df(b.ov),
+                             self.df(b.ov)))
+        amp = (
+            + einsum('idbc,jkad->ijkabc', hf.ovvv, t2)
+            + einsum('idab,jkcd->ijkabc', hf.ovvv, t2)
+            - einsum('jdab,ikcd->ijkabc', hf.ovvv, t2)
+            + einsum('kdab,ijcd->ijkabc', hf.ovvv, t2)
+            - einsum('jdbc,ikad->ijkabc', hf.ovvv, t2)
+            + einsum('kdbc,ijad->ijkabc', hf.ovvv, t2)
+            - einsum('idac,jkbd->ijkabc', hf.ovvv, t2)
+            + einsum('jdac,ikbd->ijkabc', hf.ovvv, t2)
+            - einsum('kdac,ijbd->ijkabc', hf.ovvv, t2)
+            + einsum('jkla,ilbc->ijkabc', hf.ooov, t2)
+            - einsum('ikla,jlbc->ijkabc', hf.ooov, t2)
+            + einsum('ijla,klbc->ijkabc', hf.ooov, t2)
+            + einsum('iklb,jlac->ijkabc', hf.ooov, t2)
+            - einsum('jklb,ilac->ijkabc', hf.ooov, t2)
+            - einsum('ijlb,klac->ijkabc', hf.ooov, t2)
+            - einsum('iklc,jlab->ijkabc', hf.ooov, t2)
+            + einsum('jklc,ilab->ijkabc', hf.ooov, t2)
+            + einsum('ijlc,klab->ijkabc', hf.ooov, t2)
+        )
+        return amp / denom
+
+    @cached_member_function
+    def ts3(self, space):
+        """Third order singles amplitudes"""
+        hf = self.reference_state
+        p0 = self.mp2_diffdm
+        td2 = self.td2(b.oovv)
+        tt2 = self.tt2(b.ooovvv)
+        denom = - self.df(b.ov)
+        amp = (
+            - einsum('jaib,jb->ia', hf.ovov, p0.ov)
+            + 0.5 * einsum('jkib,jkab->ia', hf.ooov, td2)
+            + 0.5 * einsum('jabc,ijbc->ia', hf.ovvv, td2)
+            + 0.25 * einsum('jkbc,ijkabc->ia', hf.oovv, tt2)
+        )
+        return amp / denom
+
+    @cached_member_function
+    def td3(self, space):
+        """Third order doubles amplitudes"""
+        hf = self.reference_state
+        p0 = self.mp2_diffdm
+        t2 = self.t2(b.oovv).evaluate()
+        td2 = self.td2(b.oovv).evaluate()
+        tt2 = self.tt2(b.ooovvv).evaluate()
+        t2eri_vv = einsum('klbd,klcd->bc', t2, hf.oovv).evaluate()
+        t2eri_oo = einsum('jlcd,klcd->jk', t2, hf.oovv).evaluate()
+        t2eri_oovv = einsum('jlbd,klcd->jkbc', t2, hf.oovv).evaluate()
+        denom = direct_sum('ia,jb->ijab', self.df(b.ov), self.df(b.ov))
+        amp = (
+            + 2 * einsum('jc,abic->ijab', p0.ov, hf.vvov).antisymmetrise(0, 1)
+            + 2 * einsum('kb,kaij->ijab', p0.ov, hf.ovoo).antisymmetrise(2, 3)
+            + 4 * einsum('ikac,kbjc->ijab', td2,
+                         hf.ovov).antisymmetrise(0, 1).antisymmetrise(2, 3)
+            - 0.5 * einsum('ijcd,abcd->ijab', td2, hf.vvvv)
+            - 0.5 * einsum('klab,klij->ijab', td2, hf.oooo)
+            + einsum('jklabc,klic->ijab', tt2, hf.ooov).antisymmetrise(0, 1)
+            + einsum('ijkbcd,kacd->ijab', tt2, hf.ovvv).antisymmetrise(2, 3)
+            - 0.25 * einsum('ijac,bc->ijab', t2, t2eri_vv)
+            - 0.25 * einsum('ijad,bd->ijab', t2, t2eri_vv)
+            + 0.25 * einsum('ijbc,ac->ijab', t2, t2eri_vv)
+            + 0.25 * einsum('ijbd,ad->ijab', t2, t2eri_vv)
+            + 0.25 * einsum('ijcd,klab,klcd->ijab', t2, t2, hf.oovv)
+            - 0.25 * einsum('ikab,jk->ijab', t2, t2eri_oo)
+            + 0.25 * einsum('ikac,jkbc->ijab', t2, t2eri_oovv)
+            + 0.25 * einsum('ikad,jkbd->ijab', t2, t2eri_oovv)
+            - 0.25 * einsum('ikbc,jkac->ijab', t2, t2eri_oovv)
+            - 0.25 * einsum('ikbd,jkad->ijab', t2, t2eri_oovv)
+            + 0.25 * einsum('il,jlab->ijab', t2eri_oo, t2)
+            - 0.25 * einsum('ilab,jl->ijab', t2, t2eri_oo)
+            + 0.25 * einsum('ilac,jlbc->ijab', t2, t2eri_oovv)
+            + 0.25 * einsum('ilad,jlbd->ijab', t2, t2eri_oovv)
+            - 0.25 * einsum('ilbc,jlac->ijab', t2, t2eri_oovv)
+            - 0.25 * einsum('ilbd,jlad->ijab', t2, t2eri_oovv)
+            + 0.25 * einsum('ik,jkab->ijab', t2eri_oo, t2)
+        )
+        return amp / denom
+
     @cached_member_function
     def t2eri(self, space, contraction):
         """
@@ -169,6 +256,34 @@ def mp2_diffdm(self):
         ret.reference_state = self.reference_state
         return evaluate(ret)
 
+    @cached_property
+    @timed_member_call(timer="timer")
+    def mp3_diffdm(self):
+        """
+        Return the MP3 difference density in the MO basis. mp2_diffdm is included
+        """
+        t2 = self.t2(b.oovv)
+        td2 = self.td2(b.oovv)
+        ts3 = self.ts3(b.ov)
+        tt2 = self.tt2(b.ooovvv)
+        p0 = self.mp2_diffdm
+        ret = self.mp2_diffdm
+
+        ret.oo = p0.oo - 0.5 * (
+            + einsum('jkab,ikab->ij', t2, td2)
+            + einsum('jkab,ikab->ij', td2, t2)
+        )
+        ret.ov = p0.ov + (
+            + ts3
+            - einsum('jb,ijab->ia', p0.ov, t2)
+            - 0.25 * einsum('jkbc,ijkabc->ia', t2, tt2)
+        )
+        ret.vv = p0.vv + 0.5 * (
+            + einsum('ijac,ijbc->ab', td2, t2)
+            + einsum('ijac,ijbc->ab', t2, td2)
+        )
+        return evaluate(ret)
+
     def density(self, level=2):
         """
         Return the MP density in the MO basis with all corrections
@@ -178,6 +293,8 @@ def density(self, level=2):
             return self.reference_state.density
         elif level == 2:
             return self.reference_state.density + self.mp2_diffdm
+        elif level == 3:
+            return self.reference_state.density + self.mp3_diffdm
         else:
             raise NotImplementedError("Only densities for level 1 and 2"
                                       " are implemented.")
@@ -191,6 +308,8 @@ def dipole_moment(self, level=2):
             return self.reference_state.dipole_moment
         elif level == 2:
             return self.mp2_dipole_moment
+        elif level == 3:
+            return self.mp3_dipole_moment
         else:
             raise NotImplementedError("Only dipole moments for level 1 and 2"
                                       " are implemented.")
@@ -274,6 +393,15 @@ def mp2_dipole_moment(self):
                              for comp in dipole_integrals])
         return refstate.dipole_moment + mp2corr
 
+    @cached_property
+    def mp3_dipole_moment(self):
+        # MP2_diffdm is included in mp3corr
+        refstate = self.reference_state
+        dipole_integrals = refstate.operators.electric_dipole
+        mp3corr = -np.array([product_trace(comp, self.mp3_diffdm)
+                            for comp in dipole_integrals])
+        return refstate.dipole_moment + mp3corr
+
 
 #
 # Register cvs_p0 intermediate

adcc/adc_pp/modified_transition_moments.py

@@ -60,6 +60,60 @@ def mtm_adc2(mp, dipop, intermediates):
     return AmplitudeVector(ph=f1, pphh=f2)
 
 
+def mtm_adc3(mp, dipop, intermediates):
+    second_order = mtm_adc2(mp, dipop, intermediates)
+    f1 = second_order.ph
+    f2 = second_order.pphh
+
+    t2 = mp.t2(b.oovv)
+    p0 = mp.mp2_diffdm
+    td2 = mp.td2(b.oovv)
+    tt2 = mp.tt2(b.ooovvv)
+    ts3 = mp.ts3(b.ov)
+    td3 = mp.td3(b.oovv)
+
+    f1 += (
+        + 1.0 * einsum('ib,ab->ia', ts3, dipop.vv)
+        - 1.0 * einsum('ja,ji->ia', ts3, dipop.oo)
+        - 1.0 * einsum('ijab,jb->ia', td3, dipop.ov)
+        + 1.0 * einsum('ijab,kb,jk->ia', t2, p0.ov, dipop.oo)
+        + 1.0 * einsum('jkab,kb,ji->ia', t2, p0.ov, dipop.oo)
+        - 1.0 * einsum('ijbc,jc,ab->ia', t2, p0.ov, dipop.vv)
+        - 1.0 * einsum('ijab,jc,cb->ia', t2, p0.ov, dipop.vv)
+        - 0.25 * einsum('ikbc,jkbc,ja->ia', t2, td2, dipop.ov)
+        - 0.25 * einsum('jkac,jkbc,ib->ia', t2, td2, dipop.ov)
+        + 0.5 * einsum('ijab,jkbc,kc->ia', t2, td2, dipop.ov)
+        - 0.25 * einsum('ikbc,jkbc,ja->ia', td2, t2, dipop.ov)
+        - 0.25 * einsum('jkac,jkbc,ib->ia', td2, t2, dipop.ov)
+        + 0.5 * einsum('ijab,jkbc,kc->ia', td2, t2, dipop.ov)
+        + 0.5 * einsum('ijkabc,jkcd,db->ia', tt2, t2, dipop.vv)
+        - 0.25 * einsum('ijkbcd,jkcd,ab->ia', tt2, t2, dipop.vv)
+        + 0.5 * einsum('ijkabc,klbc,jl->ia', tt2, t2, dipop.oo)
+        + 0.25 * einsum('jklabc,klbc,ji->ia', tt2, t2, dipop.oo)
+        + 0.25 * einsum('ijac,klbd,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        - 0.25 * einsum('ijad,klbc,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ijbd,klac,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        - 0.25 * einsum('ijcd,klab,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ikab,jlcd,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        - 0.25 * einsum('ikac,jlbd,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ikad,jlbc,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ikbc,jlad,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        - 0.25 * einsum('ikbd,jlac,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ikcd,jlab,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        - 0.25 * einsum('ilab,jkcd,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ilac,jkbd,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ilad,jkbc,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        - 0.25 * einsum('ilbc,jkad,klcd,jb->ia', t2, t2, t2, dipop.ov)
+        + 0.25 * einsum('ilbd,jkac,klcd,jb->ia', t2, t2, t2, dipop.ov)
+    )
+    f2 += (
+        - 1.0 * einsum('ijbc,ac->ijab', td2, dipop.vv).antisymmetrise(2, 3)
+        - 1.0 * einsum('ikab,kj->ijab', td2, dipop.oo).antisymmetrise(0, 1)
+        - 0.5 * einsum('ijkabc,kc->ijab', tt2, dipop.ov)
+    )
+    return AmplitudeVector(ph=f1, pphh=f2)
+
+
 def mtm_cvs_adc0(mp, dipop, intermediates):
     return AmplitudeVector(ph=dipop.cv)
 
@@ -78,6 +132,7 @@ def mtm_cvs_adc2(mp, dipop, intermediates):
     "adc0": mtm_adc0,
     "adc1": mtm_adc1,
     "adc2": mtm_adc2,
+    "adc3": mtm_adc3,
     "cvs-adc0": mtm_cvs_adc0,
     "cvs-adc1": mtm_cvs_adc0,  # Identical to CVS-ADC(0)
     "cvs-adc2": mtm_cvs_adc2,