Message Boards Message Boards

[WSS17] Construction of the ground state in the generalized Ising model

GROUPS:

GOAL OF THE PROJECT:

The goal of the project is to understand what kind of ground states are realized for a generalized Ising model with a given set of energy weights.

SUMMARY OF WORK:

  1. Starting from random distributions states with local minimal energy were obtained via Monte Carlo (MC) and Cellular automaton (CA) techniques.
  2. Exhaustive search for small field sizes was done.

The 2 D Ising model is a system of interacting spins placed in the integer points of a rectangular grid. The energy of the system s is given by the list convolution of the weight pattern and the list of spin values s

e[s_, patternNum_] := Apply[Plus, s ListConvolve[pattern, s, 2], {0, 1}]

The weight pattern is defined as

pattern = 1/2{{ w4,w1,w3},{w2,0,w2},{w3,w1,w4}} 

and the values for w's span {-2,-1,0,1,2}.

Monte Carlo

One can implement MC choosing a spin randomly and flipping it if it reduces the energy (Metropolis algorithm). To speed up one may choose a random spin within a square which is m times less than the whole field and m-periodically construct its copies to span the whole field. Then one calculates energy shifts from the (m/n)^2 spin-flips and flips only those which decrease the energy.

ClearAll[n, field, m, w1, w2, w3, w4, \[CapitalDelta]e1, position, positions, op];

n = 2^6; field = RandomInteger[1, {n, n}] /. 0 -> -1; w1 = 2; w2 = 2; w3 = 0; w4 = 0;

