Commit e8ba5e23 authored by Luuk Balkenende's avatar Luuk Balkenende

Code now capable of calculation multiple polymers

parent 80553a41
......@@ -2,7 +2,7 @@ import numpy as np
import sys
def add_bead(pos, n, T_dim, j, w_jl, w_l):
def add_bead(pos, n, T_dim, j, p_pol):
# Standard phi, KAN NAAR MAIN.PY
phi_0 = np.linspace(0, 2 * np.pi - (2 / n) * np.pi, n)
phi = phi_0 + (1 - 2 * np.random.rand(n)) * np.pi / n
......@@ -27,27 +27,9 @@ def add_bead(pos, n, T_dim, j, w_jl, w_l):
# Calculate (normalized) weights
w = np.exp(-E / T_dim)
w[np.isnan(w)] = 0
<<<<<<< HEAD
w = np.exp(-E/T_dim)
w[np.isnan(w)]=0
w_tot = np.sum(w)
P = w/w_tot
final_phi = np.random.choice(phi, p=P)
new = np.array([[pos[1,0] + np.cos(final_phi), pos[1,1] + np.sin(final_phi)]])
pos = np.r_[pos, new]
=======
w_tot = np.sum(w)
# Stop execution if beads can only be placed (almost) on top of other beads
>>>>>>> 1b70cc076f246da4e36ffb50e1759b96fb1f921b
if w_tot == 0:
print('Help, I have nowhere to go')
print('w_tot = ', w_tot)
......@@ -58,9 +40,10 @@ def add_bead(pos, n, T_dim, j, w_jl, w_l):
final_phi = np.random.choice(phi, p=P)
pos[j] = np.array([[pos[j - 1, 0] + np.cos(final_phi), pos[j - 1, 1] + np.sin(final_phi)]])
# Calculate weights of one polymer to account for deviation form the theoretical value
# Calculate weights of one polymer to account for deviation from the theoretical value
w_idx = np.where(phi == final_phi)
w_jl = w_jl * w[w_idx[0]]
w_l = w_l * w_tot
w_jl = w[w_idx[0]]
w_l = w_tot
p_pol = p_pol * w_jl / w_l
return pos, w_jl, w_l
\ No newline at end of file
return pos, p_pol
\ No newline at end of file
......@@ -9,27 +9,32 @@ from functions import add_bead
# Simulation parameters
n_beads = 100 #number of beads
n = 6 #number of angles
n_pol = 2 #number of polymers
# Physical parameters
T_dim = 0.1 #T * (Kb / eps) # Non-dimensional temperature
T_dim = 10 #T * (Kb / eps) # Non-dimensional temperature
# Preallocate arrays
pos = np.zeros([n_beads, 2])
p = np.zeros(n_pol)
# Initialization
pos[1] = np.array([1, 0]) #bead 1: pos[0, 0]. bead 2: pos[1, 0]
w_jl = 1 #start value of weight w_j
w_l = 1 #start value of weight W
# Make polymer of lenght n_beads
for j in range(2,n_beads):
# Make n polymers
for k in range(n_pol):
pos, w_jl, w_l = add_bead(pos, n, T_dim, j, w_jl, w_l)
# Make polymer of lenght n_beads
#w_jl = 1 # start value of weight w_j
#w_l = 1 # start value of weight W
p_pol = 1
for j in range(2,n_beads):
pos, p_pol = add_bead(pos, n, T_dim, j, p_pol)
P = w_jl/w_l #Weight per polymer
p[k] = p_pol #Weight per polymer
print('P = ', p)
# Plot polymer
plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
plt.show()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment