Updates to README and RTD installation instructions. Also correction to anchor location change tests.

Author: lvotapka

README.md

@@ -169,6 +169,7 @@ providing feedback, bug reports, or other comments.
 * Rommie Amaro (principal investigator)
 * Ilker Deveci (developer)
 * Hilliary Frank (contributor)
+* Sasha Heyneman (contributor)
 * Ben Jagger (developer)
 * Anand Ojha (developer)
 * Andy Stokely (developer)
@@ -194,7 +195,7 @@ You may also optionally cite the following papers:
 
 ### Copyright
 
-Copyright (c) 2023, Lane Votapka
+Copyright (c) 2024, Lane Votapka
 
 
 #### Acknowledgements

seekr2/modules/common_analyze.py

@@ -515,9 +515,9 @@ def calculate_kinetics(self):
         end_milestones = {}
         bulk_milestones = []
         MFPTs = {}
+        MFPTs_anchors_to_bulk_milestones = defaultdict(float)
         k_off = 0.0
         k_ons = {}
-        bulk_milestone = None
         for anchor in self.model.anchors:
             if anchor.endstate:
                 for milestone_id in anchor.get_ids():
@@ -532,8 +532,8 @@ def calculate_kinetics(self):
                         if anchor.alias_from_id(milestone_id) == 1: 
                             # TODO: hacky
                             continue
-                    #bulk_milestones.append(milestone_id)
-                    bulk_milestone = milestone_id
+                    bulk_milestones.append(milestone_id)
+                    #bulk_milestone = milestone_id
 
         if np.any(self.Q.sum(axis=1) > 1.E-10):
             problem_milestone = np.argmin(self.Q.T.sum(axis=1))
@@ -547,75 +547,92 @@ def calculate_kinetics(self):
                 problem_milestone)
             raise Exception(error_msg)
         B, tau = PyGT.tools.load_CTMC(self.Q.T)
-        assert len(bulk_milestones) <= 1, \
-            "There cannot be more than one bulk milestone."
-        if bulk_milestone is not None:
-            n = self.Q.shape[0]
-            for end_milestone in end_milestones:
-                if end_milestone == bulk_milestone:
+        #assert len(bulk_milestones) <= 1, \
+        #    "There cannot be more than one bulk milestone."
+        #if bulk_milestone is not None:
+        p_i_hat_normalize_per_anchor = defaultdict(float)
+        for end_milestone_src in end_milestones:
+            src_anchor_index = end_milestones[end_milestone_src]
+            p_i_hat_normalize_per_anchor[src_anchor_index] \
+                    += self.p_i[end_milestone_src]
+                
+        n = self.Q.shape[0]
+        for bulk_milestone in bulk_milestones:
+            for end_milestone_src in end_milestones:
+                
+                if end_milestone_src in bulk_milestones:
                     continue
 
                 mfpt_ij, mfpt_ji = PyGT.mfpt.compute_MFPT(
-                    end_milestone, bulk_milestone, B, tau)
+                    end_milestone_src, bulk_milestone, B, tau)
                 mfpt = mfpt_ji
-                MFPTs[(end_milestones[end_milestone], "bulk")] = mfpt
-            
-            MFPT_to_bulk = 0
+                src_anchor_index = end_milestones[end_milestone_src]
+                weighted_mfpt = mfpt * self.p_i[end_milestone_src] \
+                    / p_i_hat_normalize_per_anchor[src_anchor_index]
+                #MFPTs[(src_anchor_index, "bulk")] += weighted_mfpt
+                MFPTs_anchors_to_bulk_milestones[(src_anchor_index, bulk_milestone)] \
+                    += weighted_mfpt
+        
+        for (mfpt_key, mfpt_value) in MFPTs_anchors_to_bulk_milestones.items():
+            (src_anchor_index, bulk_milestone) = mfpt_key
+            if src_anchor_index not in MFPTs:
+                MFPTs[(src_anchor_index, "bulk")] = mfpt_value
+            elif mfpt_value < MFPTs[src_anchor]:
+                MFPTs[(src_anchor_index, "bulk")] = mfpt_value
+                
+        MFPTs_to_bulk_dict = {}
+        for bulk_milestone in bulk_milestones:
             for i in range(self.model.num_milestones):
                 mfpt_ij, mfpt_ji = PyGT.mfpt.compute_MFPT(
                     i, bulk_milestone, B, tau)
                 mfpt = mfpt_ji
+                if i not in MFPTs_to_bulk_dict:
+                    MFPTs_to_bulk_dict[i] = mfpt
+                elif mfpt < MFPTs_to_bulk_dict[i]:
+                    MFPTs_to_bulk_dict[i] = mfpt
+                
+        if len(bulk_milestones) > 0:
+            MFPT_to_bulk = 0
+            for i in range(self.model.num_milestones):
+                mfpt = MFPTs_to_bulk_dict[i]
                 MFPT_to_bulk += mfpt * self.p_i[i]
 
             # convert to 1/s
             k_off = 1.0e12 / MFPT_to_bulk
             self.k_off = k_off
 
