Bdtyktl

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?

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

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;

}

Function 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;

    }*/



    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;





        std::cout<<beta<<"n";

        for(i=0;i<=sweep;i++)

        {

        //T=0.1*i;





            //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";

    }





    for(i=0;i<N;i++)

    {

        for(j=0;j<N;j++)

        {

            std::cout<<lattice[j][i]<<"t";

        }

        std::cout<<std::endl;

    }

    outdata.close();

}

There are many, many things wrong with your code; I'll list some, and then I advise you to fix as many of them as you can (and more) and then repost a new question with the revised code.

First of all, you should compile your code with -W -Wall -Wextra (or -W4 if you use MSVC). Read the first warning diagnostic. Fix it. Recompile. Read the first diagnostic. Fix it. ...until all the warnings are completely gone.

Procedural nit: When you post code to CodeReview, please post it in one big cut-and-pasteable chunk, or at least just a couple of chunks. Your current post mixes code and commentary in a way that makes it hard for the reviewer to cut-and-paste the whole piece of code into a file for testing.

Speaking of commentary:

    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;

Is beta supposed to be 1.0 / T, or not? By commenting out that line of code, you make it invisible to the compiler. Better make it invisible to the reader, too, and just delete it! (If you're commenting things out to preserve the "history" of the code, read up on version control systems like git, and then use one. It's easy!)

Furthermore, since you don't initialize beta, you can probably just get rid of the variable entirely.

Finally, the idiomatic way to place the whitespace when you're defining a variable of type std::vector<std::vector<int>> is simply thus:

std::vector<std::vector<int>> lattice;

Notice the space between the variable's type and its name; and the lack of space anywhere else.

Populating that lattice can be done quickly and easily using vector's "filling" constructor:

std::vector<std::vector<int>> lattice(N, std::vector<int>(N, -1));

Now you don't need those nested for loops anymore!

Going back up to the top of your code:

int spin(int r) 

{

    int s;

    if(r>5)

    {

        s=+1;

    }

    else

    {

        s=-1;

    }



    return s;

}

Replace this 14-line function with a 4-line function that does the same thing more clearly and simply:

int spin(int r)

{

    return (r > 5) ? 1 : -1;

}

No local variables, no mutation, no startling use of the =+ "operator"; and perhaps most importantly for the reader, there's now 10 more lines available on my screen so I can look at other code at the same time! Vertical real estate can be important for reading comprehension.

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;

}

This is wrong in at least two ways. First, most importantly: you're constructing a std::random_device on each call to this function. That is extremely expensive. Think of std::random_device as an open file handle to /dev/urandom, because that's what it is, under the hood. That means every time you call prandom, you're opening a file, reading some bytes, and closing it again!

You should keep a global random number generator, initialized just once at the start of the program. One way to do this (not cheap, but not as expensive as opening a file on every call to prandom) would be

float prandom(int i,int N)

{

    static std::mt19937 gen = () {

        std::random_device rd;

        return std::mt19937(rd());

    }();

    std::uniform_real_distribution<float> dis(i, N);

    // ...

Notice that uniform_real_distribution<> is secretly an alias for uniform_real_distribution<double>. Your code doesn't use doubles; it uses floats. So it's (always) better to be explicit and say what type you mean — you have less chance of getting it wrong by accident!

    int t = dis(gen);

    return t;

...And then you go ahead and stuff the return value into an int anyway! So what was the point of using a uniform_real_distribution in the first place? And what's the point of returning a float from prandom? Did you mean to simply

return dis(gen);

instead?

You're also seeding your PRNG wrong, but seeding it correctly is a huge headache in C++17, so never mind that.

        if(temp>15)

        {

        H=H-0.015;

        }

Please indent your code correctly. You can use the clang-format command-line tool to automatically indent everything, or if you use a graphical IDE it almost certainly has an "Indent" option somewhere in the menus.

That's enough for one day. As I said above, I advise you to fix as much as possible (that is, fix everything I talked about here, and then fix everything else you can think of, too) and then repost.

After you fix everything, but before you repost, read your code from top to bottom one more time! Find two more things that need fixing, and fix them. Then repost.

I'll build off of Quuxplusone's answer.

You have issues with tab.

int tab[N];

tab[0] = N-1;

tab[N+1] = 0; // out of bounds

Only the elements 0 through N-1 are available, so this assignment is wrong. Your initialization loop also attempts to initialize tab[N]. And in flip, since a and b are randomly generated between 0 and N-1, getting tab[a+2] can be N+1 at maximum, also out of bounds.

While many compilers will accept it, you're not allowed to create an array with a non-const int. However, all these issues can be fixed pretty easily; just create it with a bit of extra room.

const int N = 20; // set to const

int tab[N+2];

As a side note, I would personally not use tab at all! I would simply do the wraparound logic within flip():

int val1 = lattice[a > 0 ? a-1 : N-1][b];

int val2 = lattice[a < N-1 ? a+1 : 0][b];

int val3 = lattice[a][b > 0 ? b-1 : N-1];

int val4 = lattice[a][b < N-1 ? b+1 : 0];

That way, you don't have to worry about creating it or passing it as a parameter at all.

Your flip function needlessly copies lattice each call. You should pass it as a reference.

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

                                            // ^

{

    ...

    return;

}

I've also made the function return void since you don't need to reassign lattice since its passed by reference and just call it like so:

flip(N, lattice, beta, tab, H);

Make a separate function like print_lattice. It makes it clear whats happening without even looking at it and you can use it in multiple places (like the one you've commented out).

Don't declare your loop variables outside the loop (unless you need to). And get rid of any unused variables.

inta,b,N=20,i,j,k,r,t,sweep=1500;

Use better variable names. Consider using x and y instead of a and b since they'd be more immediately understandable. And give a more descriptive name to H and M, I'm not familiar with the algorithm and these names don't help much.

for(int u=0;u<N;u++)

{

    if(i>=500)

    {M_sweep=M_sweep+lattice[t][u];}

}

This loop checks for i but never modifies it, this check can be moved outside the loop like so:

if (i >= 500)

{

    for(int u=0;u<N;u++)

    {

        M_sweep = M_sweep + lattice[t][u];

    }

}

On second look, it can moved even higher to outside the t loop.

You use lots of constants:

H=H+0.015;

...

T=2.2;

Consider giving these constants names.

const float H_STEP = 0.015f;