Glauber dynamics for the 2D Ising model

NN  =  64;
TT  =  100;

h  =  0;
J  =  1;
kT  =  
      1;      (*2.269185;*) 


energydiff[x_, 
        y_]  := 
N[-2 h  -  
          J (config[Mod[x - 1, NN], 
                  Mod[y, NN]]  + 
config[Mod[x + 1, NN], 
                  Mod[y, NN]]  + 
config[Mod[x, NN], 
                  Mod[y - 1, NN]]  +  
config[Mod[x, NN], 
                  Mod[y + 1, 
                    NN]])
];

Do[
config[x, y]  =  
          2  Random[Integer] - 1, 
{x, 0, NN - 1}, {y, 0, 
          NN - 1}
];
Show[
Graphics[
RasterArray[
Table[
Hue[(config[i, j] + 1)/
                  4], 
{i, 0, NN - 1}, {j, 0, 
                NN - 1}
]
]
], 
AspectRatio ->  
          1
];

Do[
Do[
x  =  
            Random[Integer, {0, NN - 1}]; 
y  =  
            Random[Integer, {0, NN - 1}]; 
deltaE  =  
            energydiff[x, y]; 
prob  =  
            1/(1  +  E^(deltaE/kT)); 
If[
            Random[] ≤ prob, config[x, y] = 1, 
            config[x, y] = -1], 
{s, 1, 
            NN^2}
]; 
Show[
Graphics[
RasterArray[
Table[
Hue[(config[i, j] + 1)/
                    4], 
{i, 0, NN - 1}, {j, 0, 
                  NN - 1}
]
]
], 
AspectRatio ->  
            1
], 
{t, 1, 
          TT}
];