Merge branch 'main' of https://github.com/JohnsonBiophysicsLab/ionerdss

Author: yingyue2030699

ionerdss/gillespie_simulation/adaptive_rates.py

@@ -143,6 +143,7 @@ def adaptive_bimolecular_rate_1D(
             D_tot = np.sum([diffusion_constants[i] for i in reactant_ids])
             sigma = sigmalist[reactionid]
             ka = ratelist[reactionid]
+            if b < sigma: raise ValueError(f'Too many molecules! b - sigma = {b-sigma:.2f} < 0')
             new_ratelist[reactionid] = ( 1/ka + (b-sigma)/3/D_tot )**(-1)
             reverse_reaction_id = reverse_reaction_pairs[reactionid]
             kb = ratelist[reverse_reaction_id]
@@ -157,6 +158,7 @@ def adaptive_bimolecular_rate_1D(
             D_tot = 2 * diffusion_constants[dimerization_id[0]]
             sigma = sigmalist[reactionid]
             ka = ratelist[reactionid]
+            if b < sigma: raise ValueError(f'Too many molecules! b - sigma = {b-sigma:.2f} < 0')
             new_ratelist[reactionid] = ( 1/ka + (b-sigma)/3/D_tot )**(-1)
             reverse_reaction_id = reverse_reaction_pairs[reactionid]
             kb = ratelist[reverse_reaction_id]

ionerdss/gillespie_simulation/simple_gillespie.py

@@ -135,7 +135,7 @@ def run_Gillespie(
         full_update_scheme (bool): controls if update every propensity entry in each iteration. Defualt `full_update_scheme=False`.
         rate_update_rules (callable): How to update kinetic rates. 
             It takes inputs (rate_constants, y_curr, reactant_matrix, volume)
-        avogadro=6.02214e23,
+        avogadro=6.02214e23, set `avogadro=1` if Molar is not used in units.
         seed (int): Set random seed if given. 
         excluded_volume (float): Consider excluded volume in 1D or other reasonable systems.
         excluded_species_id (int): Index of the species that creates excluded volume.
@@ -176,6 +176,8 @@ def run_Gillespie(
         raise ValueError('To update rates, macroscopic rates must be provided.')
     elif (rate_update_rules is not None) and (bool(macroscopic)):
         need_update_rates = True
+    # print('[DEBUG] rate_update_rules:', rate_update_rules)
+    # print('[DEBUG] need_update_rates:', need_update_rates)
 
     # process excluded volume
     if excluded_volume is not None:
@@ -198,15 +200,14 @@ def reduced_volume(y_curr, excluded_volume, volume):
         volume_corrected_rates = rate_constants 
     elif bool(macroscopic) ^ bool(volume_corrected):
         if macroscopic:
+            new_volume = reduced_volume(y_init, excluded_volume, volume)
             if need_update_rates:
                 volume_corrected_rates = update_rates(
-                    rate_constants, y_init, reduced_volume(y_init, excluded_volume, volume), 
-                    rate_update_rules, reactant_matrix, avogadro
+                    rate_constants, y_init, new_volume, rate_update_rules, reactant_matrix, avogadro
                 )
             else:
                 volume_corrected_rates = rate_constants_volume_correction(
-                    rate_constants, reactant_matrix, 
-                    reduced_volume(y_init, excluded_volume, volume)
+                    rate_constants, reactant_matrix, new_volume, avogadro
                 )
         else:
             volume_corrected_rates = rate_constants
@@ -232,22 +233,26 @@ def reduced_volume(y_curr, excluded_volume, volume):
     if seed is not None: np.random.seed(seed) # set random seed if given
 
     while time < max_time:  # Control simulation time scale
+
+        # update volume if needed
+        if excluded_volume is not None:
+            new_volume = reduced_volume(y, excluded_volume, volume)
+            volume_corrected_rates = rate_constants_volume_correction(
+                rate_constants, reactant_matrix, new_volume, avogadro
+            )
+        else:
+            new_volume = volume
 
+        # update rates if needed (adaptive rates)
+        if need_update_rates: volume_corrected_rates = update_rates(
+                rate_constants, y, new_volume, rate_update_rules, reactant_matrix, avogadro
+            )
+        # print('[DEBUG] volume corrected rates:', volume_corrected_rates)
         if full_update_scheme:
-            # update rates
-            if need_update_rates: volume_corrected_rates = update_rates(
-                    rate_constants, y, reduced_volume(y, excluded_volume, volume), 
-                    rate_update_rules, reactant_matrix, avogadro
-                )
             # Calculate propensity
             propensities = calculate_propensity(y, reactant_matrix, volume_corrected_rates)
 
         else:
-            # update rates
-            if need_update_rates: volume_corrected_rates = update_rates(
-                rate_constants, y, reduced_volume(y, excluded_volume, volume), 
-                rate_update_rules, reactant_matrix, avogadro
-            )
             # Calculate propensity
             propensities = calculate_propensity(y, reactant_matrix, volume_corrected_rates,
                                             propensities, is_propensity_update_needed)

proaffinity-gnn/ProAffinity_GNN_inference.py

@@ -795,6 +795,6 @@ def run_proaffinity_inference(
 
 # Example usage
 if __name__ == "__main__": 
-        run_proaffinity_inference("./pdbfiles/1i4d.pdbqt", "HL,F", verbose=True)
+    run_proaffinity_inference("./pdbfiles/1i4d.pdbqt", "D,AB", verbose=True)