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

improved end to end distance function to include lengths of subchains, added a...

improved end to end distance function to include lengths of subchains, added a least squares approximated fit line
parent 581f9ee2
......@@ -96,7 +96,7 @@ def end2end(xy):
length = np.sqrt(np.sum((xy[-1]-xy[0])**2))
return length
def chainmaker(chains = None, k = None):
def chainmaker(chains = None, k = None, flipping = True):
"""
Function used to produce k chains with a maximum length N. Calculates the
position to place the next bead based on weights determined by Lennard-Jones
......@@ -120,12 +120,16 @@ def chainmaker(chains = None, k = None):
valid = neighbourhood(i-1, xy, roptions)
roptions = roptions[valid]
if len(roptions) == 0:
if flipped == True:
if flipping == True:
if flipped == True:
xy=xy[:i]
break
flipped = True
xy[0:i] = np.flip(xy[0:i])
continue
else:
xy=xy[:i]
break
flipped = True
xy[0:i] = np.flip(xy[0:i])
continue
Ej = LJpot(roptions, i, xy)
j = MCweight(Ej)
......@@ -141,12 +145,35 @@ def chainmaker(chains = None, k = None):
chains[k, :len(xy), :] = xy
return beads, length
def data(storechains = True):
def data(storechains = True, flipping = True):
if storechains == True:
global chains
chains = np.zeros((k, N, 2))
beads, lengths = zip(*[chainmaker(chains, i) for i in range(k)])
beads, lengths = zip(*[chainmaker(chains, i, flipping) 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
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 linefit(beads, lengths):
beads = np.array(beads)
lengths = np.array(lengths)
a = np.sum(lengths*(beads-1)**1.5)/np.sum((beads-1)**3)
return a
\ No newline at end of file
......@@ -4,7 +4,8 @@ import functions
import plots
plt.close("all")
beads, lengths, chains = functions.data()
beads, lengths, chains = functions.data(storechains = True, flipping = True)
plots.lengthplot(beads, lengths)
plots.chainplot(beads, chains, np.arange(10))
\ No newline at end of file
beadstot, lengthstot = functions.e2edis(chains, beads)
plots.lengthplot(beadstot, lengthstot)
#plots.chainplot(beads, chains, np.arange(10))
\ No newline at end of file
import numpy as np
#Physical constants
kbT = 1e+1
kbT = 1e+12
epsilon = 0.25
sigma = 0.8
......@@ -12,4 +12,4 @@ theta = 2*np.pi/n*np.arange(n)
N = 1000
#iterations
k = 10
\ No newline at end of file
k = 1000
\ No newline at end of file
import numpy as np
import matplotlib.pyplot as plt
import functions
from parameters import kbT
def lengthplot(beads, lengths):
beadlengths = {}
for i in range(len(beads)):
......@@ -10,17 +12,24 @@ def lengthplot(beads, lengths):
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))
y = a*(N-1)**1.5
exp = int(np.log10(kbT))
plt.figure()
plt.title("kbT = %.2f" % (kbT))
plt.loglog()
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.ylabel("End to end length")
plt.ylabel("<R²>")
plt.xlabel("Polymer beads")
plt.show()
def chainplot(beads, chains, interval):
if interval == None:
interval = range(chains.shape[0])
# if interval == None:
# interval = range(chains.shape[0])
for i in interval:
plt.figure()
plt.title("Chain %d" % i)
......
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