\[CapitalDelta]e1 = position \[Function] ((field[[Mod[#1 - 1, n, 1], Mod[#2 + 1, n, 1]]] 
+ field[[Mod[#1 + 1, n, 1], Mod[#2 - 1, n, 1]]]) w4 
+ (field[[Mod[#1 - 1, n, 1], Mod[#2 - 1, n, 1]]] 
+ field[[Mod[#1 + 1, n, 1], Mod[#2 + 1, n, 1]]]) w3 
+ (field[[Mod[#1 - 1, n, 1], #2]] + field[[Mod[#1 + 1, n, 1], #2]])w1 
+ (field[[#1, Mod[#2 - 1, n, 1]]] + field[[#1, Mod[#2 + 1, n, 1]]])w2 ) field[[#1, #2]] &[Sequence @@ position];

 positions[m_] := With[{position = RandomInteger[{1, n/m}, 2]}, 
  If[Or @@ Flatten@Table[If[\[CapitalDelta]e1[{#1, #2}] > 0, field[[#1, #2]] *= -1; True, False] &
[Mod[position[[1]] + n/m i, n, 1], Mod[position[[2]] + n/m j, n, 1]], {i, 0, m - 1}, {j, 0, m - 1}], Sow[field]]]

Monitor[op = Reap[Do[positions[16], {monitor, 1, 500}]][[2]];//AbsoluteTiming, monitor]
Tooltip[ArrayPlot[field /. -1 -> 0, PixelConstrained -> 2], "{2,2,0,0}"]
Manipulate[ArrayPlot[op[[1, l]] /. -1 -> 0, PixelConstrained -> 2], {l,1,Length[op[[1]]], 1}, SaveDefinitions -> True]

MC500steps

Cellular Automata

One may build a cellular automaton locally reducing (not increasing) energy at each step. To this end one builds a rule based on the pattern

ClearAll[patterns]
stepNumber = 20; n = 2^6; field = RandomInteger[1, {n, n}] /. 0 -> -1;
patterns = {{{0, 1, 0},{1, 0, 1},{0, 1, 0}},  -{{0, 1, 0},{1, 0, 1},{0, 1, 0}},  {{0, 1, 0},{2, 0, 2},{0, 1, 0}},
 -{{0, 1, 0},{2, 0, 2},{0, 1, 0}},  {{0, 1, 0},{-1, 0, -1},{0, 1, 0}}, {{0, 1, 0},{-2, 0, -2},{0, 1, 0}},
 -{{0, 1, 0},{-2, 0, -2},{0, 1, 0}},  {{0, 1, 0},{0, 0, 0},{0, 1, 0}},  -{{0, 1, 0},{0, 0, 0},{0, 1, 0}}};

DecreaseEnergy[matrix_, pattern_] :=  If[matrix[[2, 2]] Plus @@ Flatten[matrix pattern] > 0, -matrix[[2, 2]], matrix[[2, 2]]]
DecreaseEnergy1[matrix_, pattern_] :=  If[matrix[[2, 2]] Plus @@ Flatten[matrix pattern] >= 0, -matrix[[2, 2]], matrix[[2, 2]]]

tuin = Tuples[{-1, 1}, {3, 3}];

In fact one can build two rules based on the decision to flip or not to flip if the energy does not change.

tuout = Table[Thread[tuin -> (DecreaseEnergy[#, patterns[[i]]] & /@ tuin)] /. -1 -> 0, {i, 1, Length[patterns]}];
tuout1 = Table[Thread[tuin -> (DecreaseEnergy1[#, patterns[[i]]] & /@ tuin)] /. -1 -> 0, {i, 1, Length[patterns]}];

rule[matrix_, a_, ruleNumber_] := matrix /. tuout[[ruleNumber]]
rule1[matrix_, a_, ruleNumber_] := matrix /. tuout1[[ruleNumber]]

caFunction0[ruleN_] := CellularAutomaton[{rule[#1, #2, ruleN] &, {}, {1, 1}}, RandomInteger[1, {n, n}], stepNumber];
caFunction1[ruleN_] := CellularAutomaton[{rule1[#1, #2, ruleN] &, {}, {1, 1}}, RandomInteger[1, {n, n}], stepNumber];

evolution = caFunction0 /@ Range[Length[patterns]];
evolution1 = caFunction1 /@ Range[Length[patterns]];

e[s_, patternNum_] := Apply[Plus, s ListConvolve[patterns[[patternNum]], s, 2], {0, 1}]

allenergies = Table[e[evolution[[j, i]], j], {j, 1, Length[patterns]}, {i, 1, stepNumber}];

Manipulate[Grid[Partition[Row[{
ArrayPlot[evolution[[#, i]], PixelConstrained -> 2, PlotLabel -> "E = " <> ToString[allenergies[[#, i]]]], Spacer[3], 
ListPlot[{allenergies[[#]], {{i, allenergies[[#, i]]}}}, PlotRange -> All, ImageSize -> 200, PlotStyle -> PointSize[Large],
AxesLabel -> {"Step"}]}] & /@ Range[Length[patterns]], 3], Frame -> All], {i, 1, stepNumber, 1},SaveDefinitions -> True]

allenergies1 = Table[e[evolution1[[j, i]], j], {j, 1, Length[patterns]}, {i, 1, stepNumber}];

Manipulate[ Grid[Partition[Row[{ArrayPlot[evolution1[[#, i]], PixelConstrained -> 2, 
PlotLabel -> "E = " <> ToString[allenergies1[[#, i]]]], Spacer[3], ListPlot[{allenergies1[[#]], 
{{i, allenergies1[[#, i]]}}}, PlotRange -> All, ImageSize -> 200, PlotStyle -> PointSize[Large]]}] & /@
Range[Length[patterns]],3], Frame -> All], {i, 1, stepNumber, 1}, SaveDefinitions -> True]

CellularAutomaton

However the second rule - to flip is pathological and increases the energy globally.

One can also find out the number of the cellular automaton built. It is the binary code of the rule in the reverse order

rulNum = FromDigits[Reverse@tuout[[1, All, -1]], 2]

and draw the CA rules

RulePlot[CellularAutomaton[{rulNum, 2, {1, 1}}], Appearance -> "Simplified"]

energyDecreasingCArules

Here the dots mean any color.

Exhaustive search 3x3 field

Since there are just 2^9 possible configurations one can build them all at once and calculate their energy

patterns = {{{0, 1, 0},{1, 0, 1},{0, 1, 0}},  -{{0, 1, 0},{1, 0, 1},{0, 1, 0}},{{0, 1, 0},{2, 0, 2},{0, 1, 0}},
-{{0, 1, 0},{2, 0, 2},{0, 1, 0}},{{0, 1, 0},{-1, 0, -1},{0, 1, 0}}, {{0, 1, 0},{-2, 0, -2},{0, 1, 0}}), 
-{{0, 1, 0},{-2, 0, -2},{0, 1, 0}}, {{0, 1, 0},{0, 0, 0},{0, 1, 0}},-{{0, 1, 0},{0, 0, 0},{0, 1, 0}}};

lattice = Tuples[{-1, 1}, {#, #}] &[3];

e[s_, patternNum_] := Apply[Plus, s ListConvolve[patterns[[patternNum]], s, 2], {0, 1}]

To take into account periodic boundary conditions and the symmety with respect to the flip of all spins, one has to exclude extra configurations

factorOutRotations[minE1_, size_] :=  Block[{minE = minE1, n = 1}, While[n < Length[minE],
With[{equivalentsets = Table[RotateLeft[RotateLeft[Partition[IntegerDigits[minE[[n]], 2, size^2]
, size], k]\[Transpose], j]\[Transpose], {k, 0, size - 1}, {j, 0, size - 1}]}, 

minE = Cases[minE,Except[Alternatives @@ Rest@Union@Flatten@Map[FromDigits[Flatten[#], 2] &, 
        Join[equivalentsets, equivalentsets /. {0 -> 1, 1 -> 0}], {2}]]]]; n++;]; minE]

exhausiveEnergies = Table[e[#, j] & /@ lattice, {j, 1, Length[patterns]}];
minEnergies = Min /@ exhausiveEnergies

energyMinimumPositions = MapThread[Position][{exhausiveEnergies, minEnergies}];
res4 = Map[FromDigits[Flatten[lattice[[#]] /. -1 -> 0], 2] &, energyMinimumPositions, {3}];
groundState = factorOutRotations[#, 3] & /@ (Flatten[res4, {3}][[1]])
 Row /@ Map[ArrayPlot[Partition[IntegerDigits[#, 2, 9], 3], PixelConstrained -> 10] &, groundState, {2}]

3x3

For bigger fields one can implement sequential algorithm

energyCompiled5 = 
 Compile[{{fieldN, _Integer}, {w1, _Integer}, {w2, _Integer}, {w3, _Integer}, {w4, _Integer}}, 
  Block[{d = Partition[N[IntegerDigits[fieldN, 2, 25]] /. 0. -> -1., 5]}, 
   w3 (d[[1, 2]] d[[2, 1]] + d[[3, 5]] d[[2, 1]] + 
   d[[1, 3]] d[[2, 2]] + d[[1, 4]] d[[2, 3]] + 
   d[[1, 5]] d[[2, 4]] + d[[1, 1]] d[[2, 5]] + 
   d[[2, 2]] d[[3, 1]] + d[[2, 3]] d[[3, 2]] + 
   d[[2, 4]] d[[3, 3]] + d[[2, 5]] d[[3, 4]] + 
   d[[3, 2]] d[[4, 1]] + d[[3, 3]] d[[4, 2]] + 
   d[[3, 4]] d[[4, 3]] + d[[3, 5]] d[[4, 4]] + 
   d[[3, 1]] d[[4, 5]] + d[[1, 5]] d[[5, 1]] + 
   d[[4, 2]] d[[5, 1]] + d[[1, 1]] d[[5, 2]] + 
   d[[4, 3]] d[[5, 2]] + d[[1, 2]] d[[5, 3]] + 
   d[[4, 4]] d[[5, 3]] + d[[1, 3]] d[[5, 4]] + 
   d[[4, 5]] d[[5, 4]] + d[[1, 4]] d[[5, 5]] + 
   d[[4, 1]] d[[5, 5]]) + 
w4 (d[[1, 5]] d[[2, 1]] + d[[3, 2]] d[[2, 1]] + 
   d[[1, 1]] d[[2, 2]] + d[[1, 2]] d[[2, 3]] + 
   d[[1, 3]] d[[2, 4]] + d[[1, 4]] d[[2, 5]] + 
   d[[2, 5]] d[[3, 1]] + d[[2, 2]] d[[3, 3]] + 
   d[[2, 3]] d[[3, 4]] + d[[2, 4]] d[[3, 5]] + 
   d[[3, 5]] d[[4, 1]] + d[[3, 1]] d[[4, 2]] + 
   d[[3, 2]] d[[4, 3]] + d[[3, 3]] d[[4, 4]] + 
   d[[3, 4]] d[[4, 5]] + d[[1, 2]] d[[5, 1]] + 
   d[[4, 5]] d[[5, 1]] + d[[1, 3]] d[[5, 2]] + 
   d[[4, 1]] d[[5, 2]] + d[[1, 4]] d[[5, 3]] + 
   d[[4, 2]] d[[5, 3]] + d[[1, 5]] d[[5, 4]] + 
   d[[4, 3]] d[[5, 4]] + d[[1, 1]] d[[5, 5]] + 
   d[[4, 4]] d[[5, 5]]) + 
w1 (d[[1, 1]] d[[2, 1]] + d[[3, 1]] d[[2, 1]] + 
   d[[1, 2]] d[[2, 2]] + d[[1, 3]] d[[2, 3]] + 
   d[[1, 4]] d[[2, 4]] + d[[1, 5]] d[[2, 5]] + 
   d[[2, 2]] d[[3, 2]] + d[[2, 3]] d[[3, 3]] + 
   d[[2, 4]] d[[3, 4]] + d[[2, 5]] d[[3, 5]] + 
   d[[3, 1]] d[[4, 1]] + d[[3, 2]] d[[4, 2]] + 
   d[[3, 3]] d[[4, 3]] + d[[3, 4]] d[[4, 4]] + 
   d[[3, 5]] d[[4, 5]] + d[[1, 1]] d[[5, 1]] + 
   d[[4, 1]] d[[5, 1]] + d[[1, 2]] d[[5, 2]] + 
   d[[4, 2]] d[[5, 2]] + d[[1, 3]] d[[5, 3]] + 
   d[[4, 3]] d[[5, 3]] + d[[1, 4]] d[[5, 4]] + 
   d[[4, 4]] d[[5, 4]] + d[[1, 5]] d[[5, 5]] + 
   d[[4, 5]] d[[5, 5]]) + 
w2 (d[[1, 1]] d[[1, 2]] + d[[1, 3]] d[[1, 2]] + 
   d[[1, 3]] d[[1, 4]] + d[[1, 1]] d[[1, 5]] + 
   d[[1, 4]] d[[1, 5]] + d[[2, 1]] d[[2, 2]] + 
   d[[2, 2]] d[[2, 3]] + d[[2, 3]] d[[2, 4]] + 
   d[[2, 1]] d[[2, 5]] + d[[2, 4]] d[[2, 5]] + 
   d[[3, 1]] d[[3, 2]] + d[[3, 2]] d[[3, 3]] + 
   d[[3, 3]] d[[3, 4]] + d[[3, 1]] d[[3, 5]] + 
   d[[3, 4]] d[[3, 5]] + d[[4, 1]] d[[4, 2]] + 
   d[[4, 2]] d[[4, 3]] + d[[4, 3]] d[[4, 4]] + 
   d[[4, 1]] d[[4, 5]] + d[[4, 4]] d[[4, 5]] + 
   d[[5, 1]] d[[5, 2]] + d[[5, 2]] d[[5, 3]] + 
   d[[5, 3]] d[[5, 4]] + d[[5, 1]] d[[5, 5]] + 
   d[[5, 4]] d[[5, 5]])], CompilationTarget -> "C", 
  CompilationOptions -> {"InlineExternalDefinitions" -> True},
  RuntimeAttributes -> {Listable}, Parallelization -> True]

res = Table[If[w2 w1 != 0, {w1, w2, 0, 0} -> factorOutRotations[Flatten[Position[#, Min[#]] - 1], 5] &
[ParallelMap[energyCompiled5[#, w1, w2, 0, 0] &, Range[0, 2^25 - 1], {1}]]], {w1, -2, 2}, {w2, -2, 2}]

res1 = Column@DeleteCases[Flatten@res, Null] /. (a_ -> b_) :> a -> Row[ (ArrayPlot[Partition[IntegerDigits[#, 2, 25], 5], 
     PixelConstrained -> 10] & /@ b)]

5x5

POSTED BY: Andrey Grabovskiy
Answer
2 months ago

Group Abstract Group Abstract