diff --git a/adcc/AdcMatrix.py b/adcc/AdcMatrix.py
index 54bec069..bdd414e4 100644
--- a/adcc/AdcMatrix.py
+++ b/adcc/AdcMatrix.py
@@ -28,7 +28,7 @@
 from .LazyMp import LazyMp
 from .adc_pp import matrix as ppmatrix
 from .timings import Timer, timed_member_call
-from .AdcMethod import AdcMethod
+from .AdcMethod import AdcMethod, IsrMethod, Method
 from .functions import ones_like
 from .Intermediates import Intermediates
 from .AmplitudeVector import AmplitudeVector
@@ -73,17 +73,22 @@ class AdcMatrixlike:
 
     _special_block_orders = {
         "adc2x": {"ph_ph": 2, "ph_pphh": 1, "pphh_ph": 1, "pphh_pphh": 1},
+        "isr1s": {"ph_ph": 1, "ph_pphh": None, "pphh_ph": None, "pphh_pphh": None},
     }
 
     @classmethod
-    def _default_block_orders(cls, method: AdcMethod) -> dict[str, int]:
+    def _default_block_orders(cls,
+                              method: Method
+                              ) -> dict[str, int]:
         """
         Determines the default block orders for the given adc method.
         """
         # check if we have a special method like adc2x
         # I guess base_method should also contain the adc_type prefix so
         # we don't need to separate different adc_types
-        block_orders = cls._special_block_orders.get(method.base_method.name, None)
+        block_orders = cls._special_block_orders.get(
+            method.base_method.name, None
+        )
         if block_orders is not None:
             return block_orders.copy()
         # otherwise assume that we have a "normal" PP/IP/...-ADC(n) method
@@ -96,9 +101,23 @@ def _default_block_orders(cls, method: AdcMethod) -> dict[str, int]:
             raise ValueError(f"Unknown adc type {method.adc_type} for method "
                              f"{method.name}. Can not determine default block "
                              "orders.")
