Commit 414c7675 authored by Luuk Balkenende's avatar Luuk Balkenende

Code to make n large polymers. Added weights.

parent 30f2728a
import numpy as np
import sys
def add_bead(pos, n, T_dim):
phi = np.linspace(0, 2*np.pi - (2/n)*np.pi, n)
def add_bead(pos, n, T_dim, j, w_jl, w_l):
phi_1 = np.linspace(0, 2 * np.pi - (2 / n) * np.pi, n)
E = np.zeros(n)
for i in range(n):
phi = phi_1 + (1-2*np.random.rand(n))*np.pi/n
pos_new = np.zeros(2)
pos_new[0] = pos[1,0] + np.cos(phi[i])
pos_new[1] = pos[1,1] + np.sin(phi[i])
rx = pos_new[0] - pos[:, 0]
ry = pos_new[1] - pos[:, 1]
r = np.sqrt(rx**2 + ry**2)
r[r==0] = np.nan
E[i] = np.sum(4*((1/r)**12 - (1/r)**6))
w = np.exp(-E/T_dim)
w[np.isnan(w)]=0
pos_new[0] = pos[j-1, 0] + np.cos(phi[i])
pos_new[1] = pos[j-1, 1] + np.sin(phi[i])
rx = pos_new[0] - pos[0:j, 0]
ry = pos_new[1] - pos[0:j, 1]
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
E[i] = np.sum(4 * ((1 / r) ** 12 - (1 / r) ** 6))
#print('E = ', E)
w = np.exp(-E / T_dim)
#if w == NaN:
print('w = ', w)
#w[np.isnan(w)] = 0
w_tot = np.sum(w)
if w_tot == 0:
print('Help, I have nowhere to go')
print('w_tot = ', w_tot)
sys.exit()
P = w/w_tot
final_phi = np.random.choice(phi, p=P)
new = np.array([[pos[1,0] + np.cos(final_phi), pos[1,1] + np.sin(final_phi)]])
P = w / w_tot
pos = np.r_[pos, new]
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)]])
w_idx = np.where(phi == final_phi)
print('w_idx = ', w_idx[0])
w_jl = w_jl*w[w_idx[0]]
w_l = w_l * w_tot
return pos
\ No newline at end of file
return pos, w_jl, w_l
\ No newline at end of file
......@@ -2,21 +2,23 @@ import numpy as np
from matplotlib import pyplot as plt
from functions import add_bead
# Simultion constants
# Simulation constants
# Initialization
n = 6 # number of phis
T_dim = 0.001
pos = np.array([[0, 0], [1, 0]])
n_beads = 30
n_beads = 100
T_dim = 1
pos = np.zeros([n_beads, 2])
pos[1] = np.array([1, 0])
w_jl = 1
w_l = 1
for j in range(2,n_beads):
for j in range(n_beads):
pos = add_bead(pos, n, T_dim)
pos, w_jl, w_l = add_bead(pos, n, T_dim, j, w_jl, w_l)
P = w_jl/w_l
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