I am doing a Monte Carlo Simulation of the properties of a 3D Ising Spin system. I want to get the free energy surface of the spin system from the simulation. It is a magnetization vs free energy curve. Free energy is defined as F(m)=-kT ln(P(m)) where m is magnetization and P is the probability of occupying the corresponding magnetization. Suppose, no external magnetic field is applied. Due to the very high energy barrier at around m=0, I have employed Umbrella sampling. But the result doesn't seem to be as expected. I don't get the stable equilibrium positions(minima at -1 and 1) with bias strength 500. For T
The expected result is like below,
What I get is that,
To find the P(m) we use the frequency histogram data from mcmove operation. (see the code below)
# -*- coding: utf-8 -*- """ Created on Sun Mar 22 19:28:53 2020 @author: Endeavour """ from numpy.random import rand import numpy as np import matplotlib.pyplot as plt import seaborn as sns Q = 10 #The bias strength def initialise(N): ''' generates a random spin spin_statesuration for initial condition''' spin_states = np.random.choice([1, -1], size=(N, N,N)) #spin_states=np.ones([N,N,N]) return spin_states def calcMag(spin_states): '''Magnetization of a given spin_statesuration''' mag = (np.sum(spin_states)) return (mag/(N*N*N)) def calcEnergy(spin_states): '''Energy of a given spin_statesuration''' energy = 0 for i in range(len(spin_states)): for j in range(len(spin_states)): for k in range(len(spin_states)) : s = spin_states[i,j,k] energy += -s*find_neighbours(spin_states,N,i,j,k) energy+=Q*(calcMag(spin_states))**2 return energy/6. def find_neighbours(spin_states,N,x,y,z): left =spin_states[x,(y-1)%N,z] right =spin_states[x,(y+1)%N,z] top =spin_states[(x-1)%N,y,z] bottom =spin_states[(x+1)%N,y,z] front =spin_states[x,y,(z+1)%N] back =spin_states[x,y,(z-1)%N] tot_spin=left+right+top+bottom+front+back return (tot_spin) def mcmove(spin_states, beta): '''Monte Carlo move using Metropolis algorithm ''' cost=calcEnergy(spin_states) #Store initial energy for x in range(len(spin_states)): for y in range(len(spin_states)): for z in range(len(spin_states)): x = np.random.randint(len(spin_states)) y = np.random.randint(len(spin_states)) z = np.random.randint(len(spin_states)) s = spin_states[x,y,z] s*=-1 cost = calcEnergy(spin_states)-cost print(cost) if rand() < np.exp(-cost*beta): #if cost<0 exp(-cost*beta) should be >1 s *= -1 spin_states[x, y,z] = s return spin_states #-------------------------Simulation Parameters---------- nt = 70 # nt>20 20 points will be between 4 &5 number of temperature points N = 16 # size of the lattice, size x size eqSteps = 10 # number of MC sweeps for equilibration mcSteps = 100 # number of MC sweeps for calculation Temp =4.6 #Temperature #------------------------------------------------------------- dat_M= #Magnetisation data for mc trials mag= #Abscissa for plotting calculated from bin values spin_states = initialise(N) iT=1.0/Temp#Inverse Temperature for beta for i in range(eqSteps): # equilibrate mcmove(spin_states, iT) # Monte Carlo moves if(i%5==0): print ('ok trial') for i in range(mcSteps): mcmove(spin_states, iT) Mag = calcMag(spin_states)/(N*N*N) # calculate the magnetisation dat_M.append(Mag) print(dat_M) divs=np.arange(-1,1,0.005) #Find out bin values sns.distplot(dat_M,kde=0,bins=divs) #Plot the histogram hist_counts,bin_edges=np.histogram(dat_M,bins=divs) #Take the count(freq) for i in range(len(divs)-1): mean=(divs[i]+divs[i+1])/2 #Aerage magnetisation for each bin mag.append(mean) W=Q*np.square(mag) #Bias potential free_energy= for y in (hist_counts): if (y!=0) : free_energy.append(-(1/iT)*np.log(y/mcSteps) ) #-kTlog(P(m)) else: free_energy.append(np.nan) free_energy=np.subtract(free_energy,W) #Correction for bias potetial #Plotting plt.figure(num=2,figsize=(10,6),dpi=80, facecolor='w', edgecolor='b') plt.xlabel('x') plt.ylabel('$-kTlogP_i$') plt.title('Free Energy surface Q=%d,Temperature=%f'%(Q,Temp)) #plt.plot(mag,free_energy,'-',label='Potential Strength(Q): %d \n Biasing strength($x^m$):%d'%(Q,m)) plt.plot(mag,free_energy,'+',label='Test') plt.legend()
I can't find out any problem with the method. Maybe there is some problem with code. Please can you point out any error? Any hint or solution will be of great help. I can't get two minimas at all. I have tried with different bias strength but no success.
Thank you in advance.
After making some corrections, the current code seems to be in an infinite loop. I have printed cost but it gives only 2 values -2.66664719581604 & 0.0.