Ttykuu

I have written this code to simulate Ising Model at one particular temperature in presence of magnetic field to observe hysteresis effect using the metropolis algorithm.

While the code runs and gave me a desired output, it is a badly written code(I feel so) because of my lack of coding experience. Could you help me out as to how could I have written this in a better way? Or what can I do to optimize it next time I write such a code? Or are there any inbuilt functions I could have called instead of writing a block?

P.S. I borrowed the random number generator directly from someone's answer to a thread on this site. I cannot find the exact thread, apologies! (I will cite it once I do)

I hope this kind of a question is allowed on this thread. Any help or resources would be appreciated.

//function to assign random spins to the lattice

int spin(int r) 

 int s;
 if(r>5)
 
 s=+1;
 
 else
 
 s=-1;
 

 return s;

////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////
//random number generator

float prandom(int i,int N)

 std::random_device rd; //Will be used to obtain a seed for the random number engine
 std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd()
 std::uniform_real_distribution<> dis(i,N);
 // Use dis to transform the random unsigned int generated by gen into a 
 // double in [1, 2). Each call to dis(gen) generates a new random double
 int t = dis(gen);
 return t;

////////////////////////////////////////////////////////////////////////////
//frunction to randomly flip the spins

I am selecting a random site and seeing how the total energy will change by calculating the energy dE for the neighbouring sites. If the energy is negative then I make the spin flip permanent if it is not negative then I gave a probability exp(-dE) by which it can flip the spin

std::vector< std::vector < int > > flip (int N,std::vector< std::vector < int > >lattice, float beta,int tab[],float H)

 int a =prandom(0,N);
 int b =prandom(0,N);
 int s=lattice[a][b];
 float dE=(2*s*H)+(2*s*(lattice[tab[a+2]][b]+lattice[tab[a]][b]+lattice[a][tab[b+2]]+lattice[a][tab[b]]));
 //std::cout<<dE<<"t"<<a<<"t"<<b<<"n";
 if(dE<0)
 
 s=-1*s;
 
 else
 
 float k = 1.0*prandom(0.0,1000)/1000;
 float H = exp(-beta*dE);
 if(k<=H)
 
 s=-1*s;
 
 else
 
 s = 1*s;
 
 
 lattice[a][b]=s;
 return lattice;

////////////////////////////////////////////////////////////////////////////

// main program///////
 int main()
 
 std::ofstream outdata;
 outdata.open("ising_model_field_final2.txt");
 int a,b,N=20,i,j,k,r,t,sweep=1500;
 float M=0,M_sweep=0,H=-0.10;
 int tab[N];
 tab[0] = N-1;
 tab[N+1] = 0;
 for (i=1;i<=N;i++)
 
 tab[i]=i-1; // this is the periodic boundary condition to make my lattice infinite (lattice site [x][0] is a neighbour of [x][N] and so on..)
 
 float T, beta;
 //beta=1.0/T; // boltzman constant is assumed to be 1.
 //creating a 2d lattice and populating it
 std::vector< std::vector < int > >lattice;

 //populate the lattice
 for (i=0; i<N; i++)
 
 std::vector< int > row; //create a row of the lattice
 for (j=0;j<N;j++)
 
 row.push_back(-1); //populate the row vector
 
 lattice.push_back(row); //populate the column vector
 
 lattice=flip(N,lattice,beta, tab,H);
 /* for(i=0;i<N;i++)
 
 for(j=0;j<N;j++)
 
 std::cout<<lattice[j][i]<<"t";
 
 std::cout<<std::endl;
 */
 ///////////////////////////////////////////////////////////////////////////
 //field control
 for(int temp=1;temp<=30;temp++)
 
 if(temp>15)
 
 H=H-0.015;
 
 else
 
 H=H+0.015;
 
 //M=0;
 T=2.2;
 beta=1.0/T;
 ///////////////////////////////////////////////////////////////////////////
 //number of sweeps
 std::cout<<beta<<"n";
 for(i=0;i<=sweep;i++)
 
 //T=0.1*i;
 ///////////////////////////////////////////////////////////////////////////
 //Number of flips
 //printf("Sweep = %dn",i);
 for(j=1;j<=N*N;j++)
 
 lattice=flip(N,lattice,beta, tab,H);
 
 M_sweep=0;
 for(t=0;t<N;t++)
 
 for(int u=0;u<N;u++)
 
 if(i>=500)
 M_sweep=M_sweep+lattice[t][u];
 
 //std::cout<<"Mag="<<M<<"t";
 
 M=M+ M_sweep/(N*N);
 //std::cout<<"Mag="<<M<<"t";
 
 M=M/(sweep-1000);
 std::cout<<T<<"n";
 outdata << M <<"t"<< H <<"n";
 
 //////////////////////////////////////////////////////////////////////////// 
 //printing the output of 2d lattice
 for(i=0;i<N;i++)
 
 for(j=0;j<N;j++)
 
 std::cout<<lattice[j][i]<<"t";
 
 std::cout<<std::endl;
 
 outdata.close();