fix: Remove the assumption that residuetypes are available

src/kimmdy/kmc.py

@@ -72,13 +72,19 @@ def rf_kmc(
     reaction_probability = []
 
     # 1. Calculate the probability for each reaction
+
     for recipe in recipe_collection.recipes:
         dt = [x[1] - x[0] for x in recipe.timespans]
-        reaction_probability.append(sum(np.multiply(dt, recipe.rates)))
+        t_interval = recipe.timespans[-1][-1] - recipe.timespans[0][0]
+        reaction_probability.append(sum(np.multiply(dt, recipe.rates)) / t_interval)
 
     # 2. Set the total rate to the sum of individual rates
     probability_cumulative = np.cumsum(reaction_probability)
     probability_sum = probability_cumulative[-1]
+    if any(np.isnan(probability_cumulative)):
+        logger.error(
+            f"nan in cumulative probability: {probability_cumulative}. Please check the timesteps and rates for this step!"
+        )
 
     # 3. Generate two independent uniform (0,1) random numbers u1,u2
     u = rng.random(2)
@@ -96,7 +102,6 @@ def rf_kmc(
     # 5. Calculate the time step associated with mu
     time_delta = np.log(1 / u[1]) / probability_sum
     logger.debug(f"Time delta: {time_delta}\nprobability {reaction_probability}")
-
     return KMCResult(
         recipe=recipe,
         reaction_probability=reaction_probability,

src/kimmdy/topology/topology.py

@@ -971,6 +971,7 @@ def update_partial_charges(self, recipe_steps: list[RecipeStep]) -> None:
         # not need to be repeated
         # Warning: resnr not unique for each residue, 'residues' can map to
         # more than one real residue
+
         residues = {}
         for atom in self.atoms.values():
             if residues.get(atom.resnr) is None:
@@ -1025,6 +1026,14 @@ def update_partial_charges(self, recipe_steps: list[RecipeStep]) -> None:
                     atom2 = self.atoms[step.atom_id_2]
                     # check whether bond exists in topology
                     residue = atom1.residue
+
+                    # only attempt update if residuetype exists
+                    if self.ff.residuetypes.get(residue) is None:
+                        logger.warning(
+                            f"Residue `{residue}` not in list of residuetypes. Not changing partial charges for Bind recipe step {step}!"
+                        )
+                        continue
+
                     bondtype = (atom1.atom, atom2.atom)
                     residue_bond_spec = self.ff.residuetypes[residue].bonds.get(
                         bondtype,
@@ -1385,8 +1394,7 @@ def bind_bond(
             if reactive_moleculetype.angles.get(key) is None:
                 reactive_moleculetype.angles[key] = Angle(key[0], key[1], key[2], "1")
 
-        # add proper and improper dihedrals
-        # add proper dihedrals and pairs
+        # add proper dihedrals
         # TODO; for now this assumes function type 9 for all dihedrals
         # and adds all possible dihedrals (and pairs)
         # later we should check the ff if there are multiple
@@ -1396,12 +1404,11 @@ def bind_bond(
         ) + reactive_moleculetype._get_atom_proper_dihedrals(atompair_nrs[1])
         for key in dihedral_keys:
             if reactive_moleculetype.proper_dihedrals.get(key) is None:
-                reactive_moleculetype.proper_dihedrals[key] = MultipleDihedrals(
-                    *key, "9", dihedrals={"": Dihedral(*key, "9")}
-                )
-            pairkey: tuple[str, str] = tuple(str(x) for x in sorted([key[0], key[3]], key=int))  # type: ignore
-            if reactive_moleculetype.pairs.get(pairkey) is None:
-                reactive_moleculetype.pairs[pairkey] = Pair(pairkey[0], pairkey[1], "1")
+                # check for three membered ring
+                if not key[0] == key[-1]:
+                    reactive_moleculetype.proper_dihedrals[key] = MultipleDihedrals(
+                        *key, "9", dihedrals={"": Dihedral(*key, "9")}
+                    )
 
         # improper dihedral
         dihedral_k_v = reactive_moleculetype._get_atom_improper_dihedrals(
@@ -1417,3 +1424,23 @@ def bind_bond(
             reactive_moleculetype.improper_dihedrals[key].dihedrals[value.c2] = (
                 Dihedral(*key, "4", c0=value.c0, c1=value.c1, periodicity=value.c2)
             )
+
+        # add pairs for the most distant excluded bonded interaction
+        scaled_interaction_keys = {
+            "1": [atompair_nrs],
+            "2": angle_keys,
+            "3": dihedral_keys,
+        }
+        try:
+            for key in scaled_interaction_keys[reactive_moleculetype.nrexcl]:
+                # check for three membered ring in propers, not harmful for bonds, angles
+                if not key[0] == key[-1]:
+                    pairkey: tuple[str, str] = tuple(str(x) for x in sorted([key[0], key[-1]], key=int))  # type: ignore
+                    if reactive_moleculetype.pairs.get(pairkey) is None:
+                        reactive_moleculetype.pairs[pairkey] = Pair(
+                            pairkey[0], pairkey[1], "1"
+                        )
+        except KeyError:
+            logger.exception(
+                f"Unexpected nrexcl: {reactive_moleculetype.nrexcl}, should be in [1,2,3]! Skipping pairs generation."
+            )