Commit 338f704b authored by Pieter van Velde's avatar Pieter van Velde

added two booleans in order to switch on/off pure random walk and PERM

parent 06cc6eb4
import numpy as np
import random
def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim):
def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim, PERM, random_walk):
# Allocate necessary arrays
pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack,\
......@@ -30,7 +30,7 @@ def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim):
end_end_stack, end_end_array, j, k, m, n_prune, n_enrich = \
perm(pos_stack, pos_array, polweight_stack, polweight_array,
end_end_stack, end_end_array, j, k, m, n, n_prune, n_enrich,
T_dim, cst_dep, alfa_up, alfa_low)
T_dim, cst_dep, alfa_up, alfa_low, PERM, random_walk)
k = k + 1
if k >= int(pos_stack.shape[1] / 2):
......@@ -72,33 +72,39 @@ def init(n_pol, pos_stack, polweight_stack):
return pos_stack, polweight_stack
def perm(pos_stack, pos_array, polweight_stack, polweight_array, end_end_stack, end_end_array,
j, k, m, n, n_prune, n_enrich, T_dim, cst_dep, alfa_up, alfa_low):
j, k, m, n, n_prune, n_enrich, T_dim, cst_dep, alfa_up, alfa_low, PERM, random_walk):
# Avoid PERM and random_walk are true together
if random_walk == True and PERM == True:
print('ERROR: random_walk and PERM cannot be True both, change this in simulations parameters.')
exit()
# Add next bead to polymer
pos_stack[:, 2 * k:2 * k + 2], polweight_stack[j, k], end_end_stack[j, k] \
= add_bead(pos_stack[:, 2 * k:2 * k + 2], n, T_dim, polweight_stack[j - 1, k], cst_dep, j)
= add_bead(pos_stack[:, 2 * k:2 * k + 2], n, T_dim, polweight_stack[j - 1, k], cst_dep, j, random_walk)
# Calculate upper and lower limits for PERM algorithm
uplim, lowlim = limits(polweight_stack, alfa_up, k, alfa_low, j)
# Conditions for pruning
if polweight_stack[j, k] < lowlim:
n_prune = n_prune + 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)
# 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)
if PERM == True:
uplim, lowlim = limits(polweight_stack, alfa_up, k, alfa_low, j)
# Conditions for pruning
if polweight_stack[j, k] < lowlim:
n_prune = n_prune + 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)
# 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
def add_bead(pos, n, T_dim, polweight, cst_dep, j):
def add_bead(pos, n, T_dim, polweight, cst_dep, j, random_walk):
"""Adds one bead to polymer, according to the Rosenblaum method
Parameters:
......@@ -141,10 +147,14 @@ def add_bead(pos, n, T_dim, polweight, cst_dep, j):
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 random_walk == True:
w_jl = np.ones(n)
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)
w_jl = np.ones(w_jl.shape) *10**-50
w_l = np.sum(w_jl)
print('fail')
polweight = polweight * w_l * cst_dep # weight per polymer off by our algorithm
......
......@@ -2,21 +2,23 @@
import numpy as np
from matplotlib import pyplot as plt
import random
from functions import calc_pol, allocate, init, perm, add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc
from functions_perm import calc_pol, allocate, init, perm, add_bead, limits, prune, enrich, maxsize, stack_init, exp_calc
# Simulation parameters
n_beads = 100 # number of beads
n = 6 # number of angles
n_pol = 10 # number of polymers
n_pol = 500 # 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
PERM = False # boolean
random_walk = False # boolean
# Physical parameters
T_dim = 10 # T * (Kb / eps) # Non-dimensional temperature
T_dim = 100 # T * (Kb / eps) # Non-dimensional temperature
# Calculate configuration of polymers and their weights
end_end_array, polweight_array = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim)
end_end_array, polweight_array = calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, cst_dep, T_dim, PERM, random_walk)
exp_value = exp_calc(polweight_array, end_end_array)
a = np.arange(2, n_beads - 1, 1)
......
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