Commit 8a07999f authored by Laurien's avatar Laurien

Few adjustments regarding error calculation and plots

parent 7dde6054
This source diff could not be displayed because it is too large. You can view the blob instead.
############################################### IMPORT MODULES ###############################################
import numpy as np
np.seterr(divide='ignore', invalid='ignore')
import matplotlib.pyplot as plt
#from anim import make_3d_animation
import scipy as sci
from scipy import integrate
from scipy import optimize
from scipy.optimize import curve_fit
import math
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import cnames
from matplotlib import animation
from IPython.display import HTML
################################# FUNCTIONS #############################################
def roulette(i, d, Ntheta, position, theta):
"""
function: determines the possible new positions of the particle, based on a distribution of theta
input: previous position, theta
output: new position, distance between the two particles(x,y), absolute distance between the particles (radius)
"""
alpha = 2*np.pi * np.random.rand(1)
possiblestep = np.array([d*np.cos(theta+alpha), d*np.sin(theta+alpha)]).T
pos_position = np.tile(position[i-1, :],(Ntheta,1)) + possiblestep
return pos_position
def lennardjones(eps, sigma, Ntheta, pos_position, position):
"""
function: calculate the energy needed to put an extra bead in each of the possible positions
input: possible positions
output: E_theta
"""
positionmatrix = position[:, np.newaxis, :]
radius = np.linalg.norm(np.repeat(positionmatrix, Ntheta, axis=1) - pos_position + 0.00001, axis=2)
E_theta = 4*eps * np.sum((radius/sigma)**-12 - (radius/sigma)**-6, axis=0)
return E_theta
def extrabead(i, k_b, T, Ntheta, E_theta, sumW_l, w_lfinal, position, pos_position):
"""
function: calculate weights corresponding to the functions
input: potential energy, zero matrices
output: weight corresponding to the chosen position, sum of all weights, chosen position of theta: n
"""
w_l = np.exp(-E_theta/(k_b*T))
sumW_l = np.sum(w_l)
probability = w_l/sumW_l
choices = np.arange(Ntheta)
n = int(np.round(np.random.choice(a=choices, p=probability)))
w_lfinal = w_l[n]
position = pos_position[n,:]
return w_lfinal, sumW_l, position
def upperlowerlimit(alpha1, alpha2, w_lfinal, w_ave, lowlimit, uplimit):
"""
function: calculate weights corresponding to the functions
input: potential energy, zero matrices
output: weight corresponding to the chosen position, sum of all weights, chosen position of theta: n
"""
w_ave=np.average(w_lfinal)
weight3=w_lfinal[2]
lowlimit=alpha1*(w_ave/weight3)
uplimit=alpha2*(w_ave/weight3)
return lowlimit, uplimit
def crossing(i, position):
cross = False
dx = np.zeros(shape=(L, L-1))
for j in range(0,i):
dx[i,j] = np.linalg.norm(position[i,:] - position[j,:], axis=0)
check = dx[i, 0:i]
if np.min(check) < treshold:
cross = True
return cross
def endtoend(position, maxbead):
e2e = position[maxbead-1,:] - position[0,:]
e2e = np.linalg.norm(e2e)
return e2e
def calc_mean(obs, weight):
weightedobs = obs * weight
meanvalue = np.sum(weightedobs, axis=1) / np.sum(weight, axis=1)
return meanvalue
def calc_std(obs, weight):
variance = calc_mean(obs**2, weight) - calc_mean(obs, weight)**2
std = np.sqrt(variance)
return std
def function(x,a,b):
return a *(x)**(b)
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