First Commit

Monte Carlo - Metropolis/autocorrelation.py

@@ -0,0 +1,29 @@
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+
+def autocorrelation(data,tenmperature,k):
+#with open("/home/abhinav/Desktop/Academics/Summer/Python/data.txt","r") as f:
+#	data = np.loadtxt(f,delimiter="\n")
+
+	n = len(data)
+
+	r = [0.0]*n
+
+	mean = np.mean(data,dtype=np.float64)
+	var = np.var(data,dtype=np.float64)
+	t = n*var
+
+	for i in range(n):
+		for j in range(n-i):
+			r[i] += (data[j] - mean)*(data[i+j] - mean)/t# ((n - i)*var		)
+
+	x = range(n)
+
+	plt.figure()
+	plt.plot(x,r)
+	plt.title(tenmperature)
+	plt.savefig('auto {0}.jpg'.format(k))
+	#plt.show()
+
+	return r

Monte Carlo - Metropolis/benchmark_ising.py

@@ -0,0 +1,81 @@
+import math
+import numpy as np
+import matplotlib.pyplot as plt
+
+from system_energy import *
+from spinenergy import *
+from magnet import *
+
+def initialize_b(m,n,p):
+	if p == 1:
+		sx = 0.0*np.random.rand(m,n)
+		sy = 0.0*np.random.rand(m,n) + 1.0
+
+		return sx,sy
+
+	if p > 1:
+		sx = 0.0*np.random.rand(m,n)
+		sy = 0.0*np.random.rand(m,n) + 1.0
+		sz = 0.0*np.random.rand(m,n) 
+
+		return sx,sy,sz
+
+def benchmark_i(jv,qv,m,n,p,d,temperature,beta):
+
+	if temperature == 0.0 and qv == 0.0:
+		return 0.0,1.0
+
+
+	total_states = 2**(m*n*p)
+
+	H = [0.0]*total_states
+	M = [0]*total_states
+	Ms = [0]*total_states
+	prob_M = 0.0
+	prob_Ms = 0.0
+
+	partition_f = 0.0
+
+	r0 = initialize_b(m,n,p)
+
+	sx = r0[0]
+	sy = r0[1]
+
+	if p > 1:
+		sz = r0[2]
+
+	for i in range(total_states):
+		state = bin(i)
+		state = state[2:]
+
+		for j in range(len(state)):
+			row = j/n
+			column = j%m
+			if state[j] == '0':
+				sy[row][column] = 1.0
+			if state[j] == '1':
+				sy[row][column] = -1.0
+		if p == 1:
+			d = '2d'
+			H[i] = system_energy_1(sx,sy,m,n,jv,qv)	
+			r1 = magnet(sx,sy,0.0,m,n,p,d)
+
+		if p > 1:
+			d = '3d'
+			H[i] = system_energy_p(sx,sy,sz,m,n,p,jv,qv)
+			r1 = magnet(sx,sy,sz,m,n,p,d)	
+
+		M[i] = r1[0]
+		Ms[i] = r1[1]
+
+		prob = np.exp(-H[i]*beta)
+
+		partition_f += prob
+
+		prob_M += M[i]*prob
+		prob_Ms += Ms[i]*prob
+
+	M = prob_M/partition_f
+	Ms = prob_Ms/partition_f
+
+	return M, Ms

Monte Carlo - Metropolis/dotproduct.py

@@ -0,0 +1,10 @@
+
+def dotproduct1(x,y,i1,j1,i2,j2):
+	r =  x[i1][j1]*x[i2][j2] + y[i1][j1]*y[i2][j2]
+	return r
+
+def dotproduct_p(x,y,z,i1,j1,k1,i2,j2,k2):
+	r =  x[k1][i1][j1]*x[k2][i2][j2] + y[k1][i1][j1]*y[k2][i2][j2] + z[k1][i1][j1]*z[k2][i2][j2]
+	return r
+
+

