Commit 4beda16e authored by Luuk Balkenende's avatar Luuk Balkenende

remove temporary files

parent 1e299ea1
# Project2
Computational physics project 2: Monte Carlo methods
\ No newline at end of file
import numpy as np
import sys
def add_bead(pos, n, T_dim, j, p_pol):
"""Adds one bead to polymer, according to the Rosenblaum method
Parameters:
----------
pos: array of size (n_beads, 2)
Position of n_beads beads in x,y direction
n: integer > 0
Number of possible angles
T_dim: float > 0
Non-dimensionalized temperature
j: integer > 1
Bead j is the j-th bead which is added to polymer
p_pol: GAAT WAARSCHIJNLIJK TOCH NOG VERANDEREN
Results:
--------
pos: array of size (n_beads, 2)
Position of n_beads beads in x,y direction
p_pol: GAAT WAARSCHIJNLIJK TOCH NOG VERANDEREN
"""
# 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
# Calculate new positions
pos_new = np.zeros([n, 2])
pos_new[:, 0] = pos[j - 1, 0] + np.cos(phi)
pos_new[:, 1] = pos[j - 1, 1] + np.sin(phi)
# Calculate distance between new bead and beads
rx = np.zeros([n, j])
ry = np.zeros([n, j])
for i in range(n):
rx[i, :] = pos_new[i, 0] - pos[0:j, 0]
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 == 0:
print('Help, I have nowhere to go')
print('w_tot = ', w_tot)
sys.exit()
# Calculate new position of next bead
P = w / w_tot
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)]])
# 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]]
w_l = w_tot
p_pol = p_pol * w_jl / w_l
return pos, p_pol
\ No newline at end of file
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
import time
from functions import add_bead
# Simulation parameters
n_beads = 100 #number of beads
n = 6 #number of angles
n_pol = 2 #number of polymers
# Physical parameters
T_dim = 10 #T * (Kb / eps) # Non-dimensional temperature
# Preallocate arrays
pos = np.zeros([n_beads, 2])
p = np.zeros(n_pol)
# Initialization
pos[1] = np.array([1, 0]) #bead 1: pos[0, 0]. bead 2: pos[1, 0]
# Make n polymers
for k in range(n_pol):
# 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):
pos, p_pol = add_bead(pos, n, T_dim, j, p_pol)
p[k] = p_pol #Weight per polymer
print('P = ', p)
# Plot polymer
plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
plt.show()
import numpy as np
def add_bead(pos, n, T_dim, L, polweight):
# 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
# Calculate new positions
pos_new = np.zeros([n, 2])
pos_new[:, 0] = pos[L - 1, 0] + np.cos(phi)
pos_new[:, 1] = pos[L - 1, 1] + np.sin(phi)
# Calculate distance between new bead and beads
rx = np.zeros([n, L])
ry = np.zeros([n, L])
for i in range(n):
rx[i, :] = pos_new[i, 0] - pos[0:L, 0]
ry[i, :] = pos_new[i, 1] - pos[0:L, 1]
r = np.sqrt(rx ** 2 + ry ** 2)
r[r == 0] = np.nan
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)
print(w_l)
#print(w_l)
polweight = polweight*w_l #weight per polymer off by our algorithm
# Calculate new position of next bead
P = w_jl / w_l
final_phi = np.random.choice(phi, p=P)
pos[L] = np.array([[pos[L - 1, 0] + np.cos(final_phi), pos[L - 1, 1] + np.sin(final_phi)]])
return pos, polweight
\ No newline at end of file
# Import libaries/functions
import numpy as np
from matplotlib import pyplot as plt
from new_functions import add_bead
# Simulation parameters
n_beads = 250 # number of beads
n = 6 # number of angles
n_pol = 10 # number of polymers
# Physical parameters
T_dim = 0.1 # T * (Kb / eps) # Non-dimensional temperature
polweight = 1
# Make n polymers
k = 0
while n_pol > k:
# 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]
#AddBead(Polymer,Polweight,L)
#Add third bead
L = 2
while L < n_beads:
#add next bead
pos, polweight = add_bead(pos, n, T_dim, L, polweight)
L = L +1
k = k+1
print(k)
plt.figure()
plt.plot(pos[:,0], pos[:,1], marker= ".", markersize=20, color="k" )
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