IP/EA-ADC up to 3rd order

adcc/AdcMatrix.py

@@ -26,6 +26,8 @@
 
 from .LazyMp import LazyMp
 from .adc_pp import matrix as ppmatrix
+from .adc_ip import matrix as ipmatrix
+from .adc_ea import matrix as eamatrix
 from .timings import Timer, timed_member_call
 from .AdcMethod import AdcMethod
 from .functions import ones_like
@@ -73,14 +75,25 @@ class AdcMatrixlike:
 
 
 class AdcMatrix(AdcMatrixlike):
-    # Default perturbation-theory orders for the matrix blocks (== standard ADC-PP).
+    # Default perturbation-theory orders for the matrix blocks
+    # (== standard ADC-PP).
     default_block_orders = {
         #             ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None),
-        "adc0":  dict(ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None),  # noqa: E501
-        "adc1":  dict(ph_ph=1, ph_pphh=None, pphh_ph=None, pphh_pphh=None),  # noqa: E501
-        "adc2":  dict(ph_ph=2, ph_pphh=1,    pphh_ph=1,    pphh_pphh=0),     # noqa: E501
-        "adc2x": dict(ph_ph=2, ph_pphh=1,    pphh_ph=1,    pphh_pphh=1),     # noqa: E501
-        "adc3":  dict(ph_ph=3, ph_pphh=2,    pphh_ph=2,    pphh_pphh=1),     # noqa: E501
+        "adc0":     dict(ph_ph=0, ph_pphh=None, pphh_ph=None, pphh_pphh=None),  # noqa: E501
+        "adc1":     dict(ph_ph=1, ph_pphh=None, pphh_ph=None, pphh_pphh=None),  # noqa: E501
+        "adc2":     dict(ph_ph=2, ph_pphh=1,    pphh_ph=1,    pphh_pphh=0),     # noqa: E501
+        "adc2x":    dict(ph_ph=2, ph_pphh=1,    pphh_ph=1,    pphh_pphh=1),     # noqa: E501
+        "adc3":     dict(ph_ph=3, ph_pphh=2,    pphh_ph=2,    pphh_pphh=1),     # noqa: E501
+        "ip_adc0":  dict(h_h=0,   h_phh=None,   phh_h=None,   phh_phh=None),
+        "ip_adc1":  dict(h_h=1,   h_phh=None,   phh_h=None,   phh_phh=None),
+        "ip_adc2":  dict(h_h=2,   h_phh=1,      phh_h=1,      phh_phh=0),
+        "ip_adc2x": dict(h_h=2,   h_phh=1,      phh_h=1,      phh_phh=1),
+        "ip_adc3":  dict(h_h=3,   h_phh=2,      phh_h=2,      phh_phh=1),
+        "ea_adc0":  dict(p_p=0,   p_pph=None,   pph_p=None,   pph_pph=None),
+        "ea_adc1":  dict(p_p=1,   p_pph=None,   pph_p=None,   pph_pph=None),
+        "ea_adc2":  dict(p_p=2,   p_pph=1,      pph_p=1,      pph_pph=0),
+        "ea_adc2x": dict(p_p=2,   p_pph=1,      pph_p=1,      pph_pph=1),
+        "ea_adc3":  dict(p_p=3,   p_pph=2,      pph_p=2,      pph_pph=1),
     }
 
     def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
@@ -96,7 +109,8 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
             HF reference or MP ground state
         block_orders : optional
             The order of perturbation theory to employ for each matrix block.
-            If not set, defaults according to the selected ADC method are chosen.
+            If not set, defaults according to the selected ADC method are
+            chosen.
         intermediates : adcc.Intermediates or NoneType
             Allows to pass intermediates to re-use to this class.
         diagonal_precomputed: adcc.AmplitudeVector
@@ -123,6 +137,12 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
 
         self.timer = Timer()
         self.method = method
+        if "ip" in method.name:
+            self.type = "ip"
+        elif "ea" in method.name:
+            self.type = "ea"
+        else:
+            self.type = "pp"
         self.ground_state = hf_or_mp
         self.reference_state = hf_or_mp.reference_state
         self.mospaces = hf_or_mp.reference_state.mospaces
@@ -143,27 +163,35 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
             block_orders = tmp_orders
 
         # Sanity checks on block_orders
