# How to get free energy surface of a 3d Ising Spin system using Monte Carlo simulation?

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.

As far as I can see, you are not applying any bias to your simulation. The only place where the bias Q comes into play is in the calcEnergy method, which is never called by your code.