You might look into the ab initio Frenkel-Davydov exciton model (AIFDEM). Here, you form direct-product basis states and then perform a nonorthogonal CI in that basis. For a dimer, that might mean matrix elements like <S0 S1 | H | S1 S0>, where either molecule #1 is excited (in the ket) or #2 is excited (in the bra). For AIFDEM, you don’t generate these monomer with constrained DFT but rather via standard monomer-based CIS or TDDFT, and you can include S2, S3, … for additional variational flexibility. There is also an option for ionized basis states |CA> (C=cation, A=anion). This is in the manual, recent example here: https://pubs.acs.org/doi/10.1021/acs.jpcc.0c07932