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,

enter image description here

What I get is that,

enter image description here

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))
    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)
    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]


    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]
                cost = calcEnergy(spin_states)-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

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

W=Q*np.square(mag)                      #Bias potential 

for y in (hist_counts):
    if (y!=0) :
        free_energy.append(-(1/iT)*np.log(y/mcSteps)  )    #-kTlog(P(m))

free_energy=np.subtract(free_energy,W)  #Correction for bias potetial

plt.figure(num=2,figsize=(10,6),dpi=80, facecolor='w', edgecolor='b')
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))

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.


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.

In general, the bias need to be used in the actual MC code when you compute the acceptance probability of your move.

  • $\begingroup$ Thank you for the point. I tested the edited code but it till not gives the correct results. Besides,I came to know that this problem can be solved without using a bias function. But my code fails to do so. Maybe there is some problem with my MC implementation. I don't have any clue. It should work without any bias. $\endgroup$ Mar 23 '20 at 14:57
  • $\begingroup$ Can you update your question then? Since you don't have any bias you don't have to subtract W from the free energy profile. Can you show us the new updated plot? $\endgroup$
    – lr1985
    Mar 23 '20 at 15:17
  • $\begingroup$ The code is now updated. I can't give you a plot as it currently produces none. It seems to be stuck in an infinite loop. I have printed cost but it gives only 2 values -2.66664719581604 & 0.0 alternatively. Bias strength is 10. $\endgroup$ Mar 24 '20 at 4:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.