Commit a9d0f582 authored by Gilliam van Oudenaren's avatar Gilliam van Oudenaren
Browse files

Added some docstrings, added some functions to store the generated chains and...

Added some docstrings, added some functions to store the generated chains and created a seperate file to plot results of chain generation.
parent 4a204e7d
import numpy as np
from parameters import epsilon, sigma, kbT
from parameters import epsilon, sigma, kbT, theta, k, N
def LJpot(roptions, i, xy):
"""
Lennard Jones potential function
......@@ -29,10 +30,19 @@ def MCweight(Ej):
return j
def radiussqrt2(index, xy, roptions):
"""
For the last bead in the chain it is checked with other beads lie within
a circle of radius sqrt(2). For each of these points, the checker function
is applied to check if the lines between adjacent points will intersect
the lines of the new possible positions. If this is the case, these
positions are removed from the set of valid potential points.
"""
xydis = xy[:index]-xy[index]
xy2 = np.sum(xydis**2, axis = 1)
close = [j for j in range(len(xydis)) if xy2[j] < 2]
close.remove(index-1)
close = [j for j in range(len(xy2)) if xy2[j] < 2]
close.remove(index-1) #remove the second to last point,
#since you cannot intersect this line
valid = list(range(len(roptions)))
for i in close:
for j in range(len(roptions)):
......@@ -50,15 +60,25 @@ def radiussqrt2(index, xy, roptions):
return valid
def direction(a, b, c):
"""
Checks what the orienation of 3 given points is.
Output correspond to colinear, anti-clockwise or clockwise.
"""
val = (b[1]-a[1])*(c[0]-b[0])-(b[0]-a[0])*(c[1]-b[1])
if val == 0:
return 1
return 1 #a, b, c are colinear
elif val < 0:
return 2
return 2 #a, b, c are anti-clockwise
else:
return 3
return 3 #a, b, c are clockwise
def checker(a,b,c,d):
"""
Checks if the linesegments between a and b, and c and d intersect. Calls
the direction function to check for orientation. Operates on the assumption
that starting/end points of the line segments are never exactly on the
other line segment, which
"""
dir1 = direction(a, b, c)
dir2 = direction(a, b, d)
dir3 = direction(c, d, a)
......@@ -67,4 +87,55 @@ def checker(a,b,c,d):
if dir1 is not dir2 and dir3 is not dir4:
return True
else:
return False
\ No newline at end of file
return False
def end2end(xy):
length = np.sqrt(np.sum((xy[-1]-xy[0])**2))
return length
def chainmaker(chains = None, k = None):
xy = np.zeros((N, 2))
xy[0] = [0, 0]
xy[1] = [1, 0]
i = 2
flipped = False
while i < len(xy):
theta0 = 2*np.pi*np.random.rand()
r0 = xy[i-1]
roptions = r0 + np.transpose(np.array([np.sin(theta+theta0), np.cos(theta+theta0)]))
if i > 2:
valid = radiussqrt2(i-1, xy, roptions)
roptions = roptions[valid]
if len(roptions) == 0:
if flipped == True:
xy=xy[:i]
break
flipped = True
xy[0:i] = np.flip(xy[0:i])
continue
Ej = LJpot(roptions, i, xy)
j = MCweight(Ej)
if type(j) == bool:
xy=xy[:i]
break
xy[i] = roptions[j]
i += 1
beads = len(xy)
length = end2end(xy)
if chains is not None:
chains[k, :len(xy), :] = xy
return beads, length
def data(storechains = True):
if storechains == True:
global chains
chains = np.zeros((k, N, 2))
beads, lengths = zip(*[chainmaker(chains, i) for i in range(k)])
return list(beads), list(lengths), chains
else:
beads, lengths = zip(*[chainmaker() for i in range(k)])
return list(beads), list(lengths)
\ No newline at end of file
import numpy as np
import matplotlib.pyplot as plt
import functions
from parameters import N, theta
import plots
plt.close("all")
plt.close("all")
beads, lengths, chains = functions.data()
xy = np.zeros((N, 2))
xy[0] = [0, 0]
xy[1] = [1, 0]
i = 2
while i < len(xy):
print("step %d" % i)
theta0 = 2*np.pi*np.random.rand()
r0 = xy[i-1]
roptions = r0 + np.transpose(np.array([np.sin(theta+theta0), np.cos(theta+theta0)]))
if i > 2:
valid = functions.radiussqrt2(i-1, xy, roptions)
roptions = roptions[valid]
if len(roptions) == 0:
print('intersection error')
xy=xy[:i]
break
Ej = functions.LJpot(roptions, i, xy)
j = functions.MCweight(Ej)
if type(j) == bool:
print('weight error')
xy=xy[:i]
break
xy[i] = roptions[j]
i += 1
txtstr = [str(i) for i in range(len(xy))]
plt.figure()
plt.title("Final chain")
plt.plot(xy[:, 0], xy[:, 1])
plt.scatter(xy[:, 0], xy[:, 1],color = 'k')
for i in range(len(xy)):
plt.text(xy[i, 0], xy[i, 1],txtstr[i])
plots.lengthplot(beads, lengths)
plots.chainplot(beads, chains)
\ No newline at end of file
......@@ -8,5 +8,8 @@ sigma = 0.8
n = 6
theta = 2*np.pi/n*np.arange(n)
#Chain length
N = 10000
#Max chain length
N = 1000
#iterations
k = 10
\ No newline at end of file
import numpy as np
import matplotlib.pyplot as plt
def lengthplot(beads, lengths):
beadlengths = {}
for i in range(len(beads)):
if beads[i] in beadlengths.keys():
beadlengths[beads[i]].append(lengths[i])
else:
beadlengths[beads[i]] = [lengths[i]]
beads = list(beadlengths.keys())
lengths = [sum(i)/len(i) for i in list(beadlengths.values())]
plt.figure()
plt.loglog()
plt.scatter(beads, lengths)
plt.ylabel("End to end length")
plt.xlabel("Polymer beads")
plt.show()
def chainplot(beads, chains):
for i in range(chains.shape[0]):
plt.figure()
plt.title("Chain %d" % i)
plt.plot(chains[i, :beads[i], 0], chains[i, :beads[i], 1])
plt.scatter(chains[i, :beads[i], 0], chains[i, :beads[i], 1],color = 'k')
txtstr = [str(i+1) for i in range(beads[i])]
for j in range(beads[i]):
plt.text(chains[i, j, 0], chains[i, j, 1], txtstr[j])
Supports Markdown
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