Monte Carlo - Metropolis/initializer.py

@@ -0,0 +1,61 @@
+import os
+import time
+import numpy as np
+import random as rand
+
+pi = np.pi
+
+def binary():
+	t = np.random.rand()
+	if t < 0.5:
+		return -1
+	else:
+		return 1
+
+def initializer(tp,d,m,n,p):
+
+	np.random.seed(int(os.urandom(4).encode('hex'), 16))
+
+	if d == '2d':
+		if tp == 'heisenberg':
+			phi = 2.0*pi*np.random.rand(m,n)
+			sx = np.cos(phi)
+			sy = np.sin(phi)
+
+		elif tp == 'ising':
+			sx = 0.0*np.random.rand(m,n) + 1
+			sy = 0.0*np.random.rand(m,n) + 1
+
+			for i in range(m):
+				for j in range(n):
+					sx[i][j] = 0.0
+					sy[i][j] = binary()
+
+		return sx,sy
+
+
+
+	if d == '3d':
+		if tp == 'heisenberg':
+			theta = pi*np.random.rand(p,m,n)
+			phi = 2.0*pi*np.random.rand(p,m,n)
+
+			sx = np.sin(theta)*np.cos(phi)
+			sy = np.sin(theta)*np.sin(phi)
+			sz = np.cos(theta)
+
+		elif tp == 'ising':
+			sx = 0.0*np.random.rand(p,m,n) + 1.0
+			sy = 0.0*np.random.rand(p,m,n) + 1.0
+			sz = 0.0*np.random.rand(p,m,n) + 1.0
+
+			for k in range(p):
+				for i in range(m):
+					for j in range(n):
+						t = binary()
+						sx[k][i][j] = 0.0
+						sy[k][i][j] = t
+						sz[k][i][j] = 0.0
+
+	return sx,sy,sz
+

Monte Carlo - Metropolis/magnet.py

@@ -0,0 +1,38 @@
+import math
+import numpy as np
+
+def magnet(sx,sy,sz,m,n,p,d):
+
+	Mx = np.sum(sx)
+	My = np.sum(sy)
+	if d == '3d':
+		Mz = np.sum(sz)
+
+	Msx = 0.0
+	Msy = 0.0
+	if d == '3d':
+		Msz = 0.0
+
+		for k in range(p):
+			for j in range(n):
+				for i in range(m):
+					Msx += math.pow(-1,(i+j+k))*sx[k][i][j]
+					Msy += math.pow(-1,(i+j+k))*sy[k][i][j]
+					Msz += math.pow(-1,(i+j+k))*sz[k][i][j]
+
+		M = (math.pow(Mx,2) + math.pow(My,2) + math.pow(Mz,2))
+		Ms = (math.pow(Msx,2) + math.pow(Msy,2) + math.pow(Msz,2))			
+
+	elif d == '2d':
+
+		for i in range(m):
+			for j in range(n):
+				Msx += math.pow(-1,(i+j))*sx[i][j]
+				Msy += math.pow(-1,(i+j))*sy[i][j]
+
+		M = math.pow(math.pow(Mx,2) + math.pow(My,2),0.5)/(m*n*p)
+		Ms = math.pow(math.pow(Msx,2) + math.pow(Msy,2),0.5)/(m*n*p)	
+
+
+
+	return M,Ms

Monte Carlo - Metropolis/main.py

