Commit 51bb0d44 authored by Pieter van Velde's avatar Pieter van Velde

import numpy as np

parent 4beda16e
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
import random
from functions_perm import add_bead, limmits, prune, enrich, maxsize
# Simulation parameters
n_beads = 100 # number of beads
n = 6 # number of angles
n_pol = 100 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 1.2 # 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
polweight_stack = np.zeros([n_beads,1])
end_end_stack = np.zeros([n_beads, 1])
pos_stack = np.zeros([n_beads, 2])
end_end_array = np.zeros([n_beads, n_pol-1])
polweight_array = np.zeros([n_beads, n_pol-1])
numprune = 0
numenrich = 0
# Make n polymers
k = 0
fails = 0
while k < n_pol-1:
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
j = 2
flag = 1
while flag == 1:
if j >= n_beads -2:
# print('pos = ', pos_stack )
polweight_array[:, k], end_end_array[:, k], polweight_stack, end_end_stack, pos_stack = maxsize(polweight_stack, end_end_stack, pos_stack)
k = k + 1
print('size stack = ', numenrich - numprune)
print('k = ', k)
# print('pol = ', polweight_stack,)
# print( 'end = ', end_end_stack)
# print('pos = ', pos_stack )
if polweight_stack.size == 0:
break
if k == n_pol - 1:
break
#add next bead
pos_stack[:,:2], polweight, end_end, fails, j = add_bead(pos_stack[:,:2], n, T_dim, polweight_stack[:,0], fails, cst_dep)
end_end_stack[j,0] = end_end
polweight_stack[j,0] = polweight
upper_lim, lower_lim = limmits(polweight_array, polweight_stack, alfa_up, j, alfa_low)
if polweight_stack[j,0] < lower_lim:
polweight_stack, end_end_stack, pos_stack = prune(polweight_stack, end_end_stack, pos_stack, j)
#print('prune')
if polweight_stack.size == 0:
flag = 0
#print('stack = 0')
numprune += 1
elif polweight_stack[j,0] > upper_lim:
if j < n_beads -1:
polweight_stack, pos_stack, end_end_stack = enrich(polweight_stack,j, pos_stack, end_end_stack)
#print('enrich')
numenrich += 1
if j == n_beads:
print('error maybe')
#print('size stack = ', numenrich - numprune)
# plt.figure()
# plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
print(end_end_array)
print(polweight_array)
print('prune = ', numprune)
print('enrich = ', numenrich)
print('size stack = ', numenrich - numprune - n_pol)
end_end_times_weight = np.multiply(end_end_array, polweight_array)
sum_end_end_times_weight = np.sum(end_end_times_weight, axis =1)
sum_weights = np.sum(polweight_array, axis=1)
exp_value = np.divide(sum_end_end_times_weight, sum_weights)
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('fails = ', fails)
#plt.figure()
# plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
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