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

made some functions more efficient, added errorbars, changed/added doc strings to various functions

parent b962629c
import numpy as np
from parameters import epsilon, sigma, kbT, theta, k, N
def LJpot(roptions, i, xy):
def LJpot(roptions, index, xy):
"""
Lennard Jones potential function
Inputs:
roptions - potential positions for next bead
index - bead number in the chain
xy - coordinates of all beads
Lennard Jones potential function. Calculates total energy of given
positions.
"""
n = len(roptions)
Vjtot = np.zeros((n, 1))
for j in range(n):
deltaxyj = xy[np.arange(i)] - roptions[j]
r = np.sqrt(np.sum(deltaxyj**2, axis = 1))
Vj = 4*epsilon*((sigma/r)**12 - (sigma/r)**6)
Vjtot = np.zeros(len(roptions))
for j in range(len(roptions)):
deltaxyj = xy[0:index] - roptions[j]
r2 = np.sum(deltaxyj**2, axis = 1)
Vj = 4*epsilon*((sigma**2/r2)**6 - (sigma**2/r2)**3)
Vjtot[j] = np.nansum(Vj)
return Vjtot
def MCweight(Ej):
"""
Attributing Monte Carlo weights to positions.
Inputs:
Ej - Lennard Jones potential energies of positions j
Attributes Monte Carlo weights to positions. Rescales energies if exponent
gets too large. Rolls a random number and assigns it into one of the bins
which correspond to the Boltzmann factor of the energies Ej.
"""
if max(Ej/kbT > 100):
Ej *= 10**(-np.log10(max(Ej/kbT)/100))
......@@ -31,16 +41,21 @@ def MCweight(Ej):
def neighbourhood(index, xy, roptions):
"""
For the last bead in the chain it is checked with other beads lie within
a circle of radius sqrt(1.25). 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.
Inputs:
index - bead number in the chain
xy - coordinates of all beads
roptions - potential positions for next bead
For the last bead in the chain it is checked which other beads lie within
or on a circle of radius sqrt(1.25). 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(xy2)) if xy2[j] < 1.25]
close = [j for j in range(len(xy2)) if xy2[j] <= 1.25]
close.remove(index-1) #remove the second to last point,
#since you cannot intersect this line
valid = list(range(len(roptions)))
......@@ -61,6 +76,8 @@ def neighbourhood(index, xy, roptions):
def direction(a, b, c):
"""
Inputs:
- points a, b, c
Checks what the orienation of 3 given points is.
Output correspond to colinear, anti-clockwise or clockwise.
"""
......@@ -72,8 +89,10 @@ def direction(a, b, c):
else:
return 3 #a, b, c are clockwise
def checker(a,b,c,d):
def checker(a, b, c, d):
"""
Inputs:
- points 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
......@@ -89,12 +108,7 @@ def checker(a,b,c,d):
else:
return False
def end2end(xy):
"""
Calculates the distance between the first and the last bead."
"""
length = np.sqrt(np.sum((xy[-1]-xy[0])**2))
return length
def chainmaker(chains = None, k = None, flipping = True):
"""
......@@ -139,7 +153,7 @@ def chainmaker(chains = None, k = None, flipping = True):
xy[i] = roptions[j]
i += 1
beads = len(xy)
length = end2end(xy)
length = np.sqrt(np.sum((xy[-1]-xy[0])**2))
if chains is not None:
chains[k, :len(xy), :] = xy
......@@ -155,24 +169,32 @@ def data(storechains = True, flipping = True):
beads, lengths = zip(*[chainmaker(flipping) for i in range(k)])
return list(beads), list(lengths)
def e2edis(chains, beadstotal):
N = chains.shape[0]
beads = []
lengths = []
for i in range(N):
for j in range(1, beadstotal[i]):
chain = chains[i, 0:j]
r0 = chain[0]
rmax = chain[-1]
R2 = np.sum((rmax - r0)**2)
beads.append(j)
lengths.append(R2)
return beads, lengths
def e2edis(chains):
"""
Inputs:
- A set of polymer chains with size N
- The total number of beads in each chain
For a given set of chains, this function calculates the end to end length
of the subchains within a chain from point 0 to j. These lengths are added
to a single list and
"""
chains = np.where(chains!=0, chains, np.nan)
chains0 = np.repeat(chains[:,0,:],chains.shape[1], axis=0).reshape(chains.shape)
R2 = np.sum((chains - chains0)**2, axis=2)
R2std = np.nanstd(R2, axis=0)
R2std = R2std[~np.isnan(R2std)]
R2mean = np.nanmean(R2, axis=0)
R2mean = R2mean[~np.isnan(R2mean)]
return R2mean[2:], R2std[2:]
def linefit(beads, lengths):
beads = np.array(beads)
lengths = np.array(lengths)
a = np.sum(lengths*(beads-1)**1.5)/np.sum((beads-1)**3)
def linefit(beads, R2mean):
"""
Calculates parameter a for the function y = a*(N-1)^1.5. Takes inputs of
number of beads N and corresponding end to end lengths y.
"""
a = np.sum(R2mean*(beads-1)**1.5, )/np.sum((beads-1)**3, dtype=float)
return a
......@@ -6,6 +6,5 @@ plt.close("all")
beads, lengths, chains = functions.data(storechains = True, flipping = True)
beadstot, lengthstot = functions.e2edis(chains, beads)
plots.lengthplot(beadstot, lengthstot)
plots.lengthplot(chains)
#plots.chainplot(beads, chains, np.arange(10))
\ No newline at end of file
import numpy as np
#Physical constants
kbT = 1e+12
kbT = 1e+2
epsilon = 0.25
sigma = 0.8
......@@ -9,7 +9,7 @@ n = 6
theta = 2*np.pi/n*np.arange(n)
#Max chain length
N = 1000
N = 256
#iterations
k = 1000
\ No newline at end of file
......@@ -3,26 +3,27 @@ import matplotlib.pyplot as plt
import functions
from parameters import kbT
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())]
a = functions.linefit(beads, lengths)
N = np.arange(1, max(beads))
def lengthplot(chains):
# 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())]
R2mean, R2std = functions.e2edis(chains)
N = np.arange(len(R2mean))+2
a = functions.linefit(N, R2mean)
y = a*(N-1)**1.5
print("a=%.2f"%a)
exp = int(np.log10(kbT))
plt.figure()
plt.title("kbT = 10^%d" % exp)
plt.loglog(N, y, c = 'r', linestyle = '--')
plt.xlim((1, 2*max(beads)))
plt.ylim((0.5, 2*max(lengths)))
plt.scatter(beads, lengths)
plt.xlim((1, 2*max(N)))
plt.ylim((0.5, 2*max(R2mean)))
plt.errorbar(N, R2mean, R2std)
plt.ylabel("<R²>")
plt.xlabel("Polymer beads")
plt.show()
......
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