-        spaces = [
-            "p" * i + min_space + "h" * i for i in range(0, (method.level // 2) + 1)
-        ]
+        # ADC matrices have first-order coupling whereas ISR matrices
+        # have a zeroth-order coupling between adjacent excitation classes
+        # see https://doi.org/10.1063/1.1752875
+        if isinstance(method, AdcMethod):
+            # First-order coupling between adjacent excitation classes.
+            spaces = [
+                "p" * i + min_space + "h" * i
+                for i in range(0, (method.level.to_int() // 2) + 1)
+            ]
+        elif isinstance(method, IsrMethod):
+            # Zeroth-order coupling between adjacent excitation classes.
+            spaces = [
+                "p" * i + min_space + "h" * i
+                for i in range(0, ((method.level.to_int() + 1) // 2) + 1)
+            ]
+        else:
+            raise ValueError(f"Invalid method: {method.name}")
         # exploit the fact that the spaces are sorted from small to high:
         # If we walk the adc matrix in any direction we always have to subtract 1!
         # Therefore, we can determine the order according to the position of the
@@ -107,26 +126,27 @@ def _default_block_orders(cls, method: AdcMethod) -> dict[str, int]:
         ret = {}
         for ((i1, bra), (i2, ket)) in \
                 itertools.product(enumerate(spaces), repeat=2):
-            order = method.level - i1 - i2
-            assert order >= 0
+            order = method.level.to_int() - i1 - i2
+            # For ISR matrices allow missing diagonal blocks.
+            order = None if order < 0 else order
             ret[f"{bra}_{ket}"] = order
         return ret
 
     @classmethod
     def _validate_block_orders(cls, block_orders: dict[str, int],
-                               method: AdcMethod,
+                               method: Method,
                                allow_missing_diagonal_blocks: bool = False) -> None:
         """
-        Validates that the given block_orders form a valid adc matrix for the given
-        adc method.
+        Validates that the given block_orders form a valid adc/isr matrix for the
+        given adc(isr) method.
 
         Parameters
         ----------
         block_orders: dict[str, int]
             The block orders to validate. Block orders should be of the form
             {'ph_ph': 2, 'ph_pphh': 1, ...}
-        method: AdcMethod
-            The adc method/adc type (PP-ADC, ...) for which to validate
+        method: Method
+            The adc/isr method/adc type (PP-ADC/ISR, ...) for which to validate
             the block_orders.
         allow_missing_diagonal_blocks: bool, optional
             If set, couplings between missing diagonal blocks are allowed, e.g.,
@@ -137,6 +157,7 @@ def _validate_block_orders(cls, block_orders: dict[str, int],
         for block, order in block_orders.items():
             if order is None:
                 continue
+            assert order >= 0
             # ensure that the block is valid for the given adc type
             bra, ket = block.split("_")
             if not cls._is_valid_space(bra, method) or \
@@ -164,11 +185,11 @@ def _validate_block_orders(cls, block_orders: dict[str, int],
                                  f"{ket_diag} are in the matrix too.")
 
     @classmethod
-    def _is_valid_space(cls, space: str, method: AdcMethod) -> bool:
+    def _is_valid_space(cls, space: str, method: Method) -> bool:
         """
         Checks whether the given space ('ph' for instance) is valid for the given
         adc method. Thereby we only verify that the space matches the adc_type of
-        adc method!
+        method!
         """
         n_particle, n_hole = space.count("p"), space.count("h")
         # ensure that the space is of the form pp...hh...
@@ -236,11 +257,7 @@ def __init__(self, method, hf_or_mp, block_orders=None, intermediates=None,
             self.intermediates = Intermediates(self.ground_state)
 
         self.block_orders = self._default_block_orders(self.method)
-        if block_orders is None:
-            if method.level > 3:
-                raise NotImplementedError("The ADC secular matrix is not "
-                                          f"implemented for method {method.name}.")
-        else:
+        if block_orders is not None:
             self.block_orders.update(block_orders)
         self._validate_block_orders(
             block_orders=self.block_orders, method=self.method,
diff --git a/adcc/AdcMethod.py b/adcc/AdcMethod.py
index d1781e28..af308d1c 100644
--- a/adcc/AdcMethod.py
+++ b/adcc/AdcMethod.py
@@ -20,76 +20,125 @@
 ## along with adcc. If not, see <http://www.gnu.org/licenses/>.
 ##
 ## ---------------------------------------------------------------------
+from typing import Optional, TypeVar
+from enum import Enum
 
+T = TypeVar("T", bound="Method")
 
-def get_valid_methods():
-    valid_prefixes = ["cvs"]
-    valid_bases = ["adc0", "adc1", "adc2", "adc2x", "adc3"]
 
-    ret = valid_bases + [p + "-" + m for p in valid_prefixes
-                         for m in valid_bases]
-    return ret
+class MethodLevel(Enum):
+    # numeric levels
+    ZERO = 0
+    ONE = 1
+    TWO = 2
+    THREE = 3
+    FOUR = 4
+    FIVE = 5
 
+    # special levels
+    TWO_X = "2x"
+    ONE_S = "1s"
+    THREE_D = "3d"
 
-class AdcMethod:
-    available_methods = get_valid_methods()
+    def to_str(self) -> str:
+        return str(self.value)
 
-    def __init__(self, method):
-        if method not in self.available_methods:
-            raise ValueError("Invalid method " + str(method) + ". Only "
-                             + ",".join(self.available_methods) + " are known.")
+    def to_int(self) -> int:
+        # numerical methods
+        if isinstance(self.value, int):
+            return self.value
+        # return base int for special methods
+        elif isinstance(self.value, str):
+            return int(self.value[0])
+        else:
+            raise ValueError
+
+
+class Method:
+    # this has to be set on the child classes
+    _method_base_name: Optional[str] = None
+    max_level: int = 0
+    special_levels: tuple[MethodLevel, ...] = tuple()
 
+    def __init__(self, method: str):
+        assert self._method_base_name is not None
+
+        # validate base method type
         split = method.split("-")
-        self.__base_method = split[-1]
+        if not split[-1].startswith(self._method_base_name):
+            raise ValueError(f"{split[-1]} is not a valid method type")
+
+        # validate method level
+        level = split[-1][len(self._method_base_name):]
+        if level.isnumeric():
+            self.level: MethodLevel = MethodLevel(int(level))
+        else:
+            self.level: MethodLevel = MethodLevel(level)
+        self._validate_level(self.level)
+
+        assert self._base_method == split[-1]
+
+        # validate prefix
         split = split[:-1]
-        self.is_core_valence_separated = "cvs" in split
+        if split and split[0] not in ["cvs"]:
+            raise ValueError(f"{split[0]} is not a valid method prefix")
+
+        self.is_core_valence_separated: bool = "cvs" in split
         # NOTE: added this to make the testdata generation ready for IP/EA
-        self.adc_type = "pp"
+        self.adc_type: str = "pp"
 
-        try:
-            if self.__base_method == "adc2x":
-                self.level = 2
-            else:
-                self.level = int(self.__base_method[-1])
-        except ValueError:
-            raise ValueError("Not a valid base method: " + self.__base_method)
+    def _validate_level(self, level: MethodLevel) -> None:
+        if isinstance(level.value, int):
+            if level.value <= self.max_level:
+                return
 
-    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))
+        # special cases
+        if level in self.special_levels:
+            return
+
+        raise NotImplementedError(f"{self._base_method} is not implemented.")
 
     @property
-    def name(self):
+    def name(self) -> str:
+        """The name of the Method as string."""
         if self.is_core_valence_separated:
-            return "cvs-" + self.__base_method
+            return "cvs-" + self._base_method
         else:
-            return self.__base_method
+            return self._base_method
 
     @property
-    def property_method(self):
-        """
-        The name of the canonical method to use for computing properties
-        for this ADC method. This only differs from the name property
-        for the ADC(2)-x family of methods.
-        """
-        if self.__base_method == "adc2x":
-            return AdcMethod(self.name.replace("adc2x", "adc2")).name
-        else:
-            return self.name
+    def _base_method(self) -> str:
+        return self._method_base_name + self.level.to_str()
 
     @property
-    def base_method(self):
+    def base_method(self: T) -> T:
         """
         The base (full) method, i.e. with all approximations such as
         CVS stripped off.
         """
-        return AdcMethod(self.__base_method)
+        return self.__class__(self._base_method)
+
+    def at_level(self: T, newlevel: int) -> T:
+        """
+        Return an equivalent method, where only the level is changed
+        (e.g. calling this on a CVS method returns a CVS method)
+        """
+        assert self._method_base_name is not None
+        if self.is_core_valence_separated:
+            return self.__class__("cvs-" + self._method_base_name + str(newlevel))
+        else:
+            return self.__class__(self._method_base_name + str(newlevel))
+
+    def as_method(self, method_cls: type[T]) -> T:
+        """
+        Return a equivalent Method with the method base name replaced
+        by the provided name.
+        """
+        assert self._method_base_name is not None
+        assert method_cls._method_base_name is not None
+        return method_cls(
+            self.name.replace(self._method_base_name, method_cls._method_base_name)
+        )
 
     def __eq__(self, other):
         return self.name == other.name
@@ -99,3 +148,15 @@ def __ne__(self, other):
 
     def __repr__(self):
         return "Method(name={})".format(self.name)
+
+
+class AdcMethod(Method):
+    _method_base_name = "adc"
+    max_level = 3
+    special_levels = (MethodLevel.TWO_X,)
+
+
+class IsrMethod(Method):
+    _method_base_name = "isr"
+    max_level = 2
+    special_levels = (MethodLevel.ONE_S,)
diff --git a/adcc/ElectronicStates.py b/adcc/ElectronicStates.py
index 0ff92f44..82d6fa59 100644
--- a/adcc/ElectronicStates.py
+++ b/adcc/ElectronicStates.py
@@ -3,7 +3,7 @@
 import warnings
 
 from .AdcMatrix import AdcMatrix
-from .AdcMethod import AdcMethod
+from .AdcMethod import AdcMethod, IsrMethod, MethodLevel
 from .AmplitudeVector import AmplitudeVector
 from .FormatDominantElements import FormatDominantElements
 from .FormatIndex import (
@@ -83,14 +83,16 @@ def __init__(self, data, method: str = None,
             self.method: AdcMethod = AdcMethod(self.method)
 
         if property_method is None:
-            if self.method.level < 3:
-                property_method = self.method
+            if self.method.level in [MethodLevel.TWO_X, MethodLevel.THREE]:
+                # Auto-select ISR(2) properties for ADC(2)-x and ADC(3) calc
+                warnings.warn(f"ISR({self.method.level.to_str()}) not implemented."
+                              f" Property method is selected as ISR(2).")
+                property_method = self.method.at_level(2).as_method(IsrMethod)
             else:
-                # Auto-select ADC(2) properties for ADC(3) calc
-                property_method = self.method.at_level(2)
-        elif not isinstance(property_method, AdcMethod):
-            property_method = AdcMethod(property_method)
-        self._property_method: AdcMethod = property_method
+                property_method = self.method.as_method(IsrMethod)
+        elif not isinstance(property_method, IsrMethod):
+            property_method = IsrMethod(property_method)
+        self._property_method: IsrMethod = property_method
 
         # Special stuff for special solvers
         if isinstance(data, EigenSolverStateBase):
@@ -138,8 +140,8 @@ def size(self) -> int:
         return self._excitation_energy.size
 
     @property
-    def property_method(self) -> AdcMethod:
-        """The method used to evaluate ADC properties"""
+    def property_method(self) -> IsrMethod:
+        """The method used to evaluate ISR properties"""
         return self._property_method
 
     @property
@@ -187,7 +189,7 @@ def state_dm(self) -> list[OneParticleDensity]:
 
     def _state_dm(self, state_n: int) -> OneParticleDensity:
         """List of state density matrices of all computed states"""
-        mp_density = self.ground_state.density(self.property_method.level)
+        mp_density = self.ground_state.density(self.property_method.level.to_int())
         diffdm = self._state_diffdm(state_n)
         return mp_density + diffdm
 
@@ -213,7 +215,9 @@ def state_dm_2p(self) -> list[TwoParticleDensity]:
 
     def _state_dm_2p(self, state_n: int) -> TwoParticleDensity:
         """List of two particle state density matrices of all computed states"""
-        mp_density = self.ground_state.density_2p(self.property_method.level)
+        mp_density = self.ground_state.density_2p(
+            self.property_method.level.to_int()
+        )
         diffdm = self._state_diffdm_2p(state_n)
         return mp_density + diffdm
 
@@ -226,10 +230,10 @@ def state_ssq(self) -> np.ndarray:
     def _state_ssq(self, state_n: int) -> float:
         """Computes the <S^2> of a single state."""
         pmethod = self.property_method
-        if pmethod.level == 0:
+        if pmethod.level.to_int() == 0:
             gs_ssq = self.reference_state.ssq
         else:
-            gs_ssq = self.ground_state.ssq(pmethod.level)
+            gs_ssq = self.ground_state.ssq(pmethod.level.to_int())
 
         ssq_1p_op = self.operators.ssq_1p
         ssq_2p_op = self.operators.ssq_2p
@@ -248,10 +252,10 @@ def state_dipole_moment(self) -> np.ndarray:
     def _state_dipole_moment(self, state_n: int) -> np.ndarray:
         """Computes the state dipole moment for a single state"""
         pmethod = self.property_method
-        if pmethod.level == 0:
+        if pmethod.level.to_int() == 0:
             gs_dip_moment = self.reference_state.dipole_moment
         else:
-            gs_dip_moment = self.ground_state.dipole_moment(pmethod.level)
+            gs_dip_moment = self.ground_state.dipole_moment(pmethod.level.to_int())
 
         dipole_integrals = self.operators.electric_dipole
         ddm = self._state_diffdm(state_n)
diff --git a/adcc/IsrMatrix.py b/adcc/IsrMatrix.py
index 82dfd09f..8c5e823a 100644
--- a/adcc/IsrMatrix.py
+++ b/adcc/IsrMatrix.py
@@ -22,8 +22,10 @@
 ## ---------------------------------------------------------------------
 import libadcc
 
+from itertools import product
+
 from .AdcMatrix import AdcMatrixlike
-from .AdcMethod import AdcMethod
+from .AdcMethod import IsrMethod
 from .adc_pp import bmatrix as ppbmatrix
 from .AmplitudeVector import AmplitudeVector
 from .LazyMp import LazyMp
@@ -37,11 +39,11 @@ class IsrMatrix(AdcMatrixlike):
     def __init__(self, method, hf_or_mp, operator, block_orders=None):
         """
         Initialise an ISR matrix of a given one-particle operator
-        for the provided ADC method.
+        for the provided ISR method.
 
         Parameters
         ----------
-        method : str or adcc.AdcMethod
+        method : str or adcc.IsrMethod
             Method to use.
         hf_or_mp : adcc.ReferenceState or adcc.LazyMp
             HF reference or MP ground state.
@@ -60,8 +62,8 @@ def __init__(self, method, hf_or_mp, operator, block_orders=None):
                             "either a LazyMp, a ReferenceState or a "
                             "HartreeFockSolution_i.")
 
-        if not isinstance(method, AdcMethod):
-            method = AdcMethod(method)
+        if not isinstance(method, IsrMethod):
+            method = IsrMethod(method)
 
         if isinstance(operator, (list, tuple)):
             self.operator = tuple(operator)
@@ -84,8 +86,7 @@ def __init__(self, method, hf_or_mp, operator, block_orders=None):
         self.block_orders = self._default_block_orders(self.method)
         if block_orders is None:
             # only implemented through PP-ADC(2)
-            if method.adc_type != "pp" or method.name.endswith("adc2x") \
-                    or method.level > 2:
+            if method.adc_type != "pp":
                 raise NotImplementedError("The B-matrix is not implemented "
                                           f"for method {method.name}.")
         else:
@@ -100,15 +101,15 @@ def __init__(self, method, hf_or_mp, operator, block_orders=None):
             variant = None
             if self.is_core_valence_separated:
                 variant = "cvs"
-            blocks = [{
+            blocks = tuple({
                 block: ppbmatrix.block(self.ground_state, op,
                                        block.split("_"), order=order,
                                        variant=variant)
                 for block, order in self.block_orders.items() if order is not None
-            } for op in self.operator]
-            self.blocks = [{
+            } for op in self.operator)
+            self.blocks = tuple({
                 b: bl[b].apply for b in bl
-            } for bl in blocks]
+            } for bl in blocks)
 
     @timed_member_call()
     def matvec(self, v):
@@ -119,9 +120,17 @@ def matvec(self, v):
         If a list of OneParticleOperator objects was passed to the class
         instantiation operator, a list of AmplitudeVector objects is returned.
         """
+        # Check which blocks are present in the AmplitudeVector.
+        # Missing blocks can be treated as zero vectors, so their
+        # contribution to the matvec is zero and does not need to be computed.
+        avail_blocks = [f"{bra}_{ket}" for bra, ket in product(v.blocks, repeat=2)]
+        filtered_blocks = [
+            {blk: v for blk, v in comp.items() if blk in avail_blocks}
+            for comp in self.blocks
+        ]
         ret = [
             sum(block(v) for block in bl_ph.values())
-            for bl_ph in self.blocks
+            for bl_ph in filtered_blocks
         ]
         if len(ret) == 1:
             return ret[0]
@@ -139,7 +148,7 @@ def rmatvec(self, v):
                 AmplitudeVector(ph=-1.0 * mv.ph, pphh=-1.0 * mv.pphh)
                 for mv in self.matvec(v)
             ]
-            # operators without any symmetry
+        # operators without any symmetry
         else:
             return NotImplemented
 
diff --git a/adcc/__init__.py b/adcc/__init__.py
index 004e66ff..c4ac3f13 100644
--- a/adcc/__init__.py
+++ b/adcc/__init__.py
@@ -29,7 +29,7 @@
 from .Symmetry import Symmetry
 from .MoSpaces import MoSpaces
 from .AdcMatrix import AdcMatrix
-from .AdcMethod import AdcMethod
+from .AdcMethod import AdcMethod, IsrMethod
 from .functions import (copy, direct_sum, dot, einsum, empty_like,
                         evaluate, lincomb, linear_combination, nosym_like,
                         ones_like, transpose, zeros_like)
@@ -55,7 +55,7 @@
 from .exceptions import InputError
 
 __all__ = ["run_adc", "InputError", "AdcMatrix",
-           "AdcMethod", "Symmetry", "ReferenceState", "MoSpaces",
+           "AdcMethod", "IsrMethod", "Symmetry", "ReferenceState", "MoSpaces",
            "einsum", "copy", "dot", "empty_like", "evaluate",
            "lincomb", "nosym_like", "ones_like", "transpose",
            "linear_combination", "zeros_like", "direct_sum",
@@ -93,7 +93,7 @@ def adc0(*args, **kwargs):
 @with_runadc_doc
 def cis(*args, **kwargs):
     state = run_adc(*args, **kwargs, method="adc1")
-    return ExcitedStates(state, property_method="adc0")
+    return ExcitedStates(state, property_method="isr0")
 
 
 @with_runadc_doc
diff --git a/adcc/adc_pp/bmatrix.py b/adcc/adc_pp/bmatrix.py
index 4177e34e..55262759 100644
--- a/adcc/adc_pp/bmatrix.py
+++ b/adcc/adc_pp/bmatrix.py
@@ -35,9 +35,9 @@
 
 """
 `apply` is a function mapping an AmplitudeVector to the contribution of this
-block to the result of applying the ADC matrix.
+block to the result of applying the ISR matrix.
 """
-AdcBlock = namedtuple("AdcBlock", ["apply"])
+IsrBlock = namedtuple("IsrBlock", ["apply"])
 
 
 def block(ground_state, operator, spaces, order, variant=None):
@@ -84,7 +84,7 @@ def apply(ampl):
             + 1.0 * einsum('ic,ac->ia', ampl.ph, op.vv)
             - 1.0 * einsum('ka,ki->ia', ampl.ph, op.oo)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
 
 
 def block_pphh_pphh_0(ground_state, op):
@@ -99,7 +99,7 @@ def apply(ampl):
                 + 2.0 * einsum('kiab,kj->ijab', ampl.pphh, op.oo)
             ).antisymmetrise(0, 1)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
 
 
 #
@@ -111,7 +111,7 @@ def apply(ampl):
             - 2.0 * einsum('ilad,ld->ia', ampl.pphh, op.ov)
             + 2.0 * einsum('ilca,lc->ia', ampl.pphh, op.ov)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
 
 
 def block_pphh_ph_0(ground_state, op):
@@ -124,7 +124,7 @@ def apply(ampl):
                 - 1.0 * einsum('jb,ai->ijab', ampl.ph, op.vo)
             ).antisymmetrise(0, 1).antisymmetrise(2, 3)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
 
 
 #
@@ -150,7 +150,7 @@ def apply(ampl):
             - 2.0 * einsum('klad,kled,ei->ia', ampl.pphh, t2, op.vo)
             - 2.0 * einsum('ilcd,nlcd,an->ia', ampl.pphh, t2, op.vo)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
 
 
 def block_pphh_ph_1(ground_state, op):
@@ -179,7 +179,7 @@ def apply(ampl):
                 + 1.0 * einsum('jc,niab,nc->ijab', ampl.ph, t2, op.ov)
             ).antisymmetrise(0, 1)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
 
 
 #
@@ -221,4 +221,4 @@ def apply(ampl):
             - 1.0 * einsum('kc,kncf,imaf,mn->ia', ampl.ph, t2, t2, op.oo)
             + 1.0 * einsum('kc,knce,inaf,ef->ia', ampl.ph, t2, t2, op.vv)
         ))
-    return AdcBlock(apply)
+    return IsrBlock(apply)
diff --git a/adcc/adc_pp/modified_transition_moments.py b/adcc/adc_pp/modified_transition_moments.py
index ade7a68a..39b86b97 100644
--- a/adcc/adc_pp/modified_transition_moments.py
+++ b/adcc/adc_pp/modified_transition_moments.py
@@ -24,28 +24,28 @@
 
 from adcc import block as b
 from adcc.LazyMp import LazyMp
-from adcc.AdcMethod import AdcMethod
+from adcc.AdcMethod import IsrMethod
 from adcc.functions import einsum, evaluate
 from adcc.Intermediates import Intermediates
 from adcc.AmplitudeVector import AmplitudeVector
 
 
-def mtm_adc0(mp, op, intermediates):
+def mtm_isr0(mp, op, intermediates):
     f1 = op.vo.transpose()
     return AmplitudeVector(ph=f1)
 
 
-def mtm_adc1(mp, op, intermediates):
-    ampl = mtm_adc0(mp, op, intermediates)
+def mtm_isr1(mp, op, intermediates):
+    ampl = mtm_isr0(mp, op, intermediates)
     f1 = - 1.0 * einsum("ijab,jb->ia", mp.t2(b.oovv), op.ov)
     return ampl + AmplitudeVector(ph=f1)
 
 
-def mtm_adc2(mp, op, intermediates):
+def mtm_isr2(mp, op, intermediates):
     t2 = mp.t2(b.oovv)
     p0 = mp.mp2_diffdm
 
-    ampl = mtm_adc1(mp, op, intermediates)
+    ampl = mtm_isr1(mp, op, intermediates)
     f1 = (
         + 0.5 * einsum("ijab,jkbc,ck->ia", t2, t2, op.vo)
         + 0.5 * einsum("ij,aj->ia", p0.oo, op.vo)
@@ -61,14 +61,14 @@ def mtm_adc2(mp, op, intermediates):
     return ampl + AmplitudeVector(ph=f1, pphh=f2)
 
 
-def mtm_cvs_adc0(mp, op, intermediates):
+def mtm_cvs_isr0(mp, op, intermediates):
     f1 = op.vc.transpose()
     return AmplitudeVector(ph=f1)
 
 
-def mtm_cvs_adc2(mp, op, intermediates):
+def mtm_cvs_isr2(mp, op, intermediates):
 
-    ampl = mtm_cvs_adc0(mp, op, intermediates)
+    ampl = mtm_cvs_isr0(mp, op, intermediates)
     f1 = (
         - 0.5 * einsum("bI,ab->Ia", op.vc, intermediates.cvs_p0.vv)
         - 1.0 * einsum("jI,ja->Ia", op.oc, intermediates.cvs_p0.ov)
@@ -78,24 +78,25 @@ def mtm_cvs_adc2(mp, op, intermediates):
 
 
 DISPATCH = {
-    "adc0": mtm_adc0,
-    "adc1": mtm_adc1,
-    "adc2": mtm_adc2,
-    "adc2x": mtm_adc2,  # Identical to ADC(2)
-    "cvs-adc0": mtm_cvs_adc0,
-    "cvs-adc1": mtm_cvs_adc0,  # Identical to CVS-ADC(0)
-    "cvs-adc2": mtm_cvs_adc2,
+    "isr0": mtm_isr0,
+    "isr1s": mtm_isr1,  # Identical to ISR(1)
+    "isr1": mtm_isr1,
+    "isr2": mtm_isr2,
+    "cvs-isr0": mtm_cvs_isr0,
+    "cvs-isr1s": mtm_cvs_isr0,  # Identical to CVS-ISR(0)
+    "cvs-isr1": mtm_cvs_isr0,  # Identical to CVS-ISR(0)
+    "cvs-isr2": mtm_cvs_isr2,
 }
 
 
 def modified_transition_moments(method, ground_state, operator=None,
                                 intermediates=None):
     """Compute the modified transition moments (MTM) for the provided
-    ADC method with reference to the passed ground state.
+    ISR method with reference to the passed ground state.
 
     Parameters
     ----------
-    method: adc.Method
+    method: adcc.IsrMethod
         Provide a method at which to compute the MTMs
     ground_state : adcc.LazyMp
         The MP ground state
@@ -109,8 +110,8 @@ def modified_transition_moments(method, ground_state, operator=None,
     -------
     adcc.AmplitudeVector or list of adcc.AmplitudeVector
     """
-    if not isinstance(method, AdcMethod):
-        method = AdcMethod(method)
+    if not isinstance(method, IsrMethod):
+        method = IsrMethod(method)
     if not isinstance(ground_state, LazyMp):
         raise TypeError("ground_state should be a LazyMp object.")
     if intermediates is None:
diff --git a/adcc/adc_pp/state2state_transition_dm.py b/adcc/adc_pp/state2state_transition_dm.py
index b9c3755e..2caf5a17 100644
--- a/adcc/adc_pp/state2state_transition_dm.py
+++ b/adcc/adc_pp/state2state_transition_dm.py
@@ -22,7 +22,7 @@
 ## ---------------------------------------------------------------------
 from adcc import block as b
 from adcc.LazyMp import LazyMp
-from adcc.AdcMethod import AdcMethod
+from adcc.AdcMethod import IsrMethod
 from adcc.functions import einsum
 from adcc.Intermediates import Intermediates
 from adcc.AmplitudeVector import AmplitudeVector
@@ -32,7 +32,7 @@
 from .util import check_doubles_amplitudes, check_singles_amplitudes
 
 
-def s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates):
+def s2s_tdm_isr0(mp, amplitude_l, amplitude_r, intermediates):
     check_singles_amplitudes([b.o, b.v], amplitude_l, amplitude_r)
     ul1 = amplitude_l.ph
     ur1 = amplitude_r.ph
@@ -43,9 +43,24 @@ def s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates):
     return dm
 
 
-def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
+def s2s_tdm_isr1(mp, amplitude_l, amplitude_r, intermediates):
+    dm = s2s_tdm_isr0(mp, amplitude_l, amplitude_r, intermediates)
+
+    try:
+        check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude_l, amplitude_r)
+        ul1, ul2 = amplitude_l.ph, amplitude_l.pphh
+        ur1, ur2 = amplitude_r.ph, amplitude_r.pphh
+        dm.ov += -2.0 * einsum("jb,ijab->ia", ul1, ur2)
+        dm.vo += -2.0 * einsum("ijab,jb->ai", ul2, ur1)
+    except ValueError:
+        pass
+
+    return dm
+
+
+def s2s_tdm_isr2(mp, amplitude_l, amplitude_r, intermediates):
     check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude_l, amplitude_r)
-    dm = s2s_tdm_adc0(mp, amplitude_l, amplitude_r, intermediates)
+    dm = s2s_tdm_isr1(mp, amplitude_l, amplitude_r, intermediates)
 
     ul1, ul2 = amplitude_l.ph, amplitude_l.pphh
     ur1, ur2 = amplitude_r.ph, amplitude_r.pphh
@@ -59,8 +74,8 @@ def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
     rul1 = einsum('ijab,jb->ia', t2, ul1).evaluate()
     rur1 = einsum('ijab,jb->ia', t2, ur1).evaluate()
 
-    dm.oo = (
-        p1_oo - 2.0 * einsum('ikab,jkab->ij', ur2, ul2)
+    dm.oo += (
+        - 2.0 * einsum('ikab,jkab->ij', ur2, ul2)
         + 0.5 * einsum('ik,kj->ij', p1_oo, p0.oo)
         + 0.5 * einsum('ik,kj->ij', p0.oo, p1_oo)
         - 0.5 * einsum('ikcd,lk,jlcd->ij', t2, p1_oo, t2)
@@ -69,8 +84,8 @@ def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
         - 0.5 * einsum('ikac,kc,ja->ij', t2, rur1, ul1)
         - 1.0 * einsum('ia,ja->ij', rul1, rur1)
     )
-    dm.vv = (
-        p1_vv + 2.0 * einsum('ijac,ijbc->ab', ul2, ur2)
+    dm.vv += (
+        + 2.0 * einsum('ijac,ijbc->ab', ul2, ur2)
         - 0.5 * einsum("ac,cb->ab", p1_vv, p0.vv)
         - 0.5 * einsum("ac,cb->ab", p0.vv, p1_vv)
         - 0.5 * einsum("klbc,klad,cd->ab", t2, t2, p1_vv)
@@ -80,19 +95,18 @@ def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
         + 1.0 * einsum("ia,ib->ab", rur1, rul1)
     )
 
+    # (TODO Move to intermediates)
     p1_ov = -2.0 * einsum("jb,ijab->ia", ul1, ur2).evaluate()
     p1_vo = -2.0 * einsum("ijab,jb->ai", ul2, ur1).evaluate()
 
-    dm.ov = (
-        p1_ov
+    dm.ov += (
         - einsum("ijab,bj->ia", t2, p1_vo)
         - einsum("ib,ba->ia", p0.ov, p1_vv)
         + einsum("ij,ja->ia", p1_oo, p0.ov)
         - einsum("ib,klca,klcb->ia", ur1, t2, ul2)
         - einsum("ikcd,jkcd,ja->ia", t2, ul2, ur1)
     )
-    dm.vo = (
-        p1_vo
+    dm.vo += (
         - einsum("ijab,jb->ai", t2, p1_ov)
         - einsum("ib,ab->ai", p0.ov, p1_vv)
         + einsum("ji,ja->ai", p1_oo, p0.ov)
@@ -103,10 +117,10 @@ def s2s_tdm_adc2(mp, amplitude_l, amplitude_r, intermediates):
 
 
 # Ref: https://doi.org/10.1080/00268976.2013.859313
-DISPATCH = {"adc0": s2s_tdm_adc0,
-            "adc1": s2s_tdm_adc0,       # same as ADC(0)
-            "adc2": s2s_tdm_adc2,
-            "adc2x": s2s_tdm_adc2,      # same as ADC(2)
+DISPATCH = {"isr0": s2s_tdm_isr0,
+            "isr1s": s2s_tdm_isr0,  # Identical to ISR(0)
+            "isr1": s2s_tdm_isr1,
+            "isr2": s2s_tdm_isr2,
             }
 
 
@@ -118,8 +132,8 @@ def state2state_transition_dm(method, ground_state, amplitude_from,
 
     Parameters
     ----------
-    method : str, AdcMethod
-        The method to use for the computation (e.g. "adc2")
+    method : str, IsrMethod
+        The method to use for the computation (e.g. "isr2")
     ground_state : LazyMp
         The ground state upon which the excitation was based
     amplitude_from : AmplitudeVector
@@ -129,8 +143,8 @@ def state2state_transition_dm(method, ground_state, amplitude_from,
     intermediates : adcc.Intermediates
         Intermediates from the ADC calculation to reuse
     """
-    if not isinstance(method, AdcMethod):
-        method = AdcMethod(method)
+    if not isinstance(method, IsrMethod):
+        method = IsrMethod(method)
     if not isinstance(ground_state, LazyMp):
         raise TypeError("ground_state should be a LazyMp object.")
     if not isinstance(amplitude_from, AmplitudeVector):
diff --git a/adcc/adc_pp/state_diffdm.py b/adcc/adc_pp/state_diffdm.py
index 83876f4e..e0f738cf 100644
--- a/adcc/adc_pp/state_diffdm.py
+++ b/adcc/adc_pp/state_diffdm.py
@@ -24,7 +24,7 @@
 
 from adcc import block as b
 from adcc.LazyMp import LazyMp
-from adcc.AdcMethod import AdcMethod
+from adcc.AdcMethod import IsrMethod
 from adcc.functions import einsum
 from adcc.Intermediates import Intermediates
 from adcc.AmplitudeVector import AmplitudeVector
@@ -34,7 +34,7 @@
 from .util import check_doubles_amplitudes, check_singles_amplitudes
 
 
-def diffdm_adc0(mp, amplitude, intermediates):
+def diffdm_isr0(mp, amplitude, intermediates):
     # C is either c(ore) or o(ccupied)
     C = b.c if mp.has_core_occupied_space else b.o
     check_singles_amplitudes([C, b.v], amplitude)
@@ -46,22 +46,36 @@ def diffdm_adc0(mp, amplitude, intermediates):
     return dm
 
 
-def diffdm_adc2(mp, amplitude, intermediates):
-    dm = diffdm_adc0(mp, amplitude, intermediates)  # Get ADC(1) result
+def diffdm_isr1(mp, amplitude, intermediates):
+    dm = diffdm_isr0(mp, amplitude, intermediates)  # Get ISR(0) result
+
+    try:
+        # ISR(1)-d
+        check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
+        u1, u2 = amplitude.ph, amplitude.pphh
+        dm.ov += -2 * einsum("jb,ijab->ia", u1, u2)
+    except ValueError:
+        # no doubles contribution
+        pass
+    return dm
+
+
+def diffdm_isr2(mp, amplitude, intermediates):
+    dm = diffdm_isr1(mp, amplitude, intermediates)  # Get ISR(1) result
     check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
     u1, u2 = amplitude.ph, amplitude.pphh
 
     t2 = mp.t2(b.oovv)
     p0 = mp.mp2_diffdm
-    p1_oo = dm.oo.evaluate()  # ADC(1) diffdm
-    p1_vv = dm.vv.evaluate()  # ADC(1) diffdm
+    p1_oo = dm.oo.evaluate()  # ISR(1) diffdm
+    p1_vv = dm.vv.evaluate()  # ISR(1) diffdm
 
     # Zeroth order doubles contributions
     p2_oo = -einsum("ikab,jkab->ij", u2, u2)
     p2_vv = einsum("ijac,ijbc->ab", u2, u2)
     p2_ov = -2 * einsum("jb,ijab->ia", u1, u2).evaluate()
 
-    # ADC(2) ISR intermediate (TODO Move to intermediates)
+    # ISR(2) intermediate (TODO Move to intermediates)
     ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()
 
     # Compute second-order contributions to the density matrix
@@ -85,8 +99,7 @@ def diffdm_adc2(mp, amplitude, intermediates):
         ).symmetrise()
     )
 
-    dm.ov = (  # adc2_p_ov
-        + p2_ov
+    dm.ov += (  # adc2_p_ov
         - einsum("ijab,jb->ia", t2, p2_ov)
         - einsum("ib,ba->ia", p0.ov, p1_vv)
         + einsum("ij,ja->ia", p1_oo, p0.ov)
@@ -96,52 +109,66 @@ def diffdm_adc2(mp, amplitude, intermediates):
     return dm
 
 
-def diffdm_cvs_adc2(mp, amplitude, intermediates):
-    dm = diffdm_adc0(mp, amplitude, intermediates)  # Get ADC(1) result
+def diffdm_cvs_isr1(mp, amplitude, intermediates):
+    dm = diffdm_isr0(mp, amplitude, intermediates)  # Get ISR(0) result
+
+    try:
+        # ISR(1)-d
+        check_doubles_amplitudes([b.o, b.c, b.v, b.v], amplitude)
+        u1, u2 = amplitude.ph, amplitude.pphh
+        p2_ov = -sqrt(2) * einsum("jb,ijab->ia", u1, u2)
+        dm.ov += p2_ov
+    except ValueError:
+        # no doubles contribution
+        pass
+    return dm
+
+
+def diffdm_cvs_isr2(mp, amplitude, intermediates):
+    dm = diffdm_cvs_isr1(mp, amplitude, intermediates)  # Get cvs-ISR(1) result
     check_doubles_amplitudes([b.o, b.c, b.v, b.v], amplitude)
     u1, u2 = amplitude.ph, amplitude.pphh
 
     t2 = mp.t2(b.oovv)
     p0 = intermediates.cvs_p0
-    p1_vv = dm.vv.evaluate()  # ADC(1) diffdm
+    p1_vv = dm.vv.evaluate()  # ISR(1) diffdm
 
     # Zeroth order doubles contributions
-    p2_ov = -sqrt(2) * einsum("jb,ijab->ia", u1, u2)
     p2_vo = -sqrt(2) * einsum("ijab,jb->ai", u2, u1)
     p2_oo = -einsum("ljab,kjab->kl", u2, u2)
     p2_vv = 2 * einsum("ijac,ijbc->ab", u2, u2)
 
     # Second order contributions
-    # cvs_adc2_dp_oo
-    dm.oo = p2_oo + einsum("ab,ikac,jkbc->ij", p1_vv, t2, t2)
+    # cvs_isr2_dp_oo
+    dm.oo += p2_oo + einsum("ab,ikac,jkbc->ij", p1_vv, t2, t2)
 
-    dm.ov = p2_ov + (  # cvs_adc2_dp_ov
+    dm.ov += (  # cvs_isr2_dp_ov
         - einsum("ka,ab->kb", p0.ov, p1_vv)
         - einsum("lkdb,dl->kb", t2, p2_vo)
         + 1 / sqrt(2) * einsum("ib,klad,liad->kb", u1, t2, u2)
     )
 
-    dm.vv = p1_vv + p2_vv - 0.5 * (  # cvs_adc2_dp_vv
+    dm.vv += p2_vv - 0.5 * (  # cvs_isr2_dp_vv
         + einsum("cb,ac->ab", p1_vv, p0.vv)
         + einsum("cb,ac->ab", p0.vv, p1_vv)
         + einsum("ijbc,ijad,cd->ab", t2, t2, p1_vv)
     )
 
-    # Add 2nd order correction to CVS-ADC(1) diffdm
+    # Add 2nd order correction to CVS-ISR(1) diffdm
     dm.cc -= einsum("kIab,kJab->IJ", u2, u2)
     return dm
 
 
 # dict controlling the dispatch of the state_diffdm function
 DISPATCH = {
-    "adc0": diffdm_adc0,
-    "adc1": diffdm_adc0,       # same as ADC(0)
-    "adc2": diffdm_adc2,
-    "adc2x": diffdm_adc2,
-    "cvs-adc0": diffdm_adc0,
-    "cvs-adc1": diffdm_adc0,   # same as ADC(0)
-    "cvs-adc2": diffdm_cvs_adc2,
-    "cvs-adc2x": diffdm_cvs_adc2,
+    "isr0": diffdm_isr0,
+    "isr1s": diffdm_isr0,   # Identical to ISR(0)
+    "isr1": diffdm_isr1,
+    "isr2": diffdm_isr2,
+    "cvs-isr0": diffdm_isr0,
+    "cvs-isr1s": diffdm_isr0,   # Identical to ISR(0)
+    "cvs-isr1": diffdm_cvs_isr1,
+    "cvs-isr2": diffdm_cvs_isr2,
 }
 
 
@@ -152,8 +179,8 @@ def state_diffdm(method, ground_state, amplitude, intermediates=None):
 
     Parameters
     ----------
-    method : str, AdcMethod
-        The method to use for the computation (e.g. "adc2")
+    method : str, IsrMethod
+        The method to use for the computation (e.g. "isr2")
     ground_state : LazyMp
         The ground state upon which the excitation was based
     amplitude : AmplitudeVector
@@ -161,8 +188,8 @@ def state_diffdm(method, ground_state, amplitude, intermediates=None):
     intermediates : adcc.Intermediates
         Intermediates from the ADC calculation to reuse
     """
-    if not isinstance(method, AdcMethod):
-        method = AdcMethod(method)
+    if not isinstance(method, IsrMethod):
+        method = IsrMethod(method)
     if not isinstance(ground_state, LazyMp):
         raise TypeError("ground_state should be a LazyMp object.")
     if not isinstance(amplitude, AmplitudeVector):
diff --git a/adcc/adc_pp/state_diffdm_2p.py b/adcc/adc_pp/state_diffdm_2p.py
index b6c31f5c..a6e41057 100644
--- a/adcc/adc_pp/state_diffdm_2p.py
+++ b/adcc/adc_pp/state_diffdm_2p.py
@@ -22,7 +22,7 @@
 ## ---------------------------------------------------------------------
 from adcc import block as b
 from adcc.LazyMp import LazyMp
-from adcc.AdcMethod import AdcMethod
+from adcc.AdcMethod import IsrMethod
 from adcc.functions import einsum, zeros_like
 from adcc.Intermediates import Intermediates
 from adcc.AmplitudeVector import AmplitudeVector
@@ -32,7 +32,7 @@
 from .util import check_doubles_amplitudes, check_singles_amplitudes
 
 
-def diffdm_adc0_2p(mp, amplitude, intermediates):
+def diffdm_isr0_2p(mp, amplitude, intermediates):
     check_singles_amplitudes([b.o, b.v], amplitude)
     u1 = amplitude.ph
 
@@ -42,7 +42,7 @@ def diffdm_adc0_2p(mp, amplitude, intermediates):
 
     dm = TwoParticleDensity(mp, symmetry=OperatorSymmetry.HERMITIAN)
 
-    # ADC(1) diffdm
+    # one-particle ISR(0) diffdm
     p1_oo = -einsum("ia,la->il", u1, u1).evaluate()
     p1_vv = einsum("ka,kb->ab", u1, u1).evaluate()
 
@@ -59,8 +59,8 @@ def diffdm_adc0_2p(mp, amplitude, intermediates):
     return dm
 
 
-def diffdm_adc1_2p(mp, amplitude, intermediates):
-    dm = diffdm_adc0_2p(mp, amplitude, intermediates)  # Get ADC(0) result
+def diffdm_isr1s_2p(mp, amplitude, intermediates):
+    dm = diffdm_isr0_2p(mp, amplitude, intermediates)  # Get ISR(0) result
     u1 = amplitude.ph
 
     hf = mp.reference_state
@@ -69,10 +69,10 @@ def diffdm_adc1_2p(mp, amplitude, intermediates):
 
     t2 = mp.t2(b.oovv)
     # TODO move to intermediates!
-    # ADC(1) diffdm
+    # one-particle ISR(0) diffdm
     p1_oo = -einsum("ia,la->il", u1, u1).evaluate()
 
-    # ADC(2) ISR intermediate (TODO Move to intermediates)
+    # ISR(2) ISR intermediate (TODO Move to intermediates)
     ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()
     # new ones
     ru1_ooov = einsum("kc,ijac->ijka", u1, t2).evaluate()
@@ -91,11 +91,65 @@ def diffdm_adc1_2p(mp, amplitude, intermediates):
             einsum("jk,ikab->ijab", p1_oo, t2)
         ).antisymmetrise(0, 1)
     )
+
     return dm
 
 
-def diffdm_adc2_2p(mp, amplitude, intermediates):
-    dm = diffdm_adc1_2p(mp, amplitude, intermediates)  # Get ADC(1) result
+def diffdm_isr1_2p(mp, amplitude, intermediates):
+    dm = diffdm_isr0_2p(mp, amplitude, intermediates)  # Get ISR(0) result
+    u1 = amplitude.ph
+
+    hf = mp.reference_state
+    d_oo = zeros_like(hf.foo)
+    d_oo.set_mask("ii", 1)
+
+    t2 = mp.t2(b.oovv)
+    # TODO move to intermediates!
+    # one-particle ISR(0) diffdm
+    p1_oo = -einsum("ia,la->il", u1, u1).evaluate()
+
+    # ISR(2) ISR intermediate (TODO Move to intermediates)
+    ru1 = einsum("ijab,jb->ia", t2, u1).evaluate()
+    # new ones
+    ru1_ooov = einsum("kc,ijac->ijka", u1, t2).evaluate()
+
+    dm.oovv += (
+        # N^4: O^2V^2 / N^4: O^2V^2
+        + 4.0 * (
+            einsum("ib,ja->ijab", ru1, u1)
+        ).antisymmetrise(0, 1).antisymmetrise(2, 3)
+        # N^5: O^3V^2 / N^4: O^2V^2
+        + 2.0 * (
+            einsum("ijka,kb->ijab", ru1_ooov, u1)
+        ).antisymmetrise(2, 3)
+        # N^5: O^3V^2 / N^4: O^2V^2
+        - 2.0 * (
+            einsum("jk,ikab->ijab", p1_oo, t2)
+        ).antisymmetrise(0, 1)
+    )
+
+    try:
+        # ISR(1)-d
+        check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
+        u2 = amplitude.pphh
+
+        dm.ooov += (
+            # N^5: O^3V^2 / N^4: O^2V^2
+            - 2.0 * einsum("kb,ijab->ijka", u1, u2)
+        )
+
+        dm.ovvv += (
+            # N^5: O^2V^3 / N^4: O^1V^3
+            - 2.0 * einsum("ja,ijbc->iabc", u1, u2)
+        )
+    except ValueError:
+        # no doubles contribution
+        pass
+    return dm
+
+
+def diffdm_isr2_2p(mp, amplitude, intermediates):
+    dm = diffdm_isr1_2p(mp, amplitude, intermediates)  # Get ISR(1) result
     check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
     u1, u2 = amplitude.ph, amplitude.pphh
     hf = mp.reference_state
@@ -160,8 +214,6 @@ def diffdm_adc2_2p(mp, amplitude, intermediates):
         ).antisymmetrise(0, 1).antisymmetrise(2, 3).symmetrise([(0, 2), (1, 3)])
     )
     dm.ooov += (
-        # N^5: O^3V^2 / N^4: O^2V^2
-        - 2.0 * einsum("kb,ijab->ijka", u1, u2)
         # N^6: O^4V^2 / N^4: O^2V^2
         + 1.0 * einsum("ijkl,la->ijka", einsum("klbc,ijbc->ijkl", u2, t2), u1)
         # N^5: O^3V^2 / N^4: O^2V^2
@@ -293,8 +345,6 @@ def diffdm_adc2_2p(mp, amplitude, intermediates):
         ).symmetrise([(0, 2), (1, 3)])
     )
     dm.ovvv += (
-        # N^5: O^2V^3 / N^4: O^1V^3
-        - 2.0 * einsum("ja,ijbc->iabc", u1, u2)
         # N^6: O^3V^3 / N^4: O^1V^3
         + 1.0 * einsum("ijka,jkbc->iabc", einsum("id,jkad->ijka", u1, u2), t2)
         # N^5: O^2V^3 / N^4: O^1V^3
@@ -340,10 +390,10 @@ def diffdm_adc2_2p(mp, amplitude, intermediates):
 
 # dict controlling the dispatch of the state_diffdm function
 DISPATCH = {
-    "adc0": diffdm_adc0_2p,
-    "adc1": diffdm_adc1_2p,
-    "adc2": diffdm_adc2_2p,
-    "adc2x": diffdm_adc2_2p,      # same as ADC(2)
+    "isr0": diffdm_isr0_2p,
+    "isr1s": diffdm_isr1s_2p,
+    "isr1": diffdm_isr1_2p,
+    "isr2": diffdm_isr2_2p,
 }
 
 
@@ -354,8 +404,8 @@ def state_diffdm_2p(method, ground_state, amplitude, intermediates=None):
 
     Parameters
     ----------
-    method : str, AdcMethod
-        The method to use for the computation (e.g. "adc2")
+    method : str, IsrMethod
+        The method to use for the computation (e.g. "isr2")
     ground_state : LazyMp
         The ground state upon which the excitation was based
     amplitude : AmplitudeVector
@@ -363,8 +413,8 @@ def state_diffdm_2p(method, ground_state, amplitude, intermediates=None):
     intermediates : adcc.Intermediates
         Intermediates from the ADC calculation to reuse
     """
-    if not isinstance(method, AdcMethod):
-        method = AdcMethod(method)
+    if not isinstance(method, IsrMethod):
+        method = IsrMethod(method)
     if not isinstance(ground_state, LazyMp):
         raise TypeError("ground_state should be a LazyMp object.")
     if not isinstance(amplitude, AmplitudeVector):
diff --git a/adcc/adc_pp/transition_dm.py b/adcc/adc_pp/transition_dm.py
index 65a40d0f..67a9629d 100644
--- a/adcc/adc_pp/transition_dm.py
+++ b/adcc/adc_pp/transition_dm.py
@@ -24,7 +24,7 @@
 
 from adcc import block as b
 from adcc.LazyMp import LazyMp
-from adcc.AdcMethod import AdcMethod
+from adcc.AdcMethod import IsrMethod
 from adcc.functions import einsum
 from adcc.Intermediates import Intermediates
 from adcc.AmplitudeVector import AmplitudeVector
@@ -34,67 +34,65 @@
 from .util import check_doubles_amplitudes, check_singles_amplitudes
 
 
-def tdm_adc0(mp, amplitude, intermediates):
+def tdm_isr0(mp, amplitude, intermediates):
     # C is either c(ore) or o(ccupied)
     C = b.c if mp.has_core_occupied_space else b.o
     check_singles_amplitudes([C, b.v], amplitude)
     u1 = amplitude.ph
 
-    # Transition density matrix for (CVS-)ADC(0)
+    # Transition density matrix for (CVS-)ISR(0)
     dm = OneParticleDensity(mp, symmetry=OperatorSymmetry.NOSYMMETRY)
     dm[b.v + C] = u1.transpose()
     return dm
 
 
-def tdm_adc1(mp, amplitude, intermediates):
-    dm = tdm_adc0(mp, amplitude, intermediates)  # Get ADC(0) result
-    # adc1_dp0_ov
+def tdm_isr1(mp, amplitude, intermediates):
+    dm = tdm_isr0(mp, amplitude, intermediates)  # Get ISR(0) result
+    # isr1_dp0_ov
     dm.ov = -einsum("ijab,jb->ia", mp.t2(b.oovv), amplitude.ph)
     return dm
 
 
-def tdm_cvs_adc2(mp, amplitude, intermediates):
-    # Get CVS-ADC(1) result (same as CVS-ADC(0))
-    dm = tdm_adc0(mp, amplitude, intermediates)
+def tdm_cvs_isr2(mp, amplitude, intermediates):
+    # Get CVS-ISR(1) result (same as CVS-ISR(0))
+    dm = tdm_isr0(mp, amplitude, intermediates)
     check_doubles_amplitudes([b.o, b.c, b.v, b.v], amplitude)
-    u1 = amplitude.ph
-    u2 = amplitude.pphh
+    u1, u2 = amplitude.ph, amplitude.pphh
 
     t2 = mp.t2(b.oovv)
     p0 = intermediates.cvs_p0
 
-    # Compute CVS-ADC(2) tdm
-    dm.oc = (  # cvs_adc2_dp0_oc
+    # Compute CVS-ISR(2) tdm
+    dm.oc = (  # cvs_isr2_dp0_oc
         - einsum("ja,Ia->jI", p0.ov, u1)
         + (1 / sqrt(2)) * einsum("kIab,jkab->jI", u2, t2)
     )
 
-    # cvs_adc2_dp0_vc
+    # cvs_isr2_dp0_vc
     dm.vc -= 0.5 * einsum("ab,Ib->aI", p0.vv, u1)
     return dm
 
 
-def tdm_adc2(mp, amplitude, intermediates):
-    dm = tdm_adc1(mp, amplitude, intermediates)  # Get ADC(1) result
+def tdm_isr2(mp, amplitude, intermediates):
+    dm = tdm_isr1(mp, amplitude, intermediates)  # Get ISR(1) result
     check_doubles_amplitudes([b.o, b.o, b.v, b.v], amplitude)
-    u1 = amplitude.ph
-    u2 = amplitude.pphh
+    u1, u2 = amplitude.ph, amplitude.pphh
 
     t2 = mp.t2(b.oovv)
     td2 = mp.td2(b.oovv)
     p0 = mp.mp2_diffdm
 
-    # Compute ADC(2) tdm
-    dm.oo = (  # adc2_dp0_oo
+    # Compute ISR(2) tdm
+    dm.oo = (  # isr2_dp0_oo
         - einsum("ia,ja->ij", p0.ov, u1)
         - einsum("ikab,jkab->ji", u2, t2)
     )
-    dm.vv = (  # adc2_dp0_vv
+    dm.vv = (  # isr2_dp0_vv
         + einsum("ia,ib->ab", u1, p0.ov)
         + einsum("ijac,ijbc->ab", u2, t2)
     )
-    dm.ov -= einsum("ijab,jb->ia", td2, u1)  # adc2_dp0_ov
-    dm.vo += 0.5 * (  # adc2_dp0_vo
+    dm.ov -= einsum("ijab,jb->ia", td2, u1)  # isr2_dp0_ov
+    dm.vo += 0.5 * (  # isr2_dp0_vo
         + einsum("ijab,jkbc,kc->ai", t2, t2, u1)
         - einsum("ab,ib->ai", p0.vv, u1)
         + einsum("ja,ij->ai", u1, p0.oo)
@@ -103,14 +101,13 @@ def tdm_adc2(mp, amplitude, intermediates):
 
 
 DISPATCH = {
-    "adc0": tdm_adc0,
-    "adc1": tdm_adc1,
-    "adc2": tdm_adc2,
-    "adc2x": tdm_adc2,
-    "cvs-adc0": tdm_adc0,
-    "cvs-adc1": tdm_adc0,  # No extra contribs for CVS-ADC(1)
-    "cvs-adc2": tdm_cvs_adc2,
-    "cvs-adc2x": tdm_cvs_adc2,
+    "isr0": tdm_isr0,
+    "isr1s": tdm_isr1,  # Identical to ISR(1)
+    "isr1": tdm_isr1,
+    "isr2": tdm_isr2,
+    "cvs-isr0": tdm_isr0,
+    "cvs-isr1": tdm_isr0,  # No extra contribs for CVS-ISR(1)
+    "cvs-isr2": tdm_cvs_isr2,
 }
 
 
@@ -121,8 +118,8 @@ def transition_dm(method, ground_state, amplitude, intermediates=None):
 
     Parameters
     ----------
-    method : str, AdcMethod
-        The method to use for the computation (e.g. "adc2")
+    method : str, IsrMethod
+        The method to use for the computation (e.g. "isr2")
     ground_state : LazyMp
         The ground state upon which the excitation was based
     amplitude : AmplitudeVector
@@ -130,8 +127,8 @@ def transition_dm(method, ground_state, amplitude, intermediates=None):
     intermediates : adcc.Intermediates
         Intermediates from the ADC calculation to reuse
     """
-    if not isinstance(method, AdcMethod):
-        method = AdcMethod(method)
+    if not isinstance(method, IsrMethod):
+        method = IsrMethod(method)
     if not isinstance(ground_state, LazyMp):
         raise TypeError("ground_state should be a LazyMp object.")
     if not isinstance(amplitude, AmplitudeVector):
diff --git a/adcc/tests/AdcMatrix_dense_export_test.py b/adcc/tests/AdcMatrix_dense_export_test.py
index b449ce26..dadf150b 100644
--- a/adcc/tests/AdcMatrix_dense_export_test.py
+++ b/adcc/tests/AdcMatrix_dense_export_test.py
@@ -48,13 +48,13 @@ def test_h2o(self, case: str, method: str):
             method = f"cvs-{method}"
         method: adcc.AdcMethod = adcc.AdcMethod(method)
         n_states = 7
-        if method.level == 1:  # only few states available
+        if method.level.to_int() == 1:  # only few states available
             n_states = 5
             if method.is_core_valence_separated:
                 n_states = 1 if "fv" in case else 2
             elif "fc" in case and "fv" in case:
                 n_states = 4
-        elif method.level == 2 and method.is_core_valence_separated:
+        elif method.level.to_int() == 2 and method.is_core_valence_separated:
             # there seems to be a difference for higher cvs-adc2 states...
             n_states = 3
 
diff --git a/adcc/tests/AdcMatrix_test.py b/adcc/tests/AdcMatrix_test.py
index 325eb2f3..c77d9c6c 100644
--- a/adcc/tests/AdcMatrix_test.py
+++ b/adcc/tests/AdcMatrix_test.py
@@ -226,7 +226,9 @@ def test_properties(self, system: str, case: str, method: str):
         # check that the blocks are correct
         blocks = matrix.axis_blocks
         if matrix.method.adc_type == "pp":
-            assert blocks == ["ph", "pphh", "ppphhh"][:matrix.method.level // 2 + 1]
+            assert blocks == (
+                ["ph", "pphh", "ppphhh"][:matrix.method.level.to_int() // 2 + 1]
+            )
         else:
             raise NotImplementedError(f"Unknown adc type {matrix.method.adc_type}.")
         assert sorted(matrix.axis_spaces.keys(), key=len) == blocks
diff --git a/adcc/tests/AdcMethod_test.py b/adcc/tests/AdcMethod_test.py
new file mode 100644
index 00000000..87a2f97d
--- /dev/null
+++ b/adcc/tests/AdcMethod_test.py
@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2026 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+import pytest
+
+import adcc
+
+adc_methods = [("adc1", None), ("adc2x", None), ("cvs-adc3", None),
+               ("adc", ValueError), ("cvs_adc2", ValueError),
+               ("xyz-adc2", ValueError), ("adc5", NotImplementedError),
+               ("isr2", ValueError)]
+
+isr_methods = [("isr1", None), ("cvs-isr2", None),
+               ("adc", ValueError), ("cvs_isr2", ValueError),
+               ("xyz-isr2", ValueError), ("isr5", NotImplementedError),
+               ("adc2", ValueError)]
+
+
+class TestAdcMethod:
+    @pytest.mark.parametrize("method, expected_exception", adc_methods)
+    def test_validate_adcmethod(self, method, expected_exception):
+        if expected_exception:
+            with pytest.raises(expected_exception):
+                adcc.AdcMethod(method)
+        else:
+            adc_method = adcc.AdcMethod(method)
+            assert adc_method.name == method
+
+    def test_adcmethod(self):
+        method = adcc.AdcMethod("adc2")
+        cvs_method = adcc.AdcMethod("cvs-adc2")
+
+        assert method.name == cvs_method.base_method.name
+
+        method_new_level = method.at_level(1)
+        assert method_new_level._method_base_name == method._method_base_name
+        assert method_new_level.level.to_int() == 1
+
+        as_isr_method = method.as_method(adcc.IsrMethod)
+        assert isinstance(as_isr_method, adcc.IsrMethod)
+
+
+class TestIsrMethod:
+    @pytest.mark.parametrize("method, expected_exception", isr_methods)
+    def test_validate_isrmethod(self, method, expected_exception):
+        if expected_exception:
+            with pytest.raises(expected_exception):
+                adcc.IsrMethod(method)
+        else:
+            isr_method = adcc.IsrMethod(method)
+            assert isr_method.name == method
+
+    def test_isrmethod(self):
+        method = adcc.IsrMethod("isr2")
+        cvs_method = adcc.IsrMethod("cvs-isr2")
+
+        assert method.name == cvs_method.base_method.name
+
+        method_new_level = method.at_level(1)
+        assert method_new_level._method_base_name == method._method_base_name
+        assert method_new_level.level.to_int() == 1
+
+        as_adc_method = method.as_method(adcc.AdcMethod)
+        assert isinstance(as_adc_method, adcc.AdcMethod)
diff --git a/adcc/tests/IsrMatrix_test.py b/adcc/tests/IsrMatrix_test.py
index 3a174fa9..f50f23bf 100644
--- a/adcc/tests/IsrMatrix_test.py
+++ b/adcc/tests/IsrMatrix_test.py
@@ -48,7 +48,7 @@ def test_matrix_vector_product(self, system: str, case: str, kind: str,
         state = testdata_cache.adcc_states(
             system=system, method=method, kind=kind, case=case
         )
-
+        method = method.replace("adc", "isr")
         n_ref = len(state.excitation_vector)
         mp = state.ground_state
         if operator_kind == "electric":  # example of a symmetric operator
@@ -96,6 +96,8 @@ def test_matvec(self):
         state = testdata_cache.adcc_states(
             system=system, method=method, kind=kind, case="gen"
         )
+
+        method = "isr2"
         assert len(state.excitation_vector) > 1
         mp = state.ground_state
         dips = state.reference_state.operators.electric_dipole
diff --git a/adcc/tests/adc_pp/modified_transition_moments_test.py b/adcc/tests/adc_pp/modified_transition_moments_test.py
index 4c05f1d3..cb34ed66 100644
--- a/adcc/tests/adc_pp/modified_transition_moments_test.py
+++ b/adcc/tests/adc_pp/modified_transition_moments_test.py
@@ -54,7 +54,7 @@ def test_modified_transition_moments(system: str, case: str, method: str, kind:
     )[kind]
 
     n_ref = len(state.excitation_vector)
-
+    method = method.replace("adc", "isr")
     if op_kind == "electric":
         dips = state.reference_state.operators.electric_dipole
         ref_tdm = ref["transition_dipole_moments"]
diff --git a/adcc/tests/adc_pp/state_diffdm_test.py b/adcc/tests/adc_pp/state_diffdm_test.py
index 24c07cd5..780fca75 100644
--- a/adcc/tests/adc_pp/state_diffdm_test.py
+++ b/adcc/tests/adc_pp/state_diffdm_test.py
@@ -24,7 +24,7 @@
 import pytest
 
 from adcc import block as b
-from adcc.AdcMethod import AdcMethod
+from adcc.AdcMethod import IsrMethod, MethodLevel
 from adcc.functions import evaluate, einsum
 from adcc.OneParticleDensity import OneParticleDensity
 from adcc.NParticleOperator import OperatorSymmetry
@@ -51,17 +51,27 @@ def calculate_adcn_excitation_energy(self, state):
         mp = state.ground_state
         n_states = len(state.excitation_energy)
         excitation_energy = np.zeros((n_states))
-        method = state.method
-        level = method.level
+        if state.method.name == "adc3":
+            # so we don't forget to switch to the actual implementation
+            with pytest.raises(NotImplementedError):
+                method = IsrMethod("isr3")
+
+            # TODO switch to ISR(3) implementation
+            method = IsrMethod("isr2")
+            method.level = MethodLevel(3)
+
+        else:
+            method = state.property_method
+        level = method.level.to_int()
 
         method_order_minus_one = None
         if level - 1 >= 0:
-            method_order_minus_one = AdcMethod("adc" + str(level - 1))
+            method_order_minus_one = IsrMethod("isr" + str(level - 1))
 
         for es in range(n_states):
             evec = state.excitation_vector[es]
             # TODO switch to ISR(3) implemntation
-            if method.level == 3:
+            if method.level.to_int() == 3:
                 # so we don't forget to switch to the actual implementation
                 with pytest.raises(NotImplementedError):
                     state_diffdm(method, mp, evec)
@@ -77,14 +87,6 @@ def calculate_adcn_excitation_energy(self, state):
             if method_order_minus_one is not None:
                 # two particle part
                 dens_2p = state_diffdm_2p(method_order_minus_one, mp, evec)
-                # go for ISR(1)-d for ADC(2)
-                if method_order_minus_one.level == 1:
-                    dens_2p.ooov += (
-                        - 2.0 * einsum("kb,ijab->ijka", evec.ph, evec.pphh)
-                    )
-                    dens_2p.ovvv += (
-                        - 2.0 * einsum("ja,ijbc->iabc", evec.ph, evec.pphh)
-                    )
                 for block in dens_2p.blocks:
                     # compute
                     # 1/4 [(1 - P_pq) (1 - P_rs) 1 / (n_occ - 1) <pi||ri> delta_qs]
diff --git a/adcc/tests/backends/backends_pcm_test.py b/adcc/tests/backends/backends_pcm_test.py
index 6ca2bb6f..85d281b6 100644
--- a/adcc/tests/backends/backends_pcm_test.py
+++ b/adcc/tests/backends/backends_pcm_test.py
@@ -119,7 +119,7 @@ def test_pcm_linear_response(self, system, case, method, backend):
             atol=1e-5
         )
 
-        state_cis = adcc.ExcitedStates(state, property_method="adc0")
+        state_cis = adcc.ExcitedStates(state, property_method="isr0")
         assert_allclose(
             state_cis.oscillator_strength,
             result["lr_osc_strength"], atol=1e-3
diff --git a/adcc/tests/functionality_test.py b/adcc/tests/functionality_test.py
index 699f1769..ba37653c 100644
--- a/adcc/tests/functionality_test.py
+++ b/adcc/tests/functionality_test.py
@@ -78,10 +78,10 @@ def base_test(self, system: testcases.TestCase, case: str, method: str,
         refmp = getattr(testdata_cache, f"{generator}_data")(
             system=system, method="mp", case=case
         )
-        if res.method.level >= 2:
+        if res.method.level.to_int() >= 2:
             assert res.ground_state.energy_correction(2) == \
                 approx(refmp["mp2"]["energy"])
-        if res.method.level >= 3:
+        if res.method.level.to_int() >= 3:
             if not res.method.is_core_valence_separated:
                 # TODO The latter check can be removed once CVS-MP3 energies
                 #      are implemented
@@ -126,7 +126,8 @@ def base_test(self, system: testcases.TestCase, case: str, method: str,
     def test_functionality(self, system: str, case: str, method: str, kind: str,
                            generator: str):
         method: adcc.AdcMethod = adcc.AdcMethod(method)
-        if generator == "adcman" and "cvs" in case and method.level == 0:
+        if generator == "adcman" and "cvs" in case \
+                and method.level.to_int() == 0:
             pytest.skip("CVS-ADC(0) adcman data is not available")
 
         system: testcases.TestCase = testcases.get_by_filename(system).pop()
@@ -135,12 +136,12 @@ def test_functionality(self, system: str, case: str, method: str, kind: str,
         kwargs = {n_states: 3}
         # only few states available for h2o sto3g
         if system.name == "h2o" and system.basis == "sto-3g":
-            if method.level < 2:  # adc0/adc1
+            if method.level.to_int() < 2:  # adc0/adc1
                 if "cvs" in case and "fv" in case:
                     kwargs[n_states] = 1
                 elif "cvs" in case:
                     kwargs[n_states] = 2
-            elif method.level < 4:  # adc2/adc3
+            elif method.level.to_int() < 4:  # adc2/adc3
                 if "cvs" in case and "fv" in case:  # only 5 states available
                     kwargs["n_guesses"] = 3
 
diff --git a/adcc/tests/state_densities_test.py b/adcc/tests/state_densities_test.py
index b2e8954d..7da0119d 100644
--- a/adcc/tests/state_densities_test.py
+++ b/adcc/tests/state_densities_test.py
@@ -54,7 +54,8 @@ class TestStateDensities:
     @pytest.mark.parametrize("system,case,kind", cases)
     def test_state_diffdm(self, system: str, case: str, kind: str, method: str,
                           generator: str):
-        if "cvs" in case and AdcMethod(method).level == 0 and generator == "adcman":
+        if ("cvs" in case and AdcMethod(method).level.to_int() == 0
+                and generator == "adcman"):
             pytest.skip("No CVS-ADC(0) adcman reference data available.")
         refdata = testdata_cache._load_data(
             system=system, method=method, case=case, source=generator
@@ -82,7 +83,8 @@ def test_state_diffdm(self, system: str, case: str, kind: str, method: str,
                              [c for c in cases if c[2] != "triplet"])
     def test_ground_to_excited_tdm(self, system: str, case: str, kind: str,
                                    method: str, generator: str):
-        if "cvs" in case and AdcMethod(method).level == 0 and generator == "adcman":
+        if ("cvs" in case and AdcMethod(method).level.to_int() == 0
+                and generator == "adcman"):
             pytest.skip("No CVS-ADC(0) adcman reference data available.")
         refdata = testdata_cache._load_data(
             system=system, method=method, case=case, source=generator