finalfunctions2.py 5.31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
############################################### 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)
Laurien's avatar
Laurien committed
40
    E_theta = 4*eps * np.sum((radius/sigma)**-12 - (radius/sigma)**-6, axis=0)
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
    
    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,:]
    
Laurien's avatar
Laurien committed
58
    return w_lfinal, sumW_l, position, n
59 60


Laurien's avatar
Laurien committed
61
def upperlowerlimit(i, k_b, T, alpha1, alpha2, w_lfinal, w_ave, lowlimit, uplimit):
62 63 64 65 66
    """
    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
    """
Laurien's avatar
Laurien committed
67
    w_ave=np.average(w_lfinal[:i])
68 69 70 71 72 73 74
    weight3=w_lfinal[2]
    lowlimit=alpha1*(w_ave/weight3)
    uplimit=alpha2*(w_ave/weight3)    
    
    return lowlimit, uplimit

def crossing(i, position):
Laurien's avatar
Laurien committed
75 76 77 78 79 80 81 82
    """
    **** Only used for checks, not implemented in the algorithm ****
    
    function: Makes sure the polymer beads are not allowed to come too close/overlap
    input: positions
    output: prints "crossing!" if polymer beads come too close
    """
    
83 84 85 86 87 88 89 90 91 92 93 94 95
    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):
Laurien's avatar
Laurien committed
96 97 98 99 100
    """
    function: calculates the end-to-end length of a polymer 
    input: maximum bead 
    output: R
    """
101 102 103 104 105 106
    e2e = position[maxbead-1,:] - position[0,:]
    e2e = np.linalg.norm(e2e) 
    
    return e2e

def calc_mean(obs, weight):
Laurien's avatar
Laurien committed
107 108 109 110 111
    """
    function: calculates the mean of an observable
    input: observable, weights
    output: mean value
    """
Laurien's avatar
Laurien committed
112
    weightedobs = np.multiply(obs, weight)
113 114 115 116 117
    meanvalue = np.sum(weightedobs, axis=1) /  np.sum(weight, axis=1)
    
    return meanvalue

def calc_std(obs, weight):
Laurien's avatar
Laurien committed
118 119 120 121 122
    """
    function: calculates the standard deviation of an observable
    input: observable, weights
    output: standard deviation
    """
123 124 125 126 127 128
    variance = calc_mean(obs**2, weight) - calc_mean(obs, weight)**2
    std = np.sqrt(variance)
    
    return std

def function(x,a,b):
Laurien's avatar
Laurien committed
129 130 131 132 133 134
    """
    function: theoretical function a * x^b to perform the least-square fit on
    input: x data, parameters a and b
    output: function y = a * x^b
    """
    return a * (x-1)**b
Laurien's avatar
Laurien committed
135

Laurien's avatar
Laurien committed
136 137 138 139 140 141 142 143
def enrich(PolWeight, NewWeight, position):
    """
    function: performs enriching of the polymer ensemble (PERM)
    input: positions, polymer weight
    output: New weight, positions duplicate
    """
    NewWeight1 = 0.5 * PolWeight
    newchain = position
Laurien's avatar
Laurien committed
144
    
Laurien's avatar
Laurien committed
145
    NewWeight = NewWeight1 
Laurien's avatar
Laurien committed
146
    
Laurien's avatar
Laurien committed
147
    return NewWeight, newchain
Laurien's avatar
Laurien committed
148

Laurien's avatar
Laurien committed
149 150 151 152 153 154
def prune(PolWeight, NewWeight, stopcondition):
    """
    function: performs pruning of the polymer ensemble (PERM)
    input: polymer weight
    output: New weight, stop condition
    """
Laurien's avatar
Laurien committed
155 156 157 158
    R=np.random.rand(1)
    
    if(R<0.5): 
        NewWeight = 2 * PolWeight
Laurien's avatar
Laurien committed
159 160 161
    else:
        stopcondition = True
    
Laurien's avatar
Laurien committed
162
            
Laurien's avatar
Laurien committed
163 164 165 166 167 168 169 170 171 172 173 174 175
    return NewWeight, stopcondition

def fitting(x, obs, min, max):
    """
    function: performs a least-squares fit on the data points, with respect to the theoretical exponential function
    Input: R
    Output: optimum values and (co)variance
    """   
    x = np.linspace(2,max-min+2,max-3)
    opt, cov = curve_fit(function, x, obs[min:max])
    perr = np.sqrt(np.diag(cov))
    
    return opt, cov, perr
Laurien's avatar
Laurien committed
176 177 178 179
    
    
    
    
180 181 182 183 184