Commit 3a1364ac authored by Pieter van Velde's avatar Pieter van Velde

Change implementation of the PERM algorithm by storing the data of low weight polymers

parent 8ef5e54a
import numpy as np
import sys
def add_bead(pos, n, T_dim, j, p_pol, fails):
import random
def add_bead(pos, n, T_dim, polweight, fails, cst_dep,j):
"""Adds one bead to polymer, according to the Rosenblaum method
Parameters:
......@@ -22,10 +21,13 @@ def add_bead(pos, n, T_dim, j, p_pol, fails):
Position of n_beads beads in x,y direction
p_pol: GAAT WAARSCHIJNLIJK TOCH NOG VERANDEREN
"""
#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
flag = 0
# Calculate new positions
pos_new = np.zeros([n, 2])
pos_new[:, 0] = pos[j - 1, 0] + np.cos(phi)
......@@ -39,35 +41,103 @@ def add_bead(pos, n, T_dim, j, p_pol, fails):
ry[i, :] = pos_new[i, 1] - pos[0:j, 1]
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
# Calculate energy of different angles
E = np.sum(4 * ((1 / r) ** 12 - (1 / r) ** 6), axis=1)
# Calculate (normalized) weights
w = np.exp(-E / T_dim)
w[np.isnan(w)] = 0
w_tot = np.sum(w)
# Stop execution if beads can only be placed (almost) on top of other beads
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
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
# Calculate new position of next bead
P = w / w_tot
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)
return pos, polweight, end_end, fails
# 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)
def enrich(polweight_stack, pos_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
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 polweight_stack, pos_stack, k
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
lower = (alfa_low * avg_weight)/weight3
return upper, lower
def prune(polweight_stack, pos_stack, k):
b = random.random()
if b < 0.5:
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
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 = 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)
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])
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])
# 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
return polweight_stack, end_end_stack, pos_stack
return pos, p_pol, flag, end_end, fails
\ No newline at end of file
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_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
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
import time
import random
from functionsperm2 import add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc
from functions import add_bead
fails =0
start = time.time()
# Simulation parameters
n_beads = 250 #number of beads
n = 6 #number of angles
n_pol = 10000 #number of polymers
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
# Physical parameters
T_dim = 100 #T * (Kb / eps) # Non-dimensional temperature
T_dim = 10 # 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):
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])
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
enrich = 0
prune = 0
flag = 1
#print('P = ', p)
# Plot polymer
#plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
fails = 0
j = 2
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)
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)
end = time.time()
elif polweight_stack[j, k] > uplim:
enrich = enrich + 1
polweight_stack, pos_stack, k = enrich(polweight_stack, pos_stack, k)
k = k +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('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" )
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.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()
print(polweight_stack.shape)
print('enrich = ', enrich)
print('prune = ', prune)
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