@@ -0,0 +1,110 @@
+import sys
+import math
+sys.path.append('/home/abhinav/Desktop/Academics/Summer/Python')
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+
+from spin import *
+from benchmark_ising import *
+
+kb = 1.0
+
+m = 4
+n = 4
+p = 1 
+
+if p > 1:
+	d = '3d'
+elif p == 1:
+	d = '2d'
+
+tp = 'ising'
+
+qv = 0.0
+jv = [0.0]*3
+jv[0] = 1.0
+
+iterations = 250000
+step = 5000
+innerloop = int(2*m*n*p)
+
+print "Model: ",tp
+print "J:", jv
+print "Q:", qv
+print "Total Iterations: ", iterations
+print "Measurement Step: ", step
+print "Inner Loop (Iterations within each Iterations): ", (innerloop/(m*n*p)),"* Total Lattice Points"
+
+temp_i = 0.0
+temp_step = 0.25
+temp_f = 6.0
+
+n0 = int((temp_f - temp_i)/temp_step)
+
+M = [0.0]*n0
+b_M = [0.0]*n0
+Ms = [0.0]*n0
+b_Ms = [0.0]*n0
+M_error = [0.0]*n0
+Ms_error = [0.0]*n0
+rj = [0.0]*n0
+Difference = [0.0]*n0
+
+start_time = time.time()
+
+for i in range(n0):
+#def simulation(i):
+	temperature = temp_i + i*temp_step
+	if temperature == 0.0:
+		beta = math.pow(10.0,50)
+	else:
+		beta = 1.0/(temperature*kb)
+
+	r1 = spin(jv,qv,iterations,m,n,p,d,temperature,tp,step,temp_step,innerloop,beta)
+
+	if tp == 'ising':
+		benchmark = benchmark_i(jv,qv,m,n,p,d,temperature,beta)
+
+	M[i] = r1[0]
+	Ms[i] = r1[1]
+
+	b_M[i] = benchmark[0]
+	b_Ms[i] = benchmark[1]
+
+	M_error[i] = r1[2]
+	Ms_error[i] = r1[3]
+
+	rj[i] = r1[4]
+	Difference[i] = M[i] - Ms[i]
+
+x = np.arange(temp_i,temp_f,temp_step)
+
+
+print "\n-----------------------------------------------\n\nFinal Values for all Temperatures:"
+print "Ferromagnetic Moment, M: ", M
+print "Anti-ferromagnetic Moment, Ms: ", Ms
+
+plt.figure()
+plt.errorbar(x, Ms, yerr=Ms_error, fmt='-o', uplims=True, lolims=True)
+plt.plot(x, b_Ms, 'r-o') 
+plt.title('Ms vs T')
+plt.savefig('Ms vs T.jpg')
+
+plt.figure()
+plt.errorbar(x, M, yerr=M_error, fmt='-o', uplims=True, lolims=True)
+plt.plot(x, b_M, 'r-o')
+plt.title('M vs T')
+plt.savefig('M vs T.jpg')
+
+plt.figure()
+plt.plot(x,rj,'b-o')
+plt.title('Rejectance')
+plt.savefig('Rejectance.jpg')
+
+plt.figure()
+plt.plot(x,Difference,'bo')
+plt.title('Difference b/w M and Ms')
+plt.savefig('Difference.jpg')
+
+print "\nRuntime: %s seconds" % (time.time() - start_time)

Monte Carlo - Metropolis/perturb.py

@@ -0,0 +1,36 @@
+import numpy as np
+import random as rand
+
+def perturb(x0,y0,z0,d,tp):
+	if d == '3d':
+		# Heisenberg Model
+		if tp == "heisenberg":
+			phi = 2.0*np.pi*np.random.rand()
+			theta = np.pi*np.random.rand()
+			x = np.cos(phi)*np.sin(theta)
+			y = np.sin(phi)*np.sin(theta)
+			z = np.cos(theta)
+
+		# Ising Model
+		if tp == "ising":
+			x = -x0
+			y = -y0
+			z = -z0
+
+		return x,y,z
+
+	elif d == '2d':
+		# Heisenberg Model
+		if tp == "heisenberg":
+			phi = 2.0*np.pi*np.random.rand()
+			x = np.cos(phi)
+			y = np.sin(phi)
+
+		# Ising Model
+		if tp == "ising":
+			x = -x0
+			y = -y0
+
+		return x,y
+
+