import scipy.spatial.distance as spd import numpy as np import matplotlib.pyplot as plt %matplotlib inline # used to compare energy before and after a move # we assumed that epsilon and alpha are 1 def local_energy(site): local_sum = 0 for i in xrange(len(system)): if site != i: tempx = system[site][0] - system[i][0] tempy = system[site][1] - system[i][1] if tempx > 100./2: tempx = 100. - tempx if tempy > 100./2: tempy = 100. - tempy #r = spd.euclidean(system[i][0],system[j][0]) dist = np.sqrt(tempx**2+tempy**2) #dist = spd.euclidean(system[site],system[i]) local_sum += -((1/dist)**12-(1/dist)**6) return local_sum # total energy of the system def system_energy(): total_sum = 0 for i in xrange(len(system)): total_sum += local_energy(i) return total_sum def radial_number_density(): rad_dist = [] for i in xrange(len(system)): for j in xrange(len(system)): if i != j: temp = spd.euclidean(system[i],system[j]) rad_dist.append(temp) return rad_dist n = 100 # number of particles # initializing the system # the atoms occupy a square of size 10X10 system = np.zeros([n,2]) for i in xrange(len(system)): system[i][0] = i%10 # np.sqrt(n) wont work as it will return a float system[i][1] = i/10 # we need an integer denominator for this to work! global_energy = [] accepted_moves = 0 plt.subplot(131) alist = radial_number_density() anarray = np.asarray(alist) plt.hist(anarray,10); initial_dist = radial_number_density() T = 1. # T = np.array([5.,5./2,1./1./2,1./10]) beta = 1./T # B = 1./T #for beta in B: for time in xrange(5000): site = np.random.choice(np.arange(n)) # atom to be moved # this atom can now move in the range (-alpha,alpha) in x & y #xmove = np.random.choice(np.array([-2.,-5./3,-4./3,-1,-2./3,-1./3,0,1./3,2./3,1.,4./3,5./3,2.])) #ymove = np.random.choice(np.array([-2.,-5./3,-4./3,-1,-2./3,-1./3,0,1./3,2./3,1.,4./3,5./3,2.])) #xmove = np.random.choice(np.array([0,1./3,2./3,1.,4./3,5./3,2.])) #ymove = np.random.choice(np.array([0,1./3,2./3,1.,4./3,5./3,2.])) xmove = np.random.choice(np.array([0.,1.,2.,3.,4.])) ymove = np.random.choice(np.array([0.,1.,2.,3.,4.])) pre_move_e = local_energy(site) system[site] = system[site] + np.array([xmove,ymove]) if system[site][0] > 100 : system[site][0] = system[site][0]%100 #if system[site][0] < 0 : # system[site][0] = system[site][0]%100 if system[site][1] > 100 : system[site][1] = system[site][1]%100 #if system[site][1] < 0 : # system[site][1] = system[site][1]%100 post_move_e = local_energy(site) if post_move_e - pre_move_e < 0: global_energy.append(system_energy()) accepted_moves += 1 else : temp = np.random.random() if temp < np.exp(-beta*(post_move_e - pre_move_e)): global_energy.append(system_energy()) accepted_moves += 1 else : system[site] = system[site] - np.array([xmove,ymove]) global_energy.append(system_energy()) # we get global_energy, accepted moves as the outputs here. plt.subplot(132) alist = radial_number_density() anarray = np.asarray(alist) plt.hist(anarray,10); plt.subplot(133) plt.plot(global_energy) print "accepted_moves", accepted_moves final_dist = radial_number_density() #print plt.hist(final_dist)[0], plt.hist(final_dist)[1] temp_array_range_initial = plt.hist(initial_dist)[1] temp_array_number_initial = plt.hist(initial_dist)[0] for i in xrange(len(temp_array_range_initial)-1): denominator = temp_array_range_initial[i+1]**2 - temp_array_range_initial[i]**2 temp_array_number_initial[i] = temp_array_number_initial[i]/denominator ####### temp_array_range_final = plt.hist(final_dist)[1] temp_array_number_final = plt.hist(final_dist)[0] for i in xrange(len(temp_array_range_final)-1): denominator = temp_array_range_final[i+1]**2 - temp_array_range_final[i]**2 temp_array_number_final[i] = temp_array_number_final[i]/denominator plt.cla() #plt.plot(temp_array_number_initial,'o',temp_array_number_final,'.') plt.plot(temp_array_number_final) x = [] y = [] for i in xrange(len(system)): x.append(system[i][0]) y.append(system[i][1]) plt.scatter(x,y); import scipy.spatial.distance as spd import numpy as np import matplotlib.pyplot as plt %matplotlib inline # used to compare energy before and after a move # we assumed that epsilon and alpha are 1 def local_energy(site): local_sum = 0 for i in xrange(len(system)): if site != i: dist = spd.euclidean(system[site],system[i]) local_sum += -((1/dist)**12-(1/dist)**6) return local_sum # total energy of the system def system_energy(): total_sum = 0 for i in xrange(len(system)): total_sum += local_energy(i) return total_sum def radial_number_density(): rad_dist = [] for i in xrange(len(system)): for j in xrange(len(system)): if i != j: temp = spd.euclidean(system[i],system[j]) rad_dist.append(temp) return rad_dist n = 100 # number of particles # initializing the system # the atoms occupy a square of size 10X10 system = np.zeros([n,2]) for i in xrange(len(system)): system[i][0] = i%10 # np.sqrt(n) wont work as it will return a float system[i][1] = i/10 # we need an integer denominator for this to work! global_energy = [] global_energy.append(system_energy()) energy_now = system_energy() #### accepted_moves = 0 plt.subplot(131) alist = radial_number_density() anarray = np.asarray(alist) plt.hist(anarray,10); T = 1. # T = np.array([5.,5./2,1./1./2,1./10]) beta = 1./T # B = 1./T #for beta in B: for time in xrange(5000): site = np.random.choice(np.arange(n)) # atom to be moved # this atom can now move in the range (-alpha,alpha) in x & y #xmove = np.random.choice(np.array([-2.,-5./3,-4./3,-1,-2./3,-1./3,0,1./3,2./3,1.,4./3,5./3,2.])) #ymove = np.random.choice(np.array([-2.,-5./3,-4./3,-1,-2./3,-1./3,0,1./3,2./3,1.,4./3,5./3,2.])) #xmove = np.random.choice(np.array([0,1./3,2./3,1.,4./3,5./3,2.])) #ymove = np.random.choice(np.array([0,1./3,2./3,1.,4./3,5./3,2.])) xmove = np.random.choice(np.array([0.,1.,2.,3.,4.])) ymove = np.random.choice(np.array([0.,1.,2.,3.,4.])) pre_move_e = local_energy(site) system[site] = system[site] + np.array([xmove,ymove]) if system[site][0] > 100 : system[site][0] = system[site][0]%100 #if system[site][0] < 0 : # system[site][0] = system[site][0]%100 if system[site][1] > 100 : system[site][1] = system[site][1]%100 #if system[site][1] < 0 : # system[site][1] = system[site][1]%100 post_move_e = local_energy(site) diff = post_move_e - pre_move_e ### if post_move_e - pre_move_e < 0: energy_now += diff #### global_energy.append(energy_now) #### accepted_moves += 1 else : temp = np.random.random() if temp < np.exp(-beta*(post_move_e - pre_move_e)): energy_now += diff #### global_energy.append(energy_now) #### accepted_moves += 1 else : system[site] = system[site] - np.array([xmove,ymove]) global_energy.append(energy_now) #### # we get global_energy, accepted moves as the outputs here. plt.subplot(132) alist = radial_number_density() anarray = np.asarray(alist) plt.hist(anarray,10); plt.subplot(133) plt.plot(global_energy) print "accepted_moves", accepted_moves x = [] y = [] for i in xrange(len(system)): x.append(system[i][0]) y.append(system[i][1]) plt.scatter(x,y); def local_energy(site): local_sum = 0 for i in xrange(len(system)): if site != i: tempx = system[site][0] - system[i][0] tempy = system[site][1] - system[i][1] if tempx > 50: if temp y > 50 : dist = spd.euclidean(system[site],system[i]) local_sum += -((1/dist)**12-(1/dist)**6) return local_sum