Ising model simulation

Question

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 = %d\n",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();
}

Answer 1

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.

Answer 2

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;

Answer 3

Building on the answers of Quuxplusone and kmdreko, here's a few more issues I noticed. I'll start with a couple of serious ones, before moving on to stylistic issues and potential optimizations:

Dead code

The only thing seriously wrong with your spin() function is that you never call it. It's generally not very useful to ask people to review code you're not using.

(The other, minor thing wrong with it is that it seems to be written with the assumption that the input is a random integer from 1 to 10. Or maybe it's meant to be from 0 to 9, in which case the r > 5 should be r >= 5? But either way is needlessly complicated: since you're only using the input value to make a binary choice, just have it be 0 or 1. Or you could have the input be a float between 0.0 and 1.0, if you'd like to have the ability to bias the initial magnetization state. But of course, none of that makes any difference if you don't actually use the function.)

Use of uninitialized variables

On line 29 of your main function (10 lines below //populate the lattice), you're calling flip() and passing beta as an argument to it. But beta has not been initialized at that point, since you've commented out the line that would do it!

That's a bad thing, and you shouldn't do it. Not only can the value passed to flip() be anything (and possibly vary between runs), leading to unpredictable results, but using the value of an uninitialized variable makes the behavior of your program formally undefined, which means that the compiler is allowed to make it do anything, up to and including making demons fly out of your nose when you run it. A slightly more likely (but still protentially catastrophic) result is that the compiler might simply decide to optimize out your entire main function, since it can prove that there are no code paths through it that don't have undefined behavior!

(Of course, the out-of-bounds array access pointed out by kmdreko also leads to undefined behavior, so just setting an initial value for beta isn't enough to fix your code.)

Confusing variable naming

In flip(), you use H both as an argument to the function and as the name of a local variable. While that's technically allowed, it's very confusing to read. Rename one of them.

Also, your use of temp as the name of a temporary(?) variable in the main function is potentially confusing too, given that the surrounding code deals with (thermodynamic) temperatures. For that matter, temp isn't a particularly informative or appropriate variable name anyway. If you really can't think of a good name for a variable, call it foo or something so that readers can at least tell at a glance that the name means nothing. But in this case, something like iteration would be a decently meaningful name.

Optimization opportunities

First of all, note that your code spends almost all of its time calling flip() repeatedly, so making that function run fast should be your first optimization priority.

Getting rid of the indirect indexing via the tab array, as suggested by kmdreko, is a good start. In the same vein, you might want to turn lattice from a vector of vectors into a proper two-dimensional array, which will eliminated another layer of indirection.

(The down side of using two-dimensional arrays is that you'll have to know the dimensions at compile time, but in your code that's the case anyway. If you want to allow the size of the lattice to be specified at runtime, one option is to dynamically allocate a one-dimensional array of the appropriate size and treat it as a pseudo-multidimensional array by indexing it e.g. as lattice[a*N + b] instead of lattice[a][b].)

You can also trivially save some memory by using signed char instead of int as the type of the lattice elements. You might even want to consider using std::bitset (or implementing your own bit-packed lattice representation) to save even more memory, but this would require you to represent your lattice spin states as 0 and 1 instead of -1 and +1 (which in itself is not necessarily a bad idea at all), and likely comes at the cost of some minor speed reduction and added code complexity. For relatively small values of N, it's probably not worth it.

---

Also, calculating exp(-beta*dE) inside flip() looks like it could potentially be somewhat slow (although, of course, you really ought to profile the code first before spending too much effort on such optimizations). It would be nice to move the exponential calculation out of the inner loop, if we can. Can we?

Expanding the definition of dE earlier in the code (and using the definitions of the neighbor site values val1 to val4 from kmdreko's answer), the argument to exp() is -beta * 2 * s * (H + val1 + val2 + val3 + val4). Here, beta and H do not change within the inner loop, while s and val1 to val4 are always either +1 or -1.

The addition of the external field parameter H to the summed spin of the neighbors makes things a bit more complicated, but we can deal with it by distributing the common -beta * 2 * s term over the sum, giving -beta*2*s*H - beta*2*s*val1 - beta*2*s*val2 - beta*2*s*val3 - beta*2*s*val4. Now, depending on the signs of s and val1...val4, each of these terms can only take one of two values: ±beta*2*H for the first term, and ±beta*2 for the others. If we precalculate the exponentials of each of these values outside the flip() function (and the inner loop that calls it), we can calculate exp(-beta*dE) simply by multiplying these precalculated terms together, e.g. like this:

void flip (float exp_2_beta, float exp_m2_beta, float exp_2_beta_H, float exp_m2_beta_H)
{
    int a = (int)prandom(0, N);
    int b = (int)prandom(0, N);
    int s = lattice[a][b];

    // calculate exp(-beta * 2 * s * (H + sum(neighbors))) based on precalculated
    // values of exp(2*beta), exp(-2*beta), exp(2*beta*H) and exp(-2*beta*H):
    float prob = (s != 1 ? exp_2_beta_H : exp_m2_beta_H);
    prob *= (lattice[a > 0 ? a-1 : N-1][b] != s ? exp_2_beta : exp_m2_beta);
    prob *= (lattice[a < N-1 ? a+1 : 0][b] != s ? exp_2_beta : exp_m2_beta);
    prob *= (lattice[a][b > 0 ? b-1 : N-1] != s ? exp_2_beta : exp_m2_beta);
    prob *= (lattice[a][b < N-1 ? b+1 : 0] != s ? exp_2_beta : exp_m2_beta);

    // flip spin of this site with probability min(prob, 1.0)
    if (prob >= 1 || prob >= prandom(0, 1))
    {
        lattice[a][b] = -s;
    }
}

This makes use of short-circuit evaluation to avoid calling prandom(0, 1) when it's not needed. Rewriting the code like this also avoids making an unnecessary write to lattice[a][b] if the spin doesn't change. (Also, I'm assuming that N and lattice are global constants, since there's really no good reason for them not to be, and that you've fixed prandom() to actually return a float.)

Now, passing all these precalculated factors as parameters to flip() may feel a bit ugly, since they're really just an internal optimization detail. Fortunately, there's a simple fix to that: just move the inner loop (and the precalculation of those values) into the function, e.g. like this:

void do_flips (int count, float beta, float H)
{
    // precalculate some useful exponential terms
    float exp_2_beta = exp(2 * beta), exp_m2_beta = exp(-2 * beta);
    float exp_2_beta_H = exp(2 * beta * H), exp_m2_beta_H = exp(-2 * beta * H);

    // do up to (count) spin flips
    for (int i = 0; i < count; i++) {
        int a = (int)prandom(0, N);
        int b = (int)prandom(0, N);
        int s = lattice[a][b];
    
        // calculate prob = exp(-beta * 2 * s * (H + sum(neighbors)))
        float prob = (s != 1 ? exp_2_beta_H : exp_m2_beta_H);
        prob *= (lattice[a > 0 ? a-1 : N-1][b] != s ? exp_2_beta : exp_m2_beta);
        prob *= (lattice[a < N-1 ? a+1 : 0][b] != s ? exp_2_beta : exp_m2_beta);
        prob *= (lattice[a][b > 0 ? b-1 : N-1] != s ? exp_2_beta : exp_m2_beta);
        prob *= (lattice[a][b < N-1 ? b+1 : 0] != s ? exp_2_beta : exp_m2_beta);
    
        // flip spin of this site with probability min(prob, 1.0)
        if (prob >= 1 || prob >= prandom(0, 1))
        {
            lattice[a][b] = -s;
        }
    }
}

Finally, instead of "sweeping" the lattice periodically to calculate its overall magnetization state, it may be more efficient to just keep track of the sum of the spin states as you update them. For example, you could take the code above and replace the last if statement with:

        // flip spin of this site with probability min(prob, 1.0)
        if (prob >= 1 || prob >= prandom(0, 1))
        {
            lattice[a][b] = -s;
            lattice_sum += -2*s;  // keep track of the sum of all the spins
        }

where lattice_sum is either a global — or, better yet, a local copy of a global variable — or passed into the function as a parameter and returned from it at the end.

Whether such real-time tracking of the total magnetization is faster or slower than periodic sweeping depends on how often you want to sample the magnetization state. With one sample per transition per lattice site, as your current code effectively does, I'd expect the real-time tracking to be somewhat faster, if only because it has better cache locality. Especially so if you make lattice_sum a local variable in do_flips(), which ensures that the compiler will know that it can't be aliased and can be safely stored in a CPU register during the loop.

(Also, the way you're updating the time-averaged magnetization state M looks wonky, and I'm pretty sure you have a bug in it. To test this, try fixing M_sweep to a non-zero constant value and see if M correctly evaluates to M_sweep/(N*N). The way your code is currently written, it doesn't look like it will. Maybe you changed one of your "magic numbers" from 1000 to 500 — or vice versa — and forgot to update the other? One more reason to prefer named constants...)

Answer 4

One thing to do, is to create a lookup table for the float H = exp(-beta*dE) . As your change in energy can only take 4 different values, it is more efficient to store the result of the above exponential for each of those 4 values, and then access it via an index instead of calculating it everytime. (You will actually only need 2 of those values.)

Answer 5

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

As noted by Quuxplusone you should store the mt19937 engine rather than reinstantiate it every time by taking data from std::random_device. You also have the issue that you are initialising a random number generator with a state size of 19937 bits with the unsigned int returned by rd() (likely to be 32 bits on most common architectures). This is clearly an insufficient amount of randomness. You should actually be using a std::seed_seq for this purpose.

std::seed_seq::result_type is std::uint_least32_t. Annoyingly the std::random_device::result_type is an alias of unsigned int, which could potentially be a 16-bit type, so this requires an intermediate distribution, otherwise we run the risk of having zero padding in the supplied integers.

The next question is how many values we need: this can be calculated by multiplying the state size of the generator std::mt19937::state_size by the word size std::mt19937::word_size and dividing by 32. (For std::mt19937 the word size is 32, so you'd be able to just use the state size, but if you decided to replace it with std::mt19937_64 this would become relevant).

Putting that all together gives something along the lines of:

std::mt19937 getengine()
{
    std::random_device rd;

    // uniform int distribution to ensure we take 32-bit values from random_device
    // in the case that std::random_device::result_type is narrower than 32 bits
    std::uniform_int_distribution<std::uint_least32_t> seed_dist{ 0, 0xffffffffu };

    // create an array of enough 32-bit integers to initialise the generator state
    constexpr std::size_t seed_size = std::mt19937::state_size * std::mt19937::word_size / 32;
    std::array<std::uint_least32_t, seed_size> seed_data;
    std::generate_n(seed_data.data(), seed_data.size(), [&rd, &seed_dist]() { return seed_dist(rd); });

    // use the seed data to initialise the Mersenne Twister
    std::seed_seq seq(seed_data.begin(), seed_data.end());
    return std::mt19937{ seq };
}

This should take enough data from std::random_device to seed the generator. Unfortunately there are implementations where std::random_device is deterministic (MinGW is the main culprit here) and there is no reliable way to detect this. The entropy() method on std::random_device should return zero for a deterministic implementation, unfortunately there are non-deterministic implementations that also return zero (e.g. LLVM libc++), which is rather unhelpful.

There are also various subtle flaws in std::seed_seq to be aware of, see Melissa O'Neill's article on C++ Seeding Surprises for details: basically it turns out that even using an (apparently) sufficient amount of data, there are still various states that cannot be output.