+        valid_blocks = {"pp": ("ph_ph", "ph_pphh", "pphh_ph", "pphh_pphh"),
+                        "ip": ("h_h", "h_phh", "phh_h", "phh_phh"),
+                        "ea": ("p_p", "p_pph", "pph_p", "pph_pph")}
         for block in block_orders.keys():
-            if block not in ("ph_ph", "ph_pphh", "pphh_ph", "pphh_pphh"):
+            if block not in valid_blocks[self.type]:
                 raise ValueError(f"Invalid block order key: {block}")
-        if block_orders["ph_pphh"] != block_orders["pphh_ph"]:
-            raise ValueError("ph_pphh and pphh_ph should always have "
-                             "the same order")
-        if block_orders["ph_pphh"] is not None \
-           and block_orders["pphh_pphh"] is None:
-            raise ValueError("pphh_pphh cannot be None if ph_pphh isn't.")
+        if (block_orders[valid_blocks[self.type][1]]
+                != block_orders[valid_blocks[self.type][2]]):
+            raise ValueError(f"{valid_blocks[self.type][1]} and "
+                             f"{valid_blocks[self.type][2]} should always "
+                             "have the same order")
+        if block_orders[valid_blocks[self.type][1]] is not None \
+                and block_orders[valid_blocks[self.type][3]] is None:
+            raise ValueError(f"{valid_blocks[self.type][3]} cannot be None if "
+                             f"{valid_blocks[self.type][1]} isn't.")
         self.block_orders = block_orders
 
         # Build the blocks and diagonals
         with self.timer.record("build"):
             variant = None
             if self.is_core_valence_separated:
                 variant = "cvs"
+            matrix = {"pp": ppmatrix, "ip": ipmatrix, "ea": eamatrix}
             blocks = {
-                block: ppmatrix.block(self.ground_state, block.split("_"),
-                                      order=order, intermediates=self.intermediates,
-                                      variant=variant)
-                for block, order in self.block_orders.items() if order is not None
+                block: matrix[self.type].block(
+                    self.ground_state, block.split("_"), order=order,
+                    intermediates=self.intermediates, variant=variant)
+                for block, order in self.block_orders.items()
+                if order is not None
             }
             # TODO Rename to self.block in 0.16.0
             self.blocks_ph = {bl: blocks[bl].apply for bl in blocks}

adcc/AdcMethod.py

@@ -24,10 +24,13 @@
 
 def get_valid_methods():
     valid_prefixes = ["cvs"]
-    valid_bases = ["adc0", "adc1", "adc2", "adc2x", "adc3"]
+    valid_bases = ["adc0", "adc1", "adc2", "adc2x", "adc3",
+                   "ip_adc0", "ip_adc1", "ip_adc2", "ip_adc2x", "ip_adc3",
+                   "ea_adc0", "ea_adc1", "ea_adc2", "ea_adc2x", "ea_adc3"]
 
+    # CVS-IP calculations not yet implemented
     ret = valid_bases + [p + "-" + m for p in valid_prefixes
-                         for m in valid_bases]
+                         for m in valid_bases[:5]]
     return ret
 
 
@@ -37,15 +40,15 @@ class AdcMethod:
     def __init__(self, method):
         if method not in self.available_methods:
             raise ValueError("Invalid method " + str(method) + ". Only "
-                             + ",".join(self.available_methods) + " are known.")
+                             + ",".join(self.available_methods)
+                             + " are known.")
 
         split = method.split("-")
-        self.__base_method = split[-1]
-        split = split[:-1]
+        self.__base_method = split.pop()
         self.is_core_valence_separated = "cvs" in split
 
         try:
-            if self.__base_method == "adc2x":
+            if self.__base_method in ["adc2x", "ip_adc2x", "ea_adc2x"]:
                 self.level = 2
             else:
                 self.level = int(self.__base_method[-1])
@@ -57,10 +60,13 @@ def at_level(self, newlevel):
         Return an equivalent method, where only the level is changed
         (e.g. calling this on a CVS method returns a CVS method)
         """
-        if self.is_core_valence_separated:
-            return AdcMethod("cvs-adc" + str(newlevel))
-        else:
-            return AdcMethod("adc" + str(newlevel))
+        try:
+            int(self.name[-1])
+        except ValueError:
+            raise NotImplementedError("This method is not implemented for ADC "
+                                      "methods that do not end with a number "
+                                      "like 'adc2x'.")
+        return AdcMethod(self.name[:-1] + str(newlevel))
 
     @property
     def name(self):

adcc/AmplitudeVector.py

@@ -27,13 +27,15 @@ class AmplitudeVector(dict):
     def __init__(self, *args, **kwargs):
         """
         Construct an AmplitudeVector. Typical use cases are
