Ising model, proposed by Ernst Ising, desribes an ordered array of atoms with magnetic dipole moments called spins. Hamiltonian of such system is $$H =-\sum_{\left\langle i,j\right\rangle }J_{ij}\sigma_{i}\sigma_{j}-h\sum_{i}\sigma_{i} $$ where $\sigma_{i}=\pm1$ is an Ising spin variable, $h$ is a magnetic field, $J_{ij}$ is an interaction constant between spins and $\left\langle i,j\right\rangle $ denotes sum over nearest neighbours only. This model can be solved analyticaly in one and two dimensions. At zero magnetic field, in two dimensions and further, the system undergoes a phase transition from phase with finite non-zero magnetization to zero magnetization, as temperature is increased above critical value $T_c$.

Spin glass model based on Ising model in the way that the coupling constant $J_{ij}$ varies from site to site, introducing frustration to the lattice. Disorder in interaction makes this system to be more interesting and even harder to solve. Common propertie of glassy system are it's slow dynamics and complicated energy landscape. The system may be trapped in local valley of minimal energy. Specific model simulated here is Edwards-Anderson model, where interaction is drawn with equal probability to be +1 or -1. When $J_{ij}=+1$, spins tend to align and magnetzation bulds up in that direction. When $J_{ij}=-1$, energeticly favorable to anti-align the neighbouring spins with no total magnetization.

Using statistical mechanics approach, the probability to find the system in specific configuration $\left\{ \mu\right\} $ is defined by Boltzmann distribution: $$P_{\left\{ \mu\right\} }=\frac{1}{z}e^{-\beta H_{\left\{ \mu\right\} }}$$ where $z$ is the partition function $$z=\sum_{\left\{ \mu\right\} }e^{-\beta H_{\left\{ \mu\right\} }}$$ One can study the system by calculating any observables like, magnetization or energy per spin using the definition $$\left\langle m\right\rangle =\frac{\left\langle M\right\rangle }{N}=\frac{1}{N}\sum_{\left\{ \mu\right\} }M_{\left\{ \mu\right\} }P_{\left\{ \mu\right\} }$$ or $$\frac{\left\langle E\right\rangle }{N}=\frac{1}{N}\sum_{\left\{ \mu\right\} }E_{\left\{ \mu\right\} }P_{\left\{ \mu\right\} }$$ Here $$M_{\left\{ \mu\right\} }=\sum_{i\in\left\{ \mu\right\} }\sigma_{i}$$ sum over certain configuration of the lattice and $$E_{\left\{ \mu\right\} }=-\sum_{\left\langle i,j\right\rangle \in\left\{ \mu\right\} }J_{ij}\sigma_{i}\sigma_{j}- h\sum_{i\in\left\{ \mu\right\} }\sigma_{i}$$ energy of certain configuration. The problem, however, is that it is extremely difficult to calculate $z$. The size of configurations space is exponentially large with system size. That's where computer simulation gets in.

Approximation can be made in order to calculate an observable, say magnetization $$\left\langle m\right\rangle \approx\bar{m}=\sum_{i=1}^{N}M_{i}$$ where $M_{i}$ is a magnetization of some configuration, that was drawn from the Boltzmann distribution. This insured by law of large numbers, namely $N$ must be taken to be large ( e.g. $N>100000$). Method for obtaining a configuration from specific distribution used here is Markov Chain. Markov chain describes system with states that transform from one to another with some transition matrix. Those states are probability distributions of the system. After many transitions, the system will converge to a unique stationary distribution. This distribution is the Boltzmann distribution, and can be obtained when detailed balance condition is satisfied. Briefly, if $P_{1}$ and $P_{2}$ are some probability distributions, $\Pi_{12}$ and $\Pi_{21}$ are transition matrices for $P_{1}\rightarrow P_{2}$ and $P_{2}\rightarrow P_{1}$, respectively. Then detailed balance condition is $$P_{1}\Pi_{12}=P_{2}\Pi_{21}$$ From next equation can be seen that the partition function is canceled if the condition is satified $$\frac{\Pi_{12}}{\Pi_{21}}=\frac{P_{2}}{P_{1}}=\frac{\frac{1}{z}e^{-\beta E_{1}}}{\frac{1}{z}e^{-\beta E_{2}}}=e^{-\beta\left(E_{1}-E_{2}\right)}$$

To utilize this knowledge on how to exclude the partition function, we apply here the Metropolisâ€“Hastings algorithm that satisfies detailed balance condition and generates states according to Boltzmann coefficient.

1. Choose random spin on the lattice and suggest to flip it.

2. Calculate $A_{12}$ from $e^{-\beta\left(E_{1}-E_{2}\right)}$.

3. Draw u fom uniform [0,1] and check if u<$A_{12}$ then accept the move, flip the spin state. Else do not flip.

4. Measure observable (like magnetization).

5. Repeat from 1 for desired number of runs.

Concerning the demand of Markov chain process to many transitions before the system reaches stationary distribution, the simulation executes "burn runs" . Those runs are excluded from calculation of observables in order to give time to the system to reach the Boltzmann distribution.

Also in order to avoid correlations between measurements, simulation performs a certain number of Monte Carlo moves (MC move). Only after $L\times L\times L$ MC moves the observavle is measured.

Initial code was weritten by Dr. Snir Gazit for Statistical mechanis II course. The code was modified by me for 3D Sping glass.