Commit 3a520a3a authored by Pieter van Velde's avatar Pieter van Velde

Implemented fit function in proccesfunctions.py

parent ce176b3a
# Project2
Computational physics project 2: Monte Carlo methods
\ No newline at end of file
Computational physics project 2: Monte Carlo simulations of polymers interacting according the Lennard-Jones potential
Made by:
Luuk Balkenende
Pieter van Velde
15-5-2019
\ No newline at end of file
......@@ -3,14 +3,14 @@ from simfunctions import calc_pol, allocate, init, perm, add_bead, limits, prune
from processfunctions import exp_calc, bootstrap
# Simulation parameters
n_beads = 250 # number of beads
n_beads = 100 # number of beads
n = 6 # number of angles
n_pol = 300 # number of polymers
n_pol = 50 # number of polymers
alfa_up = 2 # alfa to compute upper limit
alfa_low = 0.7 # alfa to compute lower limit
# Choose algorithm
PERM = True # Put on True if PERM algorithm is desired
PERM = False # Put on True if PERM algorithm is desired
random_walk = False # Put on True if random walk is desired
# Physical parameters
......
import numpy as np
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline
def exp_calc(polweight_array, end_end_array, n_beads):
"""Calculates the end to end distances with taking into account their weights, and
......@@ -19,39 +21,43 @@ def exp_calc(polweight_array, end_end_array, n_beads):
Figure: end to end distance against number of beads
"""
end_end_times_weight = np.multiply(end_end_array, polweight_array)
sum_end_end_times_weight = np.sum(end_end_times_weight, axis=1)
sum_weights = np.sum(polweight_array, axis=1)
exp_value = np.divide(sum_end_end_times_weight, sum_weights)
#foute manier van berekenen std, moet ook sws met kwadraten ofzoiets
std_weight_end = np.std(end_end_times_weight, axis=1)
std_weight = np.std(polweight_array, axis=1)
exp_value_std2 = np.divide(std_weight_end, std_weight)
exp_value_mean, exp_value_std = bootstrap(polweight_array, end_end_array, n_beads)
#print('bootstrap = ', exp_value_std)
#print('normal std = ', exp_value_std2)
#print('bootstrap = ', exp_value_std.shape)
#print('normal std = ', exp_value_std2.shape)
a = np.arange(0, n_beads - 1, 1)
b = np.arange(0, n_beads, 1)
exp_value[0:2] = [0, 1]
exp_value_mean[0:2] = [0, 1]
font = {'family' : 'DejaVu Sans',
'weight' : 'normal',
'size' : 17}
matplotlib.rc('font', **font)
q = np.polyfit(np.log(a[1::]), np.log(exp_value_mean[1::]**2), 1)
print(q)
xx = np.linspace(0, 250 , 1000)
yy = np.exp(q[1])*xx**q[0]
plt.figure(1)
plt.title('bootstrap')
plt.xscale("log")
plt.yscale("log")
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, xerr=0, yerr=2*exp_value_std[:]**2, color="k")
plt.ylabel('R^2 (end to end distance)')
plt.xlabel('N (Number of beads)')
plt.grid(True, which="both", ls="-")
plt.ylim(1, 10000)
plt.plot( xx, yy, 'r--', label="Fit", linewidth =3)
plt.xlabel(r'$N_{beads}$')
plt.ylabel(r'$R^2$ $[\sigma^2]$')
plt.xscale("log")
plt.yscale("log")
plt.errorbar(a, exp_value_mean[0:n_beads - 1] ** 2, fmt='x', xerr=None, yerr=exp_value_std[:]**2, label='Data',color="k", capsize=3, capthick=1, markersize=3)
plt.grid(True, which="both", ls=":")
plt.legend(loc='best')
plt.xlim(1, 250)
plt.ylim(bottom=1)
plt.ylim(top=10000)
plt.savefig('joe', bbox_inches='tight', dpi=300)
plt.figure(3)
plt.title('Population size')
......
import numpy as np
import random
def calc_pol(n_beads, n, n_pol, alfa_low, alfa_up, T_dim, PERM, random_walk):
"""Calculates configuration of n_pol polymers up to a length of n_beads
......
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