-        ``AmplitudeVector(ph=tensor_singles, pphh=tensor_doubles)``.
+        ``AmplitudeVector(ph=tensor_singles, pphh=tensor_doubles)``. For IP-ADC
+        ``AmplitudeVector(h=tensor_singles, phh=tensor_doubles)``, and for
+        EA-ADC ``AmplitudeVector(p=tensor_singles, pph=tensor_doubles)``
         """
         if args:
             warnings.warn("Using the list interface of AmplitudeVector is "
-                          "deprecated and will be removed in version 0.16.0. Use "
-                          "AmplitudeVector(ph=tensor_singles, pphh=tensor_doubles) "
-                          "instead.")
+                          "deprecated and will be removed in version 0.16.0. "
+                          "Use AmplitudeVector(ph=tensor_singles, "
+                          "pphh=tensor_doubles) instead.")
             if len(args) == 1:
                 super().__init__(ph=args[0])
             elif len(args) == 2:
@@ -54,11 +56,11 @@ def __setattr__(self, key, item):
     @property
     def blocks(self):
         warnings.warn("The blocks attribute will change behaviour in 0.16.0.")
-        if sorted(self.blocks_ph) == ["ph", "pphh"]:
+        if sorted(self.blocks_ph) == ["ph", "pphh"] or sorted(self.blocks_ph) == ["h", "phh"]:
             return ["s", "d"]
-        if sorted(self.blocks_ph) == ["pphh"]:
+        elif sorted(self.blocks_ph) == ["pphh"] or sorted(self.blocks_ph) == ["phh"]:
             return ["d"]
-        elif sorted(self.blocks_ph) == ["ph"]:
+        elif sorted(self.blocks_ph) == ["ph"] or sorted(self.blocks_ph) == ["h"]:
             return ["s"]
         elif sorted(self.blocks_ph) == []:
             return []
@@ -76,8 +78,8 @@ def blocks_ph(self):
     def __getitem__(self, index):
         if index in (0, 1, "s", "d"):
             warnings.warn("Using the list interface of AmplitudeVector is "
-                          "deprecated and will be removed in version 0.16.0. Use "
-                          "block labels like 'ph', 'pphh' instead.")
+                          "deprecated and will be removed in version 0.16.0. "
+                          "Use block labels like 'ph', 'pphh' instead.")
             if index in (0, "s"):
                 return self.__getitem__("ph")
             elif index in (1, "d"):
@@ -90,8 +92,8 @@ def __getitem__(self, index):
     def __setitem__(self, index, item):
         if index in (0, 1, "s", "d"):
             warnings.warn("Using the list interface of AmplitudeVector is "
-                          "deprecated and will be removed in version 0.16.0. Use "
-                          "block labels like 'ph', 'pphh' instead.")
+                          "deprecated and will be removed in version 0.16.0. "
+                          "Use block labels like 'ph', 'pphh' instead.")
             if index in (0, "s"):
                 return self.__setitem__("ph", item)
             elif index in (1, "d"):
@@ -144,7 +146,8 @@ def dot(self, other):
         if isinstance(other, list):
             # Make a list where the first index is all singles parts,
             # the second is all doubles parts and so on
-            return sum(self[b].dot([av[b] for av in other]) for b in self.keys())
+            return sum(self[b].dot([av[b] for av in other])
+                       for b in self.keys())
         else:
             return sum(self[b].dot(other[b]) for b in self.keys())
 
@@ -164,7 +167,8 @@ def __forward_to_blocks(self, fname, other):
             ret = {k: getattr(tensor, fname)(other[k])
                    for k, tensor in self.items()}
         else:
-            ret = {k: getattr(tensor, fname)(other) for k, tensor in self.items()}
+            ret = {k: getattr(tensor, fname)(other)
+                   for k, tensor in self.items()}
         if any(r == NotImplemented for r in ret.values()):
             return NotImplemented
         else:

adcc/ElectronicTransition.py

@@ -229,7 +229,8 @@ def rotatory_strength(self):
     @property
     @mark_excitation_property()
     def cross_section(self):
-        """List of one-photon absorption cross sections of all computed states"""
+        """List of one-photon absorption cross sections of all computed
+        states"""
         # TODO Source?
         fine_structure = constants.fine_structure
         fine_structure_au = 1 / fine_structure