From 11ae29dea8b663d7c54c79ff0199cfaa6204d47e Mon Sep 17 00:00:00 2001
From: Pim Nelissen <pim.nelissen@protonmail.com>
Date: Mon, 27 Oct 2025 15:40:31 +0100
Subject: [PATCH] update RubberBand to add biasing feature

---
 RubberBand.py | 33 +++++++++++++++++++++------------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/RubberBand.py b/RubberBand.py
index 7ee6e76..6fc9f07 100644
--- a/RubberBand.py
+++ b/RubberBand.py
@@ -4,39 +4,48 @@ from Link import Link
 
 class RubberBand:
     # class attributes
-    k_B = 1
-    
-    def __init__(self, N, a=1., T=1.):
+ 
+    def __init__(self, N, a=1., kBT=1., biased=False, force=0.):
         """
         RubberBand is a class that can simulate a rubber band with N
         incompressible links of length a.
 
-        Parameters:
-            N (int)     The number of links.
-            a (float)   The fixed length of links.
-            T (float)   Temperature of the system.
+        Parameters
+        ----------
+        N : int
+            The number of links.
+        a : float, optional
+            The fixed length of the links. The default is 1..
+        kBT : float, optional
+            Temperature of the system (arb. units). The default is 1..
+        force : float, optional
+            Force applied to the RubberBand.
+            The default is 0, which is an unbiased RubberBand.
         """
 
         self.N = N
         self.a = a
-        self.T = T
-        self.links = self.__sample_links()
+        self.kBT = kBT
+        self.links = self.__sample_links(force=force)
 
     def w(self, f, l):
-        beta = (self.k_B*self.T)**(-1)
+        beta = (self.kBT)**(-1)
         return np.exp(beta*f*l)
     
     @property
     def length(self):
         return self.a * np.sum([l.direction for l in self.links])
     
-    def __sample_links(self):
+    def __sample_links(self, force):
         """
         Sample N Link objects with a random direction.
         """
 
         samples = np.random.random_sample((self.N,))
         directions = np.ones(samples.shape)
-        directions[samples < 0.5] = -1.
+        
+        p_negative = 1 - 0.5*(1 + np.tanh(force*self.a/self.kBT))
+
+        directions[samples < p_negative] = -1.
         directions = directions.astype(int)
         return [Link(d) for d in directions]
\ No newline at end of file