-        p_i_hat_normalize = defaultdict(float)
-        
-        for end_milestone_dest in end_milestones:
-            if end_milestone_dest == bulk_milestone:
-                continue
-            for end_milestone_src in end_milestones:
-                if end_milestones[end_milestone_dest] \
-                        == end_milestones[end_milestone_src]:
-                    continue
-                if end_milestone_src == bulk_milestone:
-                    continue
-                p_i_hat_normalize[end_milestone_src] \
-                    += self.p_i[end_milestone_src]
-        
         # Next, compute the MFPTs between different states
         possible_MFPTs = defaultdict(float)
         for end_milestone_dest in end_milestones:
-            if end_milestone_dest == bulk_milestone:
+            dest_anchor_index = end_milestones[end_milestone_dest]
+            if end_milestone_dest in bulk_milestones:
                 continue
 
             for end_milestone_src in end_milestones:
-                if end_milestones[end_milestone_dest] \
-                        == end_milestones[end_milestone_src]:
+                src_anchor_index = end_milestones[end_milestone_src]
+                if dest_anchor_index == src_anchor_index:
                     # don't get the MFPT from a milestone to itself
                     continue
-                if end_milestone_src == bulk_milestone:
+                if end_milestone_src in bulk_milestones:
                     # a bulk milestone will never be a source
                     continue
 
                 mfpt_ij, mfpt_ji = PyGT.mfpt.compute_MFPT(
                     end_milestone_src, end_milestone_dest, B, tau)
                 mfpt = mfpt_ji
-                possible_MFPTs[(end_milestones[end_milestone_src], 
-                       end_milestone_dest)] += mfpt \
+                possible_MFPTs[(src_anchor_index, end_milestone_dest)] += mfpt \
                        * self.p_i[end_milestone_src] \
-                       / p_i_hat_normalize[end_milestone_src]
+                       / p_i_hat_normalize_per_anchor[src_anchor_index]
 
         # TODO: verify if correct: this finds the destination with the 
         # fastest MFPT and just uses that, although the sources are
         # statistically weighed.
         for end_milestone_src in end_milestones:
+            src_anchor_index = end_milestones[end_milestone_src]
             for end_milestone_dest in end_milestones:
-                key = (end_milestones[end_milestone_src], end_milestone_dest)
-                full_key = (end_milestones[end_milestone_src], 
-                            end_milestones[end_milestone_dest])
+                dest_anchor_index = end_milestones[end_milestone_dest]
+                key = (src_anchor_index, end_milestone_dest)
+                full_key = (src_anchor_index, dest_anchor_index)
                 if key in possible_MFPTs:
                     if full_key in MFPTs:
                         # Get the minimum MFPT

seekr2/modules/common_sim_browndye2.py

@@ -88,7 +88,7 @@ def __init__(self):
         self.dielectric = 78.0
         self.relative_viscosity = 1.0
         self.kT = -1.0
-        self.desolvation_parameter = 0.0
+        self.desolvation_parameter = 1.0
         self.ions = []
         return
 

seekr2/modules/mmvt_analyze.py

@@ -700,6 +700,8 @@ def __init__(self, model, N_alpha_beta=None, k_alpha_beta=None,
         self.pi_alpha = None
         self.N_ij = {}
         self.R_i = {}
+        self.N_ij_unmodified = {}
+        self.R_i_unmodified = {}
         self.T = None
         self.Q = None
         self.K = None
@@ -906,6 +908,8 @@ def fill_out_data_quantities(self):
                 raise Exception("R_i_alpha should always be set.")
 
         # Now we need to try and correct the zero entries
+        self.N_ij_unmodified = deepcopy(self.N_ij)
+        self.R_i_unmodified = deepcopy(self.R_i)
         for key in self.N_ij:
             if self.N_ij[key] == 0.0:
                 i = key[0]
@@ -1051,9 +1055,10 @@ def calculate_extra_thermodynamics(self):
         Use this data sample's statistics to construct the 
         anchor free energy profile.
         """
-        self.free_energy_anchors = np.zeros(self.pi_alpha.shape[0])
+        self.free_energy_anchors = np.zeros(self.model.num_anchors)
         highest_pi_alpha = max(self.pi_alpha)
-        for i, pi_alpha_val in enumerate(self.pi_alpha):
+        for i in range(self.model.num_anchors):
+            pi_alpha_val = self.pi_alpha[i, 0]
             free_energy = -GAS_CONSTANT*self.model.temperature*np.log(
                 pi_alpha_val / highest_pi_alpha)
             self.free_energy_anchors[i] = free_energy

seekr2/tests/test_common_prepare.py

@@ -101,10 +101,14 @@ def test_move_add_delete_input_anchors(tmp_path, host_guest_mmvt_model_input):
         .upper_milestone_radius = 1.65
     host_guest_mmvt_model_input.cv_inputs[0].input_anchors[-2]\
         .lower_milestone_radius = 1.65
+    host_guest_mmvt_model_input.cv_inputs[0].input_anchors[-2]\
+        .radius = 1.75
     host_guest_mmvt_model_input.cv_inputs[0].input_anchors[-2]\
         .upper_milestone_radius = 1.85
     host_guest_mmvt_model_input.cv_inputs[0].input_anchors[-1]\
         .lower_milestone_radius = 1.85
+    host_guest_mmvt_model_input.cv_inputs[0].input_anchors[-1]\
+        .radius = 1.95
     host_guest_mmvt_model_input.root_directory = os.path.join(
         tmp_path, "host_guest_mmvt_5")
     os.chdir(TEST_DIRECTORY)