Commit 22916a33 authored by Luuk Balkenende's avatar Luuk Balkenende

Radius of gyration calculation added to processfunctions.py

parent 3a520a3a
......@@ -17,8 +17,8 @@ random_walk = False # Put on True if random walk is desired
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
# Simulate configuration of polymers and their weights
end_end_array, polweight_array, pol_size = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk)
end_end_array, polweight_array, pos_array = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk)
# Process configuration of polymers and their weights into figures (end2end distance against n_beads)
# and preform bootstrap to calculate errors.
exp_calc(polweight_array, end_end_array, n_beads, pol_size)
exp_calc(polweight_array, end_end_array, n_beads, pos_array)
......@@ -3,7 +3,7 @@ import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline
def exp_calc(polweight_array, end_end_array, n_beads):
def exp_calc(polweight_array, end_end_array, n_beads, pos_array):
"""Calculates the end to end distances with taking into account their weights, and
makes a figure of the mean values of those distances against number of beads.
......@@ -22,13 +22,17 @@ def exp_calc(polweight_array, end_end_array, n_beads):
"""
exp_value_mean, exp_value_std = bootstrap(polweight_array, end_end_array, n_beads)
rad_gyr = radiusgyration(pos_array)
rad_gyr_mean, rad_gyr_std = bootstrap(polweight_array, rad_gyr, n_beads)
popsize = polsize(end_end_array)
a = np.arange(0, n_beads - 1, 1)
b = np.arange(0, n_beads, 1)
exp_value[0:2] = [0, 1]
exp_value_mean[0:2] = [0, 1]
rad_gyr_mean[0] = 0
font = {'family' : 'DejaVu Sans',
'weight' : 'normal',
......@@ -116,4 +120,48 @@ def polsize(end_end_array):
popsize = np.count_nonzero(end_end_array, axis=1)
print(popsize)
print(popsize.shape)
return popsize
\ No newline at end of file
return popsize
def radiusgyration(pos_array):
pos_array = np.delete(pos_array, [0, 1], axis=1)
#pos_array = np.delete(pos_array, 0, axis=0)
[nx, ny] = pos_array.shape
print('nx',nx)
#print('ny',ny)
#print('pos_array', pos_array)
imax = nx
jmax = ny
rad_gyr = np.zeros([imax, int(jmax/2)])
j = 0
while j < jmax:
i = 1
imax = np.count_nonzero(pos_array[:, j], axis=0)
#print('imax = ', imax)
#print('j = ', j)
while i < imax:
N = i + 1
xk = pos_array[:i+1, j]
#print(xk)
#exit()
yk = pos_array[:i+1, j+1]
xmean = np.sum(xk)/N
ymean = np.sum(yk)/N
xr = xk - xmean
yr = yk - ymean
rad_gyr[i, int(j/2)] = 1/N * np.sum(np.power(xr, 2) + np.power(yr, 2))
#print('rad_gyr ij', rad_gyr[i, int(j/2)])
#print('i = ', i, ' j = ', j)
i = i + 1
j = j + 2
print('shape gyros = ', rad_gyr.shape)
#rad_gyr = np.insert(rad_gyr, 0, rad_gyr[0, :], axis=0)
rad_gyr_zeros = np.zeros(rad_gyr.shape[0])
rad_gyr = np.insert(rad_gyr, 0, rad_gyr_zeros, axis=1)
print('shape gyros = ', rad_gyr.shape)
return rad_gyr
\ No newline at end of file
......@@ -71,17 +71,17 @@ def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk):
print('n_pol = ', polweight_stack.shape[1])
print(j)
pol_size[j] = polweight_stack.shape[1]
j = j + 1
# Append the two end2end matrices, and their weights
end_end_array = np.c_[end_end_array, end_end_stack]
polweight_array = np.c_[polweight_array, polweight_stack]
pos_array = np.c_[pos_array, pos_stack]
print('n_enrich = ', n_enrich)
print('n_prune = ', n_prune)
return end_end_array, polweight_array, pol_size
return end_end_array, polweight_array, pos_array
def allocate(n_beads, n_pol):
"""Allocates necessary arrays to store positions and weights of beads of polymers
......@@ -222,6 +222,12 @@ def perm(pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack,
uplim, lowlim = limits(polweight_stack, alfa_up, k, alfa_low, j)
# Conditions for enriching
if polweight_stack[j, k] > uplim:
n_enrich = n_enrich + 1
polweight_stack, pos_stack, end_end_stack, k \
= enrich(polweight_stack, pos_stack, end_end_stack, k)
# Conditions for pruning
if polweight_stack[j, k] < lowlim:
n_prune = n_prune + 1 # (prune1 + 1) / 2
......@@ -229,12 +235,6 @@ def perm(pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack,
= prune(polweight_stack, pos_stack, pos_array, polweight_array,
end_end_stack, end_end_array, k, m)
# Conditions for enriching
elif polweight_stack[j, k] > uplim:
n_enrich = n_enrich + 1
polweight_stack, pos_stack, end_end_stack, k \
= enrich(polweight_stack, pos_stack, end_end_stack, k)
return pos_stack, pos_array, polweight_stack, polweight_array, \
end_end_stack, end_end_array, j, k, m, n_prune, n_enrich
......@@ -263,7 +263,7 @@ def add_bead(pos, n, T_dim, polweight, j, random_walk):
end_end: float > 0
Value which gives the distance from first bead to newly added bead
"""
flag = 0
# Calculate (random) possible angles of next bead
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
......@@ -292,10 +292,13 @@ def add_bead(pos, n, T_dim, polweight, j, random_walk):
if w_l == 0:
w_jl = np.ones(w_jl.shape)*10**(-50)
w_l = np.sum(w_jl)
flag = 1
cst_dep = 1 / (0.75 * n) # multiplication constant for L dependence perm
polweight = polweight * w_l * cst_dep # weight per polymer off by our algorithm
if flag == 1:
polweight = 0
# Calculate new position of next bead
P = w_jl / w_l
final_phi = np.random.choice(phi, p=P)
......
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