diff --git a/finalfunctions2.py b/finalfunctions2.py index 9b7e5b3e340b03cbdacfbec7f8edbe01d05b55e3..2690006be54b17e8a31abe3233e40ad4c0b243c6 100644 --- a/finalfunctions2.py +++ b/finalfunctions2.py @@ -37,7 +37,7 @@ def lennardjones(eps, sigma, Ntheta, pos_position, position): """ 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) + E_theta = 4*eps * np.sum((radius/sigma)**-12 - (radius/sigma)**-6, axis=0) return E_theta @@ -55,16 +55,16 @@ def extrabead(i, k_b, T, Ntheta, E_theta, sumW_l, w_lfinal, position, pos_positi w_lfinal = w_l[n] position = pos_position[n,:] - return w_lfinal, sumW_l, position + return w_lfinal, sumW_l, position, n -def upperlowerlimit(alpha1, alpha2, w_lfinal, w_ave, lowlimit, uplimit): +def upperlowerlimit(i, k_b, T, 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) + w_ave=np.average(w_lfinal[:i]) weight3=w_lfinal[2] lowlimit=alpha1*(w_ave/weight3) uplimit=alpha2*(w_ave/weight3) @@ -93,43 +93,86 @@ def crossing(i, position): return cross def endtoend(position, maxbead): + """ + function: calculates the end-to-end length of a polymer + input: maximum bead + output: R + """ e2e = position[maxbead-1,:] - position[0,:] e2e = np.linalg.norm(e2e) return e2e def calc_mean(obs, weight): + """ + function: calculates the mean of an observable + input: observable, weights + output: mean value + """ weightedobs = np.multiply(obs, weight) meanvalue = np.sum(weightedobs, axis=1) / np.sum(weight, axis=1) return meanvalue def calc_std(obs, weight): - + """ + function: calculates the standard deviation of an observable + input: observable, weights + output: standard deviation + """ 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 + """ + 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 -def enrich(PolWeight, NewWeight): +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 - NewWeight = 0.5 * PolWeight - #print("Weight is halved!") + NewWeight = NewWeight1 - return NewWeight + return NewWeight, newchain -def prune(PolWeight, NewWeight): - +def prune(PolWeight, NewWeight, stopcondition): + """ + function: performs pruning of the polymer ensemble (PERM) + input: polymer weight + output: New weight, stop condition + """ R=np.random.rand(1) if(R<0.5): NewWeight = 2 * PolWeight + else: + stopcondition = True + - return NewWeight + 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