# Ising model Monte Carlo simulations in 4D and 5D

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?

• In 3d there are 3 perpendicular planes of 2d nearest neighbors. In 4d, there would be 4 cubes of 3d nearest neighbors. And so on.
– d_b
Mar 9, 2015 at 5:26
• To prove user37496's statement, consider choosing three axes from your set of four axes. Three defines a volume, so the number of (N-1)-hyperplanes mutually perpendicular passing through a point is $N$ choose $N-1$, that is, $N$.
– user12029
Mar 9, 2015 at 6:25

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