diff --git a/fanpy/wfn/cc/pccd_ap1rog.py b/fanpy/wfn/cc/pccd_ap1rog.py
index 0f37949a..485f2a4f 100644
--- a/fanpy/wfn/cc/pccd_ap1rog.py
+++ b/fanpy/wfn/cc/pccd_ap1rog.py
@@ -187,12 +187,13 @@ def __init__(
             annihilation to the creation operators.
 
         """
+        # indices are used for assign_exops, without any modification 
+        # we can pass it here to BaseCC as assign_exops gets overwritten by PCCD and its child classes
         super().__init__(
-            nelec, nspin, memory=memory, params=params, exop_combinations=exop_combinations, refresh_exops=refresh_exops
+            nelec, nspin, memory=memory, params=params, exop_combinations=exop_combinations, indices=indices, refresh_exops=refresh_exops
         )
         self.assign_s_type(s_type=s_type)
         self.assign_ranks(ranks=ranks)
-        self.assign_exops(indices=indices)
         self.assign_refwfn(refwfn=refwfn)
 
         # mapping from indices to exops
diff --git a/tests/test_wfn_cc_ap1rog_generalized.py b/tests/test_wfn_cc_ap1rog_generalized.py
new file mode 100644
index 00000000..a7b3b2f2
--- /dev/null
+++ b/tests/test_wfn_cc_ap1rog_generalized.py
@@ -0,0 +1,57 @@
+""" Test for AP1roGSD """
+
+import numpy as np
+import pytest
+from fanpy.tools.sd_list import sd_list
+from fanpy.tools.slater import excite, ground
+
+from fanpy.wfn.cc.ap1rog_generalized import AP1roGSDGeneralized
+
+def test_assign_ranks():
+    # only default rank allowed
+    with pytest.raises(TypeError):
+        AP1roGSDGeneralized(2, 4, ranks=2)
+
+    with pytest.raises(TypeError):
+        AP1roGSDGeneralized(2, 4, ranks=[1, 2, 3])
+
+    # AP1roGSD only has single and double excitations
+    wfn = AP1roGSDGeneralized(2, 4)
+    assert wfn.ranks == [1, 2]
+
+def test_assign_expos():
+    # only default exops allowed
+    with pytest.raises(TypeError):
+        AP1roGSDGeneralized(2, 4, indices=[[0, 2]])
+    # build double excitations for PCCD:
+    n_elec = 6
+    n_spin = 14
+    wfn = AP1roGSDGeneralized(n_elec, n_spin)
+
+    # build single and double excitations manually
+    doubles = sd_list(n_elec, n_spin, exc_orders=[2], spin=0, seniority=0)
+    singles = sd_list(n_elec, n_spin, exc_orders=[1])
+    exc_sds = []
+    # sd list adds ground state as well because it satisfies the spin and sen restrictions
+    # ground is added as the first element, we jump over it here
+    exc_sds.extend(doubles[1:])
+    exc_sds.extend(singles[1:])
+
+    # generate excited sds from exc ops in wfn 
+    wfn_exc_sds = []
+    ground_state = ground(n_elec, n_spin)
+    for exc_orders in wfn.exops.keys():
+        sd = excite(ground_state, *exc_orders)
+        wfn_exc_sds.append(sd)
+
+    # convert arrays to numpy so that comparison easier
+    exc_sds_np = np.asarray(exc_sds, dtype=int)
+    exc_sds_np = np.sort(exc_sds_np)
+    # sort arrays 
+    wfn_exc_sds_np = np.asarray(wfn_exc_sds, dtype=int)
+    wfn_exc_sds_np = np.sort(wfn_exc_sds_np)
+
+    # compare sds
+    np.testing.assert_array_equal(exc_sds_np, wfn_exc_sds_np)
+
+
diff --git a/tests/test_wfn_cc_ap1rog_spin.py b/tests/test_wfn_cc_ap1rog_spin.py
new file mode 100644
index 00000000..c3844b12
--- /dev/null
+++ b/tests/test_wfn_cc_ap1rog_spin.py
@@ -0,0 +1,53 @@
+"""tests for fanpy.wfn.cc.ap1rog_spin"""
+
+import pytest
+import numpy as np
+
+from fanpy.wfn.cc.ap1rog_spin import AP1roGSDSpin
+from fanpy.tools.slater import excite, ground
+from fanpy.tools.sd_list import sd_list
+
+def test_assign_ranks():
+    # only default rank allowed
+    with pytest.raises(TypeError):
+        AP1roGSDSpin(2, 4, ranks=2)
+
+    with pytest.raises(TypeError):
+        AP1roGSDSpin(2, 4, ranks=[1, 2, 3])
+
+    # AP1roGSD only has single and double excitations
+    wfn = AP1roGSDSpin(2, 4)
+    assert wfn.ranks == [1, 2]
+
+def test_assign_expos():
+    # only default exops allowed
+    with pytest.raises(TypeError):
+        AP1roGSDSpin(2, 4, indices=[[0, 2]])
+    # build double excitations for PCCD:
+    n_elec = 6
+    n_spin = 14
+    wfn = AP1roGSDSpin(n_elec, n_spin)
+
+    # build single and double excitations manually
+    doubles = sd_list(n_elec, n_spin, exc_orders=[2], spin=0, seniority=0)
+    singles = sd_list(n_elec, n_spin, exc_orders=[1], spin=0)
+    exc_sds = []
+    exc_sds.extend(doubles[1:])
+    exc_sds.extend(singles[1:])
+
+    # generate excited sds from exc ops in wfn 
+    wfn_exc_sds = []
+    ground_state = ground(n_elec, n_spin)
+    for exc_orders in wfn.exops.keys(): # keys are excitation inds
+        sd = excite(ground_state, *exc_orders)
+        wfn_exc_sds.append(sd)
+
+    # convert arrays to numpy so that comparison easier
+    exc_sds_np = np.asarray(exc_sds, dtype=int)
+    exc_sds_np = np.sort(exc_sds_np)
+    # sort arrays 
+    wfn_exc_sds_np = np.asarray(wfn_exc_sds, dtype=int)
+    wfn_exc_sds_np = np.sort(wfn_exc_sds_np)
+
+    # compare sds
+    np.testing.assert_array_equal(exc_sds_np, wfn_exc_sds_np)
\ No newline at end of file
diff --git a/tests/test_wfn_cc_apset1gro_d.py b/tests/test_wfn_cc_apset1gro_d.py
index 8762b1af..46a5ffa7 100644
--- a/tests/test_wfn_cc_apset1gro_d.py
+++ b/tests/test_wfn_cc_apset1gro_d.py
@@ -1,6 +1,7 @@
 """Test fanpy.wavefunction.cc.apset1rog_d."""
 import pytest
 from fanpy.wfn.cc.apset1rog_d import APset1roGD
+from fanpy.tools.slater import ground, vir_indices
 
 
 class TempAPset1roGD(APset1roGD):
@@ -10,7 +11,7 @@ def __init__(self):
         self.exop_combinations = {}
 
 
-def test_assign_exops():
+def test_assign_exops_default_inds():
     """Test APset1roGD.assign_exops."""
     test = TempAPset1roGD()
     test.assign_nelec(4)
@@ -21,3 +22,52 @@ def test_assign_exops():
     test.assign_exops()
     assert test.exops == {(0, 4, 2, 6): 0, (0, 4, 2, 7): 1, (0, 4, 3, 6): 2, (0, 4, 3, 7): 3,
                           (1, 5, 2, 6): 4, (1, 5, 2, 7): 5, (1, 5, 3, 6): 6, (1, 5, 3, 7): 7}
+    
+def test_assign_exops_non_default_inds():
+    test = TempAPset1roGD()
+    nelec = 6
+    nspin = 14
+    test.assign_nelec(nelec)
+    test.assign_nspin(nspin)
+    test.assign_refwfn()
+    ground_state = ground(nelec, nspin)
+    vir_idx = vir_indices(ground_state, nspin).tolist()
+    alpha_idx = vir_idx[:len(vir_idx)//2]
+    beta_idx = vir_idx[len(vir_idx)//2:]
+    # trigger non-default case, but with same alpha beta subsets
+    # as for the default case. This allows easier testing. 
+    test.assign_exops([alpha_idx, beta_idx])
+
+    test_default = TempAPset1roGD()
+    test_default.assign_nelec(nelec)
+    test_default.assign_nspin(nspin)
+    test_default.assign_refwfn()
+    test_default.assign_exops()
+
+    assert test_default.exops == test.exops
+
+def test_assign_exops_errors():
+    test = TempAPset1roGD()
+    test.assign_nelec(4)
+    test.assign_nspin(8)
+    test.assign_refwfn()
+
+    # not having two sets:
+    with pytest.raises(TypeError, match="`indices` must have exactly 2 elements"):
+        test.assign_exops([[ 2, 3 ], [6], [7, 5]])
+
+    # non int inds: 
+    with pytest.raises(TypeError, match="The elements of `indices` must be lists of non-negative ints"):
+        test.assign_exops([[ 2, 3 ], [ 6, 7.5]])
+
+    # negative inds:
+    with pytest.raises(ValueError, match="All `indices` must be lists of non-negative ints"):
+        test.assign_exops([[ -2, 3 ], [ 6, 7]])
+
+    # set has occ orbitals
+    with pytest.raises(ValueError, match="`indices` cannot correspond to occupied spin-orbitals"):
+        test.assign_exops([[0, 1, 4, 5], [2, 3, 6, 7]])
+
+    # not a disjoint set
+    with pytest.raises(ValueError, match="The sets of creation operators must be disjoint"):
+        test.assign_exops([[ 3, 6], [ 6, 7]])
\ No newline at end of file
diff --git a/tests/test_wfn_cc_apset1gro_sd.py b/tests/test_wfn_cc_apset1gro_sd.py
index 676c2f15..03ed413d 100644
--- a/tests/test_wfn_cc_apset1gro_sd.py
+++ b/tests/test_wfn_cc_apset1gro_sd.py
@@ -87,3 +87,56 @@ def test_generate_possible_exops():
         test.exop_combinations[(0, 1, 4, 2, 3, 6)][2][1],
         [test.get_ind((1, 3)), test.get_ind((0, 4, 2, 6)), sign if sign == 1 else 0],
     )
+
+def test_assign_exops_errors():
+    test = TempAPset1GroSD()
+    test.assign_nelec(4)
+    test.assign_nspin(8)
+    test.assign_refwfn()
+    # wrong type
+    with pytest.raises(TypeError, match="The elements of `indices` must be lists of non-negative ints" ):
+        test.assign_exops(["not a list", "not a list"])
+
+    # not having two sets:
+    with pytest.raises(TypeError, match="`indices` must have exactly 2 elements"):
+        test.assign_exops([[ 2, 3 ], [6], [7, 5]])
+
+    # non int inds: 
+    with pytest.raises(TypeError, match="The elements of `indices` must be lists of non-negative ints"):
+        test.assign_exops([[ 2, 3 ], [ 6, 7.5]])
+
+    # negative inds:
+    with pytest.raises(ValueError, match="All `indices` must be lists of non-negative ints"):
+        test.assign_exops([[ -2, 3 ], [ 6, 7]])
+
+    # set has occ orbitals
+    with pytest.raises(ValueError, match="`indices` cannot correspond to occupied spin-orbitals"):
+        test.assign_exops([[0, 1, 4, 5], [2, 3, 6, 7]])
+
+    # not a disjoint set
+    with pytest.raises(ValueError, match="The sets of annihilation operators must be disjoint"):
+        test.assign_exops([[ 3, 6], [ 6, 7]])
+
+def test_assign_exops_non_default_inds():
+    test = TempAPset1GroSD()
+    nelec = 6
+    nspin = 14
+    test.assign_nelec(nelec)
+    test.assign_nspin(nspin)
+    test.assign_refwfn()
+    ground_state = slater.ground(nelec, nspin)
+    vir_idx = slater.vir_indices(ground_state, nspin).tolist()
+    alpha_idx = vir_idx[:len(vir_idx)//2]
+    beta_idx = vir_idx[len(vir_idx)//2:]
+    # trigger non-default case, but with same alpha beta subsets
+    # as for the default case. This allows easier testing.
+    # print("DEBUG >>> test.exop_combs: ", test.exops)
+    test.assign_exops([alpha_idx, beta_idx])
+
+    test_default = TempAPset1GroSD()
+    test_default.assign_nelec(nelec)
+    test_default.assign_nspin(nspin)
+    test_default.assign_refwfn()
+    test_default.assign_exops()
+
+    assert test_default.exops == test.exops
\ No newline at end of file
diff --git a/tests/test_wfn_cc_pccd_ap1rog.py b/tests/test_wfn_cc_pccd_ap1rog.py
index 66a295b2..4c961c92 100644
--- a/tests/test_wfn_cc_pccd_ap1rog.py
+++ b/tests/test_wfn_cc_pccd_ap1rog.py
@@ -3,6 +3,8 @@
 import numpy as np
 from fanpy.tools import slater
 from fanpy.wfn.cc.pccd_ap1rog import PCCD
+from fanpy.wfn.cc.ap1rog_generalized import AP1roGSDGeneralized
+from fanpy.wfn.cc.apset1rog_sd import APset1roGSD
 
 
 class TempPCCD(PCCD):
@@ -11,6 +13,7 @@ def __init__(self):
         self._cache_fns = {}
         self.exop_combinations = {}
 
+######### ASSIGN METHODS #########
 
 def test_assign_nelec():
     """Test PCCD.assign_nelec."""
@@ -110,6 +113,7 @@ def test_init_s_type():
     test = PCCD(4, 8, s_type="sen-v")
     assert test.s_type == "sen-v"
 
+######### OVERLAP #########
 
 def tests_type_effct_on_pCCD_overlap():
     """Test s_type effect on pCCD overlap.
@@ -120,13 +124,248 @@ def tests_type_effct_on_pCCD_overlap():
     sd = 0b10100011
 
     test_free = PCCD(4, 8, s_type="free")
+    params = np.random.rand(test_free.nparams)
+    test_free.assign_params(params)
     test_seno = PCCD(4, 8, s_type="sen-o")
+    test_seno.assign_params(params)
 
     olp_free = test_free.get_overlap(sd)
     olp_seno = test_seno.get_overlap(sd)
 
     assert olp_free == olp_seno
 
+def test_olp_identity_is_one():
+    """Test that the reference determinant has unit overlap."""
+    test = PCCD(4, 8)
+
+    assert test._olp(test.refwfn) == pytest.approx(1.0)
+
+def test_olp_pair_excitation_matches_parameter():
+    """Test that a single pair excitation returns the matching amplitude."""
+    test = PCCD(4, 8)
+    test.assign_params(np.array([0.25, -0.5, 0.75, -1.25]))
+
+    sd = slater.excite(test.refwfn, 0, 4, 2, 6)
+    sign = slater.sign_excite(test.refwfn, [0, 4], [2, 6])
+
+    assert test._olp(sd) == pytest.approx(sign * test.params[test.get_ind((0, 4, 2, 6))])
+
+def test_olp_broken_pair_is_zero():
+    """Test that seniority-breaking determinants have zero overlap in pCCD."""
+    test = PCCD(4, 8)
+    test.assign_params(np.array([0.25, -0.5, 0.75, -1.25]))
+
+    sd = slater.excite(test.refwfn, 0, 2)
+
+    assert test._olp(sd) == pytest.approx(0.0)
+
+def test_olp_single_excitation_matches_parameter():
+    """Test that a single pair excitation returns the matching amplitude."""
+    test = AP1roGSDGeneralized(4, 8)
+    test.assign_params(np.random.rand(test.nparams))
+
+    sd = slater.excite(test.refwfn, 0, 2)
+    sign = slater.sign_excite(test.refwfn, [0], [2])
+
+    assert test._olp(sd) == pytest.approx(sign * test.params[test.get_ind((0, 2))])
+
+@pytest.mark.parametrize("s_type", ["sen-o", "sen-v", "sen-ov"])
+def test_olp_filtering(monkeypatch, s_type):
+    """Check if single excitations get filtered for s_types o, v, ov
+    Note: this does not check if the filtering logic is working properly
+    """
+
+    wfn_olp = APset1roGSD(
+        nelec=4,
+        nspin=8,
+        s_type=s_type,
+    )
+    wfn_indices_multi = APset1roGSD(
+        nelec=4,
+        nspin=8,
+        s_type=s_type,
+    )
+    #
+    # reference:
+    refsd = slater.ground(4, 8)
+
+    wfn_olp.refwfn = refsd
+    wfn_indices_multi.rewfwfn = refsd
+
+    # Build excited determinant
+    # annihilate:
+    #   pair from spatial orbital 0
+    #   alpha electron from spatial orbital 1
+    #
+    # create:
+    #   pair in spatial orbital 2
+    #   alpha electron in spatial orbital 3
+    #
+
+    excited_sd = slater.excite(refsd, 0, 1, 4, 2, 3, 6)
+
+    # Generate all possible exc combinations
+    a_inds, c_inds = slater.diff_orbs(refsd, excited_sd)
+    wfn_indices_multi.generate_possible_exops(a_inds, c_inds)
+
+    original_indices = wfn_indices_multi.exop_combinations[tuple(a_inds + c_inds)]
+
+    # Capture filtered indices_multi
+
+    captured = {}
+
+    def fake_product_amplitudes_multi(indices_multi):
+        captured["indices_multi"] = indices_multi
+        return 1.0
+
+    monkeypatch.setattr(
+        wfn_olp,
+        "product_amplitudes_multi",
+        fake_product_amplitudes_multi,
+    )
+
+    # Run overlap
+    wfn_olp._olp(excited_sd)
+
+    olp_indices = captured["indices_multi"]
+
+    # all sen types need to remove the case, where the excited SD is generated
+    # from single excitations only
+
+    removed_any = False
+
+    for exc_order in original_indices:
+        if len(olp_indices[exc_order]) < len(original_indices[exc_order]):
+            removed_any = True
+
+    assert removed_any
+
+    # make sure some excitations survived
+
+    surviving = any(
+        len(olp_indices[exc_order]) > 0
+        for exc_order in olp_indices
+    )
+
+    assert surviving
+
+def test_olp_no_filtering_for_free(monkeypatch):
+    """
+    Make sure that seniority type free does not filter single excitations. 
+    """
+
+    wfn_olp = APset1roGSD(
+        nelec=4,
+        nspin=8,
+        s_type="free",
+    )
+    wfn_indices_multi = APset1roGSD(
+        nelec=4,
+        nspin=8,
+        s_type="free",
+    )
+    refsd = slater.create(0, 0, 1, 4, 5)
+
+    wfn_olp.refwfn = refsd
+    wfn_indices_multi.rewfwfn = refsd
+
+    # Build excited determinant
+    # annihilate:
+    #   pair from spatial orbital 0
+    #   alpha electron from spatial orbital 1
+    #
+    # create:
+    #   pair in spatial orbital 2
+    #   alpha electron in spatial orbital 3
+    #
+
+    excited_sd = refsd
+
+    excited_sd = slater.excite(refsd, 0, 1, 4, 2, 3, 6)
+
+    # Generate all excitation combinations
+
+    a_inds, c_inds = slater.diff_orbs(refsd, excited_sd)
+    wfn_indices_multi.generate_possible_exops(a_inds, c_inds)
+
+    original_indices = wfn_indices_multi.exop_combinations[tuple(a_inds + c_inds)]
+
+    # Capture filtered indices_multi
+
+    captured = {}
+
+    def fake_product_amplitudes_multi(indices_multi):
+        captured["indices_multi"] = indices_multi
+        return 1.0
+
+    monkeypatch.setattr(
+        wfn_olp,
+        "product_amplitudes_multi",
+        fake_product_amplitudes_multi,
+    )
+
+    # Run overlap
+
+    wfn_olp._olp(excited_sd)
+
+    olp_indices = captured["indices_multi"]
+
+    removed_any = False
+
+    for exc_order in original_indices:
+        if len(olp_indices[exc_order]) < len(original_indices[exc_order]):
+            removed_any = True
+
+    assert not removed_any # we should not be removing exc with s_type free
+
+    # make sure there are exc operators present
+    # this is a bit redundant, but keeping it here just in case
+
+    surviving = any(
+        len(olp_indices[exc_order]) > 0
+        for exc_order in olp_indices
+    )
+
+    assert surviving
+
+######### DERIV ##########
+
+@pytest.mark.parametrize("sd,s_type", [[0b10100101, "sen-o"], 
+                                       [0b11000101, "sen-v"], 
+                                       [0b10010101, "sen-ov"], 
+                                       [0b11000011, "free"], 
+                                       [0b00110011, "sen-o"]])
+def test_olp_derivative_finite_difference(sd, s_type):
+    """Test overlap Hessian with finite differences."""
+
+    test = AP1roGSDGeneralized(4, 8, s_type = s_type)
+    test.assign_params(np.random.rand(test.nparams))
+
+    h = 1e-7
+
+    analytic = test._olp_deriv(sd)
+
+    numerical = np.zeros_like(analytic)
+
+    orig = test.params.copy()
+
+    for j in range(test.nparams):
+
+        test.params = orig.copy()
+        test.params[j] += h
+        plus = test._olp(sd)
+
+        test.params = orig.copy()
+        test.params[j] -= h
+        minus = test._olp(sd)
+
+        numerical[j] = (plus - minus) / (2 * h)
+
+    test.params = orig
+
+    assert np.allclose(analytic, numerical, atol=1e-5)
+
+######### DOUBLE DERIV #########
 
 def test_olp_double_derivative_shape():
     """Test overlap Hessian shape."""
@@ -142,7 +381,7 @@ def test_olp_double_derivative_symmetric():
     """Test overlap Hessian symmetry."""
 
     test = PCCD(4, 8)
-
+    test.assign_params(np.random.rand(test.nparams))
     sd = 0b11001100
 
     hess = test._olp_double_derivative(sd)
@@ -154,20 +393,25 @@ def test_olp_double_derivative_zero_diagonal():
     """Test Hessian diagonal vanishes."""
 
     test = PCCD(4, 8)
-
+    test.assign_params(np.random.rand(test.nparams))
     sd = 0b11001100
 
     hess = test._olp_double_derivative(sd)
 
     assert np.allclose(np.diag(hess), 0)
 
+# Test double deriv with multiple s types
 
-def test_olp_double_derivative_finite_difference():
+@pytest.mark.parametrize("sd,s_type", [[0b10100101, "sen-o"], 
+                                       [0b11000101, "sen-v"], 
+                                       [0b10010101, "sen-ov"], 
+                                       [0b11000011, "free"], 
+                                       [0b00110011, "sen-o"]])
+def test_olp_double_derivative_finite_difference(sd, s_type):
     """Test overlap Hessian with finite differences."""
 
-    test = PCCD(4, 8)
-
-    sd = 0b11001100
+    test = AP1roGSDGeneralized(4, 8, s_type = s_type)
+    test.assign_params(np.random.rand(test.nparams))
 
     h = 1e-7