Commit 8f21d227 authored by Luuk Balkenende's avatar Luuk Balkenende

Implemented right way of calculating end to end distance

parent 3a1364ac
import numpy as np
import random
def add_bead(pos, n, T_dim, polweight, fails, cst_dep,j):
def add_bead(pos, n, T_dim, polweight, fails, cst_dep, j):
"""Adds one bead to polymer, according to the Rosenblaum method
Parameters:
......@@ -22,8 +24,7 @@ def add_bead(pos, n, T_dim, polweight, fails, cst_dep,j):
p_pol: GAAT WAARSCHIJNLIJK TOCH NOG VERANDEREN
"""
#polweight = polweight[j-1]
# polweight = polweight[j-1]
# 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
......@@ -41,103 +42,107 @@ def add_bead(pos, n, T_dim, polweight, fails, cst_dep,j):
ry[i, :] = pos_new[i, 1] - pos[0:j, 1]
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
E = np.sum(4 * 0.25*((0.8 / r) ** 12 - (0.8 / r) ** 6), axis=1)
#print(L)
w_jl = np.exp(-E / T_dim) #small weights (per angle)
w_l = np.sum(w_jl) #big weight per added bead (sum of angle weights)
E = np.sum(4 * 0.25 * ((0.8 / r) ** 12 - (0.8 / r) ** 6), axis=1)
# print(L)
w_jl = np.exp(-E / T_dim) # small weights (per angle)
w_l = np.sum(w_jl) # big weight per added bead (sum of angle weights)
if w_l == 0:
w_jl = np.ones(w_jl.shape)
w_l = np.sum(w_jl)
fails = fails+1
#print(w_l)
#print(w_l)
polweight = polweight*w_l *cst_dep #weight per polymer off by our algorithm
fails = fails + 1
# print(w_l)
# print(w_l)
polweight = polweight * w_l * cst_dep # weight per polymer off by our algorithm
# Calculate new position of next bead
P = 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)]])
end_end = np.sqrt(pos[j, 0]**2 + pos[j,1]**2)
end_end = np.sqrt(pos[j, 0] ** 2 + pos[j, 1] ** 2)
return pos, polweight, end_end, fails
def enrich(polweight_stack, pos_stack, k):
def enrich(polweight_stack, pos_stack, end_end_stack, k):
polweight_stack[:, k] = 0.5 * polweight_stack[:, k]
polweight_stack = np.insert(polweight_stack, k, polweight_stack[:,k], axis =1)
pos_stack = np.insert(pos_stack, 2*k, pos_stack[:,2*k+1], axis =1)
pos_stack = np.insert(pos_stack, 2*k, pos_stack[:,2*k+1], axis =1)
k = k + 1
return polweight_stack, pos_stack, k
polweight_stack = np.insert(polweight_stack, k, polweight_stack[:, k], axis=1)
pos_stack = np.insert(pos_stack, 2 * k, pos_stack[:, 2 * k + 1], axis=1)
pos_stack = np.insert(pos_stack, 2 * k, pos_stack[:, 2 * k + 1], axis=1)
end_end_stack = np.insert(end_end_stack, k, end_end_stack[:,k], axis=1)
k = k + 1
return polweight_stack, pos_stack, end_end_stack, k
def limits(polweight_stack, alfa_up, k, alfa_low, j):
def limits(polweight_stack, alfa_up, k, alfa_low, j):
# find all nonzero elements of stack & pw_matrix, of bead j, sum and average:
avg_weight = np.average(polweight_stack[j, :k+1])
weight3 = polweight_stack[2,k]
upper = (alfa_up * avg_weight)/weight3
avg_weight = np.average(polweight_stack[j, :k + 1])
weight3 = polweight_stack[2, k]
upper = (alfa_up * avg_weight) / weight3
lower = (alfa_low * avg_weight)/weight3
lower = (alfa_low * avg_weight) / weight3
return upper, lower
def prune(polweight_stack, pos_stack, k):
def prune(polweight_stack, pos_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m):
b = random.random()
# print('b =', b)
if b < 0.5:
polweight_stack[:,k] = polweight_stack[:,k] * 2
polweight_stack[:, k] = polweight_stack[:, k] * 2
else:
polweight_stack = np.delete(polweight_stack, (k), axis = 1)
pos_stack = np.delete(pos_stack, [2*k, 2*k+1], axis = 1)
k = k-1
polweight_array = np.c_[polweight_array, polweight_stack[:, k]]
#polweight_array[:, m] = polweight_stack[:, k]
polweight_stack = np.delete(polweight_stack, k, axis=1)
pos_array = np.c_[pos_array, pos_stack[:, 2*k]]
pos_array = np.c_[pos_array, pos_stack[:, 2*k+1]]
# pos_array[:, 2*m] = pos_stack[:, 2*k]
# pos_array[:, 2*m+1] = pos_stack[:, 2*k+1]
pos_stack = np.delete(pos_stack, [2 * k, 2 * k + 1], axis=1)
end_end_array = np.c_[end_end_array, end_end_stack[:, k]]
end_end_stack = np.delete(end_end_stack, k, axis=1)
k = k - 1
m = m + 1
return pos_stack, polweight_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m
return pos_stack, polweight_stack, k
def maxsize(polweight_stack, end_end_stack, pos_stack):
a = polweight_stack[:, 0]
a = a.flatten()
b = end_end_stack[:,0]
b = end_end_stack[:, 0]
b = b.flatten()
polweight_stack = np.delete(polweight_stack, (0), axis = 1)
end_end_stack = np.delete(end_end_stack, (0), axis = 1)
pos_stack = np.delete(pos_stack, (0), axis = 1)
pos_stack = np.delete(pos_stack, (0), axis = 1)
polweight_stack = np.delete(polweight_stack, (0), axis=1)
end_end_stack = np.delete(end_end_stack, (0), axis=1)
pos_stack = np.delete(pos_stack, (0), axis=1)
pos_stack = np.delete(pos_stack, (0), axis=1)
return a, b, polweight_stack, end_end_stack, pos_stack
def stack_init(n_beads):
# initilize stacks to for storing data
polweight_stack = np.zeros([n_beads,1])
# initilize stacks to for storing data
polweight_stack = np.zeros([n_beads, 1])
end_end_stack = np.zeros([n_beads, 1])
pos_stack = np.zeros([n_beads, 2])
polweight_stack[0], polweight_stack[1] = 1, 1
# Preallocate arrays
#pos = np.zeros([n_beads, 2])
# pos = np.zeros([n_beads, 2])
# Initialization
pos_stack[1] = np.array([1, 0]) # bead 1: pos[0, 0]. bead 2: pos[1, 0]
#AddBead(Polymer,Polweight,L)
#Add third bead
# AddBead(Polymer,Polweight,L)
# Add third bead
return polweight_stack, end_end_stack, pos_stack
def exp_calc(weigths, exp_values):
end_end_times_weight = np.multiply(exp_values, weigths)
sum_end_end_times_weight = np.sum(end_end_times_weight, axis =1)
sum_end_end_times_weight = np.sum(end_end_times_weight, axis=1)
sum_weights = np.sum(weigths, axis=1)
exp_value = np.divide(sum_end_end_times_weight, sum_weights)
return exp_value
\ No newline at end of file
......@@ -2,65 +2,94 @@
import numpy as np
from matplotlib import pyplot as plt
import random
from functionsperm2 import add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc
from functions import add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc
# Simulation parameters
n_beads = 250 # number of beads
n = 6 # number of angles
n_pol = 100 # number of polymers
alfa_up = 100 # alfa to compute upper limit
alfa_low = 0 # alfa to compute lwoer limit
cst_dep = 1/(0.75*n) # multiplication constant for L dependence perm
n_pol = 1000 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 1 # alfa to compute lwoer limit
cst_dep = 1 / (0.75 * n) # multiplication constant for L dependence perm
# Physical parameters
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
m = np.arange(0, 2*n_pol, 2)
pos_stack = np.zeros([n_beads, n_pol*2])
m = np.arange(0, 2 * n_pol, 2)
pos_stack = np.zeros([n_beads, n_pol * 2])
pos_stack[1, m] = 1
#print(pos_stack)
polweight_stack = np.zeros([n_beads,n_pol])
polweight_stack[0:2,:] = 1
end_end_stack = np.zeros([n_beads,n_pol])
# print(pos_stack)
polweight_stack = np.zeros([n_beads, n_pol])
polweight_stack[0:2, :] = 1
polweight_array = np.zeros([n_beads])
end_end_stack = np.zeros([n_beads, n_pol])
end_end_array = np.zeros([n_beads, 1])
pos_array = np.zeros([n_beads, 2])
enrich = 0
prune = 0
enrich1 = 0
prune1 = 0
flag = 1
fails = 0
j = 2
m = 0
while flag == 1:
k=0
while k < int(pos_stack.shape[1]/2) :
pos_stack[:, 2*k:2*k+2], polweight_stack[j, k], end_end_stack[j,k], fails = add_bead(pos_stack[:, 2*k:2*k+2], n, T_dim, polweight_stack[j-1,k], fails, cst_dep,j)
k = 0
flag2 = 1
while flag2 == 1:
#print('end =', end_end_stack)
#print(k)
#print('prune = ', prune1)
pos_stack[:, 2 * k:2 * k + 2], polweight_stack[j, k], end_end_stack[j, k], fails = add_bead(
pos_stack[:, 2 * k:2 * k + 2], n, T_dim, polweight_stack[j - 1, k], fails, cst_dep, j)
uplim, lowlim = limits(polweight_stack, alfa_up, k, alfa_low, j)
if polweight_stack[j, k] < lowlim:
prune = (prune + 1)/2
pos_stack, polweight_stack, k = prune(polweight_stack, pos_stack, k)
prune1 = prune1+1#(prune1 + 1) / 2
pos_stack, polweight_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m = prune(
polweight_stack, pos_stack, pos_array, polweight_array, end_end_stack, end_end_array, k, m)
elif polweight_stack[j, k] > uplim:
enrich = enrich + 1
polweight_stack, pos_stack, k = enrich(polweight_stack, pos_stack, k)
k = k +1
enrich1 = enrich1 + 1
polweight_stack, pos_stack, end_end_stack, k = enrich(polweight_stack, pos_stack, end_end_stack, k)
k = k + 1
if k >= int(pos_stack.shape[1] / 2):
flag2 = 0
#print('shape = ', pos_stack.shape[1]/2)
#print('shape = ', polweight_stack.shape[1] / 2)
#print('k = ', k)
#print('enrich = ', enrich1)
if pos_stack.size == 0 or j== n_beads-1 :
if pos_stack.size == 0 or j == n_beads - 1:
flag = 0
print('stack equal to zero at number of beads = ', j )
j=j+1
print('stack equal to zero at number of beads = ', j)
j = j + 1
print(j)
end_end_array = np.c_[end_end_array, end_end_stack]
polweight_array = np.c_[polweight_array, polweight_stack]
exp_value = exp_calc(polweight_array, end_end_array)
a = np.arange(2, n_beads - 1, 1)
# if np.count_nonzero(end_end_stack) != 0:
# print('sukkel')
# #print('end', end_end_stack)
# exp_value = exp_calc(polweight_stack, end_end_stack)
# a = np.arange(2, n_beads - 1, 1)
# print(np.count_nonzero(end_end_array))
# if np.count_nonzero(end_end_array) != 0:
# print('wiehoee')
# times_weight = np.multiply(end_end_array, polweight_array)
#
# print('times_weight = ', times_weight)
exp_value = exp_calc(polweight_stack, end_end_stack)
a = np.arange(2,n_beads-1,1)
plt.figure()
plt.loglog(a, exp_value[2:n_beads-1]**2, marker= "x", markersize=7, color="k" )
plt.loglog(a, exp_value[2:n_beads - 1] ** 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.grid(True, which="both", ls="-")
plt.ylim(1, 10000)
plt.xlim(1, 250)
print(polweight_stack.shape)
print('enrich = ', enrich)
print('prune = ', prune)
print('enrich = ', enrich1)
print('prune = ', prune1)
plt.show()
\ No newline at end of file
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