Ising model Monte Carlo simulations in 4D and 5D

Question

I'm going to be simulating the Ising Model in 4D and above to calculate spin-spin correlations and critical exponents and am wondering how to tackle this algorithmically.

For example, in 1D, use an array to store spins. In 2D, use a matrix. In 3D, use a 2D matrix, or a three-tuple. One way I think about it is that in any dimension N, you can backtrace to the 3rd dimension to help understand nearest neighbor terms.

For example, in 2D, the spins next a particular spin in the center are just 2 1D nearest neighbors. In 3D, the spins next to a particular spin in the center are just 2 2D nearest neighbors (ie the planes perpendicular to each other).

In 4D, can I think of it as the 2 3D nearest neighbors? Similarly, how can I think of the 5D algorithmically?

Or, should I refrain from using the standard container classes and just use tuples everywhere?

Answer 1

The best way to store them may be to use a 1D array and use a separate list of bonds which keeps track of the connections between the sites. In the approach, all you need to do to change your code from 1D to 2D to 3D is to modify the bonds array.

As an example, in 2D, you could have a list of the spin values called spins[N] in 2D there are 2N 'bonds' (each site is connected to four others, but we have to divide by 2 to avoid overcounting). In this segment of C++ code, we build the bonds array (assuming we already have an empty array: bonds[2*N][2]

for (int i = 0; i < nsites; i++) {
        //calculate the x and y coordinates of site i
        int xi = i % lx; 
        int yi = i/lx;
        //x direction bonds are even numbered
        bonds[0][2*i] = i;
        bonds[1][2*i] = yi*lx + ((xi+1) %lx);
        //y direction bonds are odd numbered
        bonds[0][2*i+1] = i; // same
        bonds[1][2*i+1] = ((yi+1)%ly)*lx + xi; 
    }