Problem concerning a part of a simulation for the Ising Model

Question

I was wondering if anyone could help me with a simulation in Python (there is a small part of the code where I have a question) about the Ising model. First, I wanted to get the magnetization $M$ which is defined as: $$ M = \sum_i \sigma ^ z_i $$

And the part of the code to obtain $M$ was successful and was the following thing :

sigmaxop = []  
sites = []  
for i in range (L):  
    # \ sum_i sigma ^ z (i)  
    sigmaxop.append ((sigmaz) .ToList ())  
    sites.append ([i])  

where L is the length of the network (ie the number of spins).

But now I would like to get the average fluctuations ie:
$$ \langle M^2\rangle = \sum_ {i, j} \langle\sigma_i \sigma_j\rangle $$

And to get this average fluctuation, I tried to rewrite/modify the part of the code above (that allowed me to have $M$) to obtain instead the fluctuations $\langle M^2\rangle$ but I am wrong every time I try something: I do not get these fluctuations $\langle M^2\rangle$. I feel lost.

Someone (paradoxy) gave me some welcome informations to help: it was that I have to note that $m_j=\sum_k \sigma_k$ (for every lattice) and $<M>=\sum_j m_j$ which is the average of magnetization. And so $<M^2>=\sum_j m^2_j$. But finally after trying again I am still stuck with my code. I didn't succeed yet.

Would anyone have any idea how I could modify this part of code to get the fluctuations $\langle M^2\rangle$, please?

Answer 1

It is very difficult to actually help you if I can't see how you are calculating anything... but I can try to give you some advice. I think you are on the right track. This is a great problem to learn how Monte-Carlo simulations work, especially for more complicated QFT calculations. There is a lot of good documentation I suggest you search through to familiarize yourself. I don't know how well you know the Ising Model, but I assume your are learning this for the first time. Before talking about $\langle M \rangle$ and $ \langle M^2 \rangle $ and their implementations, make sure you know the physical system you are trying to program! You should be able to tell me why I should care about the observables you are trying to implement. Wikipedia/the internet is your friend here. Here is a nice resource. After you do your research, come back to your code and draft an outline of the problem with your notes.

Some helpful things to consider, I don't know if you are doing these, but it may help identify issues.

1. Appending to lists works, but I suggest you familiarize yourself with the numpy array. It is much clearer and guaranteed to be faster. Look at any sophisticated numerical code still written in python and you will see it. Declaring zero arrays with np.zeros() at the start, diagnosing problems down the line is much easier as you will immediately notice indexing errors due to size mismatches.
2. Make sure to make your program modular, that is, take advantage of being able to define independent functions. It will help keep variables under control and make it easy to test new solutions. This means you should have a function to do your measurements (in this case $\langle M \rangle$ and $ \langle M^2 \rangle $ but there can be an entire battery of complicated observables here!). I would try to measure $\langle E \rangle, \langle E^2 \rangle, \langle M \rangle, \langle M^2 \rangle, $and$ \langle |M| \rangle$.

In the case of these observables you mentioned, you are close. You are confusing magnetization and energy though, and your notation is incorrect. These observables are a single number we get from a whole field configuration. You are correct that each lattice point has a single value, and then we sum all of the lattice points. If we want to compare different lattice sizes though, we are going to have problems with comparing them, since 10 sites with spin +1 will give 10, but 20 sites with spin +1 will give 20. It is common to measure the 'Observable-Per-Site'. This idea applies to more complex calculations and observables as well. Let's clear up some things.

• The Ising Hamiltonian without a background field and constant coupling term J looks like $H(\sigma)=J \sum_{\langle i,j\rangle} \sigma_i \sigma_j $. What you defined as $\langle M \rangle$ is actually $H$, no brackets, which is the energy of a configuration. I usually define in $d$ dimensions the Energy-Per-Site as $E_s = \frac{1}{L^d} \sum_{\langle i,j\rangle} J_{i j} \sigma_i \sigma_j $ Which assuming a constant $J$ in 1D can be written as $E_s =\frac{J}{L^d} \sum_i ( \sigma_{i} \sigma_{i+1} + \sigma_{i} \sigma_{i-1} )$. There is sometimes a factor of 2 in the denominator as well, I don't have it in my notes, but it is present in some sources but not others. Energy is a good measurement to make, because it is the function we are trying to minimize, so it's good for visualization, but also used in calculating heat capacity for example.
• Taking advantage of numpy array operations you can usei=np.arange(0,L) and then use this array as an index, Sigma[i]*(Sigma[i+1]+Sigma[i-1]). In 2D (which you should find slightly more interesting from a physical perspective, the np.meshgrid() function can be used in a similar way. I encourage you to experiment with this.
• Now, here is the heart of the analysis. This may be generalized, but lets focus on magnetization, and keep it simple. The magnetization per site $M$ is defined as $M=\frac{1}{L^d}\sum_i \sigma_i$ is what we measure as the program runs, but at the end of the program (or you can do this as a "running average" if you want) we calculate the expectation value $\langle M \rangle = \frac{1}{N}\sum_n M_n$. Also, you should see that it takes some time to approach the mean value, this is typically called "burn in" time or "thermalization" in the QFT context, and for now you want to disregard those values in your analysis. If you were to run a more complex simulation, you would have to account for the correlation time between measurements and being careful about tuning these parameters. Additionally, your measurements need error constraints. The simplest way to do this is the standard deviation perhaps, but for more complex simulations you might use the jackknife method.

Now more specifically, you should ask why we also want a magnitude like $\langle M^2 \rangle$ rather than the $\pm$ number that is $\langle M \rangle$. This should be quite simply implemented by just squaring your array M**2. These however are defined straightforwardly as $ M = \sum_i \sigma_i$ and $ M^2 = \sum_i \sigma_i^2$

The heart of my Ising model code is simply the following, I use S as the lattice of spins, a numpy ndarray. this can be any dimension here.

def run_model(S,Energy,maxit,Evec,Evec2,Mag,Mag2.B): for i in range(maxit): S=bnc_periodic(S) S,Energy=markov_ising(S,Energy,B) Evec[i] = Energy Evec2[i] = Energy**2 Mag[i] = np.sum(S) Mag2[i] = np.sum(S**2) return(M,Evec,Evec2,Mag,Mag2)

My function markov_ising() is just the MCMC step and calculates the energy as I described. In fact, if you are clever, you can see that there is no need to recalculate the energy over the whole field if you just change one point (This is very important in larger simulations as it takes a very long time to measure the energy). So you can see that it's a relatively quick and clean measurement. This kind of coding makes it much easier to extend to higher dimensions, or even other models. And debugging is straightforward. Then I calculate expectation values in an analysis function, by taking the mean of the measured values.
I pulled most of this from my notes, but it is identical to what is on wikipedia and the first few results of "Ising model" in google. G## Heading ##ood luck!