The difference in magnetisation is therefore only given by the difference \( s_l^1-s_{l}^2=\pm 2 \), or in a more compact way as $$ \begin{equation} M_2 = M_1+2s_l^2, \tag{4} \end{equation} $$ where \( M_1 \) and \( M_2 \) are the magnetizations before and after the spin flip, respectively. Eqs. (3) and (4) are implemented in the function metropolis shown here
void Metropolis(int n_spins, long& idum, int **spin_matrix, double& E, double&M, double *w)
{
// loop over all spins
for(int y =0; y < n_spins; y++) {
for (int x= 0; x < n_spins; x++){
// Find random position
int ix = (int) (ran1(&idum)*(double)n_spins);
int iy = (int) (ran1(&idum)*(double)n_spins);
int deltaE = 2*spin_matrix[iy][ix]*
(spin_matrix[iy][periodic(ix,n_spins,-1)]+
spin_matrix[periodic(iy,n_spins,-1)][ix] +
spin_matrix[iy][periodic(ix,n_spins,1)] +
spin_matrix[periodic(iy,n_spins,1)][ix]);
// Here we perform the Metropolis test
if ( ran1(&idum) <= w[deltaE+8] ) {
spin_matrix[iy][ix] *= -1; // flip one spin and accept new spin config
// update energy and magnetization
M += (double) 2*spin_matrix[iy][ix];
E += (double) deltaE;
}
}
}
} // end of Metropolis sampling over spins