Combining previous commits

Author: Angelika Genya Hirsch

chromo/mc/__init__.py

@@ -109,7 +109,7 @@ def _polymer_in_field(
         (default = None)
     """
     np.random.seed(random_seed)
-    print("output init "+ str(output_dir))
+    #print("output init "+ str(output_dir))
     #print(log_directory)
     if mc_move_controllers is None:
         mc_move_controllers = all_moves(
@@ -143,6 +143,8 @@ def _polymer_in_field(
                 temperature_adjust_factor = temp_schedule.logarithmic_decrease(mc_count, num_saves)
             elif temperature_schedule == "decreasing stepwise":
                 temperature_adjust_factor = temp_schedule.decreasing_stepwise(mc_count, num_saves)
+            elif temperature_schedule == "sin decrease":
+                temperature_adjust_factor = temp_schedule.sin_decrease(mc_count, num_saves)
             elif temperature_schedule == "no schedule":
                 temperature_adjust_factor = temp_schedule.no_schedule()
             else:
@@ -153,10 +155,14 @@ def _polymer_in_field(
         if lt_schedule is not None:
             if lt_schedule == "logarithmic increase":
                 lt_change = twist_schedule.logarithmic_increase(mc_count, num_saves)
+            elif lt_schedule == "exponential increase":
+                lt_change = twist_schedule.exponential_increase(mc_count, num_saves)
             elif lt_schedule == "linear increase":
                 lt_change = twist_schedule.linear_increase(mc_count, num_saves)
             elif lt_schedule == "increasing stepwise":
-                lt_change = twist_schedule.increasing_stepwise(mc_count, num_saves)
+                lt_change = twist_schedule.step_wise_increase(mc_count, num_saves)
+            elif lt_schedule == "increasing sawtooth":
+                lt_change = twist_schedule.increasing_sawtooth(mc_count, num_saves)
             elif lt_schedule == "no schedule":
                 lt_change = twist_schedule.no_schedule()
             else:

chromo/mc/mc_sim.pyx

@@ -76,7 +76,6 @@ cpdef void mc_sim(
     cdef long i, j, k, n_polymers
     cdef list active_fields
     cdef poly
-    #cdef list accept_reject
     np.random.seed(random_seed)
     if field.confine_type == "":
         active_fields = [poly.is_field_active() for poly in polymers]
@@ -166,6 +165,7 @@ cpdef void mc_step(
             exp_dE = 1
 
     if (<double>rand() / RAND_MAX) < exp_dE/temperature_adjust_factor: #temperature variation
+        #print(<double>rand() / RAND_MAX)
         adaptible_move.accept(
             poly, dE, inds, n_inds, log_move=True, log_update=True
         )

chromo/polymers.pyx

@@ -1162,8 +1162,8 @@ cdef class SSWLC(PolymerBase):
 
     cdef double bead_pair_dE_poly_forward(
         self,
-        double[:] r_0, #posoiton of first bead in move
-        double[:] r_1, # position of the escond bead in the move
+        double[:] r_0, #positon of first bead in move
+        double[:] r_1, # position of the second bead in the move
         double[:] test_r_1,
         double[:] t3_0,
         double[:] t3_1,
@@ -2030,7 +2030,6 @@ cdef class SSTWLC(SSWLC):
             0.5 * self.eps_perp[ind] * vec_dot3(dr_perp, dr_perp) +
             0.5 * self.eps_twist[ind] * omega**2
         )
-        #print("does E_pair_with_twist get used")
         return E
 
     cdef double bead_pair_dE_poly_forward_with_twist(

chromo/run_simulation.py

@@ -111,7 +111,7 @@
     controller=ctrl.SimpleControl
 )
 
-num_snapshots = 560
+num_snapshots = 200
 # num_snapshots = 1000 # try 1000 and average for each set of 100, depending on pre-equilibration steps
 # count number of accepted moves for different conditions
 mc_steps_per_snapshot = 40000
@@ -127,11 +127,10 @@
     num_saves = num_snapshots,
     bead_amp_bounds = amp_bead_bounds,
     move_amp_bounds = amp_move_bounds,
-    #output_dir = '/scratch/users/ahirsch1', # for running on sherlock
     output_dir = 'output',
     mc_move_controllers = moves_to_use,
-    temperature_schedule = "no schedule",
-    lt_schedule = "linear increase"
+    temperature_schedule = "logarithmic decrease",
+    lt_schedule = "no schedule"
 )
 
 output_files = os.listdir(latest_sim)

chromo/twist_schedule.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 
 class Schedule:
     """Class representation of the simulated annealing schedule.
@@ -16,25 +18,90 @@ def to_file(self, path):
         with open(path, 'w') as f:
             f.write(self.name)
 
-def increasing_stepwise(current_step, total_steps):
-    step_size = (1 - 0.1) / total_steps
-    value = max(0.1, 1 - step_size * current_step)
 
-    return value
+# dynamic structure factor
+# density correlation function at different length-scales
+# fourier transform of density-density correlation function
+# how much largest length-scale of polymer is changing
+# correlation function of radius of gyration
+# how density fluctuation within polymer correlate in time
+# plot rmsd for different twist values
+
+
 def logarithmic_increase(current_step, total_steps):
-    #value = 5
-    value =  (current_step/total_steps) * 100
-    return value
+    if current_step == 0:
+        return 0
+    ratio = current_step / total_steps
+    result = np.log(ratio / 2) / 5.5 * 100 + 113
+    if result > 100:
+        result = 100
+    if result < 0:
+        result = 0
+    return result
 
+def exponential_increase(current_step, total_steps):
+    result = np.exp(current_step/total_steps) * 58
+    return result - 58
+# for lp
 def linear_increase(current_step, total_steps):
-    value =  (current_step/total_steps) * 100
+    max_value = 100
+    slope = max_value / total_steps
+    value = slope * current_step
     return value
+
+
+# for lp
+def step_wise_increase(current_step, total_steps):
+    num_blocks = 25
+    max_height = 100
+    step_height = max_height / num_blocks  # each step has a height of 10
+    step_length = total_steps / num_blocks  # so 20 is the length
+    division = current_step / step_length  # so if we are at snapshot 105 we get 5.11
+    ceiling = np.ceil(division)  # we are currently on the 6th step
+    result = step_height * ceiling
+    return result
+
+
+def increasing_sawtooth(current_step, total_steps):
+    num_blocks = 10
+    max_height = 100
+    step_height = max_height / num_blocks  # each step has a height of 10
+    step_length = total_steps / num_blocks  # so 20 is the length
+    division = current_step / step_length  # so if we are at snapshot 105 we get 5.11
+    ceiling = np.ceil(division)  # we are currently on the 6th step
+    result = step_height * ceiling  # this is the height we are at in the 6th block
+
+    num_sections = num_blocks * 2  # if we are at 105
+    section_division = np.floor(
+        current_step / (total_steps / num_sections))  # 10.5 is the section from 20 total sections
+    if section_division % 2 == 0:
+        result = result + (current_step % num_sections) * 4 - section_division - 14
+    else:
+        result = result - (current_step % num_sections) * 2.5 - section_division + 46
+    if ceiling == num_blocks and section_division % 2 != 0:
+        result = max_height
+    if result > max_height:
+        result = max_height
+    if result < 0:
+        result = 0
+    return result
+
+def increasing_sawtooth(current_step, total_steps):
+    num_blocks = 10
+    max_height = 100
+    step_height = max_height/num_blocks # each step has a height of 10
+    step_length = total_steps / num_blocks # so 20 is the length
+    division = current_step/step_length # so if we are at snapshot 105 we get 5.11
+    ceiling = np.ceil(division) # we are currently on the 6th step
+    result = step_height * ceiling # this is the height we are at in the 6th block
+
+    num_sections = num_blocks * 2 #if we are at 105
+    section_division = np.floor(current_step/(total_steps/num_sections)) # 10.5 is the section from 20 total sections
+    if section_division%2 == 0:
+        result = result + (current_step % num_sections) * 2 - section_division - 20
+    else:
+        result = result - (current_step%num_sections)*2 - section_division + 20
+    return result
+
 def no_schedule():
     return 100
-#dynamic structure factor
-#density correlation function at differnet lengthscales
-#fouirer transform of density density correlation function
-#how much largest lengthscale of polymer is changing
-#correleation function of radius of gyration
-#how density fluctluation within polymer correlate in time
-#plot rmsd for different twist values

chromo/util/mc_stat.py

@@ -100,8 +100,8 @@ def create_log_file(self, ind: int):
             Index with which to append log file name
         """
         log_path = self.log_dir+"/"+self.log_file_prefix+str(ind)+".csv"
-        print("log dir in mc_stat "+str(self.log_dir))
-        print("prefix in mcstat " + str(self.log_file_prefix))
+        #print("log dir in mc_stat "+str(self.log_dir))
+        #print("prefix in mcstat " + str(self.log_file_prefix))
         column_labels = [
             "snapshot",
             "iteration",

chromo/util/temperature_schedule.py

@@ -18,15 +18,34 @@ def to_file(self, path):
         with open(path, 'w') as f:
             f.write(self.name)
 
+def sin_decrease(current_step, total_steps):
+    slope = (-1.3) / total_steps
+    linear_value = 1 + slope * current_step
+    result = np.sin(current_step/total_steps *100) 
+    result = result *.1 + linear_value *.6 +.32
+    if result > 1:
+        result = 1
+    if result < 0.1:
+        result = 0.1
+    return result * 10
+
 
 def decreasing_stepwise(current_step, total_steps):
-    step_size = (1 - 0.1) / total_steps
-    value = max(0.1, 1 - step_size * current_step)
+    raw_fraction = 1 - current_step/(total_steps + 150)
+    fraction = raw_fraction * 10
+    step_fraction = np.floor(fraction)/10
+    return step_fraction
 
-    return value
 def logarithmic_decrease(current_step, total_steps):
-    value = 1 - (0.9 ** (current_step / total_steps))
-    return value
+    if current_step == 0:
+    	return 1
+    ratio = current_step/total_steps
+    result = -1 * np.log(ratio/2)/5.5 -0.1
+    if result > 1:
+        result = 1
+    if result < 0:
+        result = 0
+    return result * 10
 
 def linear_decrease(current_step, total_steps):
     slope = (0.1 - 1) / total_steps