Commit ae10accb authored by Pieter van Velde's avatar Pieter van Velde

Added compuatation for end to end length for each polymer length in main.py and functions.py

parent d0ef7aa4
import numpy as np
import sys
def add_bead(pos, n, T_dim, j, p_pol):
def add_bead(pos, n, T_dim, j, p_pol, fails):
"""Adds one bead to polymer, according to the Rosenblaum method
Parameters:
......@@ -25,7 +25,7 @@ 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
flag = 0
# Calculate new positions
pos_new = np.zeros([n, 2])
pos_new[:, 0] = pos[j - 1, 0] + np.cos(phi)
......@@ -49,10 +49,11 @@ def add_bead(pos, n, T_dim, j, p_pol):
w_tot = np.sum(w)
# Stop execution if beads can only be placed (almost) on top of other beads
if w_tot == 0:
print('Help, I have nowhere to go')
print('w_tot = ', w_tot)
sys.exit()
if w_tot < 1e-3: # BEETJE ME GEPIELD
#print('Help, I have nowhere to go')
#print('w_tot = ', w_tot)
flag = 1
fails = fails + 1
# Calculate new position of next bead
P = w / w_tot
......@@ -62,7 +63,11 @@ def add_bead(pos, n, T_dim, j, p_pol):
# Calculate weights of one polymer to account for deviation from the theoretical value
w_idx = np.where(phi == final_phi)
w_jl = w[w_idx[0]]
# print(w_jl)
w_l = w_tot
#print(w_tot)
p_pol = p_pol * w_jl / w_l
end_end = np.sqrt(pos[j, 0]**2 + pos[j,1]**2)
return pos, p_pol
\ No newline at end of file
return pos, p_pol, flag, end_end, fails
\ No newline at end of file
......@@ -4,42 +4,61 @@ from matplotlib import pyplot as plt
import time
from functions import add_bead
fails =0
start = time.time()
# Simulation parameters
n_beads = 100 #number of beads
n_beads = 250 #number of beads
n = 6 #number of angles
n_pol = 2 #number of polymers
n_pol = 10000 #number of polymers
# Physical parameters
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]
T_dim = 100 #T * (Kb / eps) # Non-dimensional temperature
p = np.zeros([n_beads, n_pol])
end_end_matrix = np.zeros([n_beads, n_pol])
# Make n polymers
for k in range(n_pol):
# Preallocate arrays
pos = np.zeros([n_beads, 2])
# Initialization
pos[1] = np.array([1, 0]) #bead 1: pos[0, 0]. bead 2: pos[1, 0]
# 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[k] = p_pol #Weight per polymer
print('P = ', p)
pos, p_pol, flag, end_end, fails = add_bead(pos, n, T_dim, j, p_pol, fails)
if flag == 1:
#print('unsucesfull polymer')
#p[:, k] = 0
k = k - 1 # throw polymer away
break
if flag == 0: # keep polymer
p[j,k] = p_pol #Weight per polymer
end_end_matrix[j,k] = end_end
#print('P = ', p)
# Plot polymer
plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
plt.show()
#plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
end = time.time()
print('fails = ', fails)
average_end_end = np.average(end_end_matrix[2::], axis=1, weights=p[2::])
a = np.arange(2,n_beads,1)
plt.loglog(a, average_end_end**2, marker= "x", markersize=7, color="k" )
plt.ylabel('R^2 (end to end distance)')
plt.xlabel('N (Number of beads)')
plt.grid(True,which="both",ls="-")
plt.ylim(1, 10000)
plt.xlim(1, 250)
print('time = ', end-start)
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