finalfunctions2.py 5.31 KB
 Laurien committed May 07, 2019 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 committed May 15, 2019 40 `````` E_theta = 4*eps * np.sum((radius/sigma)**-12 - (radius/sigma)**-6, axis=0) `````` Laurien committed May 07, 2019 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 committed May 15, 2019 58 `````` return w_lfinal, sumW_l, position, n `````` Laurien committed May 07, 2019 59 60 `````` `````` Laurien committed May 15, 2019 61 ``````def upperlowerlimit(i, k_b, T, alpha1, alpha2, w_lfinal, w_ave, lowlimit, uplimit): `````` Laurien committed May 07, 2019 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 committed May 15, 2019 67 `````` w_ave=np.average(w_lfinal[:i]) `````` Laurien committed May 07, 2019 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 committed May 08, 2019 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 """ `````` Laurien committed May 07, 2019 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 committed May 15, 2019 96 97 98 99 100 `````` """ function: calculates the end-to-end length of a polymer input: maximum bead output: R """ `````` Laurien committed May 07, 2019 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 committed May 15, 2019 107 108 109 110 111 `````` """ function: calculates the mean of an observable input: observable, weights output: mean value """ `````` Laurien committed May 08, 2019 112 `````` weightedobs = np.multiply(obs, weight) `````` Laurien committed May 07, 2019 113 114 115 116 117 `````` meanvalue = np.sum(weightedobs, axis=1) / np.sum(weight, axis=1) return meanvalue def calc_std(obs, weight): `````` Laurien committed May 15, 2019 118 119 120 121 122 `````` """ function: calculates the standard deviation of an observable input: observable, weights output: standard deviation """ `````` Laurien committed May 07, 2019 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 committed May 15, 2019 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 committed May 08, 2019 135 `````` `````` Laurien committed May 15, 2019 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 committed May 08, 2019 144 `````` `````` Laurien committed May 15, 2019 145 `````` NewWeight = NewWeight1 `````` Laurien committed May 08, 2019 146 `````` `````` Laurien committed May 15, 2019 147 `````` return NewWeight, newchain `````` Laurien committed May 08, 2019 148 `````` `````` Laurien committed May 15, 2019 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 committed May 08, 2019 155 156 157 158 `````` R=np.random.rand(1) if(R<0.5): NewWeight = 2 * PolWeight `````` Laurien committed May 15, 2019 159 160 161 `````` else: stopcondition = True `````` Laurien committed May 08, 2019 162 `````` `````` Laurien committed May 15, 2019 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 committed May 08, 2019 176 177 178 179 `````` `````` Laurien committed May 07, 2019 180 181 182 183 184 `````` ``````