1

I have been working on a Monte Carlo simulation on the two-dimensional Ising model for quite a while now. The problem I have is purely computational and does not require any knowledge about the Ising model.

I have been able to implement the model completely. However, the code is not running "fast enough" as I want to simulate quite large systems. Running a suitable number of Monte Carlo sweeps, the run time is around half an hour using the largest system I have tried.

The main reason for the high runtime is that I'm not able to keep the functions running on a "numpy form", I have to use an ordinary Python for-loop which, for me, seems as the largest problem. The two functions I need help with are the following.

def helix(S, i, N, L, beta, r, B = 0, J_v = 1, J_p = 1):
    nn = i + 1
    if nn >= N:
        nn -= N
    sum_nn = S[nn]
    nn = i - 1
    if nn < 0:
        nn += N
    sum_nn += S[nn]
    nn = i + L
    if nn >= N:
        nn -= N
    sum_nn += S[nn]
    nn = i - L
    if nn < 0:
        nn += N
    sum_nn += S[nn]
    del_E = 2 * S[i] * sum_nn
    if del_E <= 0:
        S[i] = - S[i]
        return del_E, 2 * S[i] / N
    elif np.exp(-beta * sum_nn) > r:
        S[i] = - S[i]
        return del_E, 2 * S[i] / N
    else:
        return 0, 0


def random_sweep(S, N, L, beta, B = 0, J_v = 1, J_p = 1):
    del_m_sweep = 0
    del_He_sweep = 0
    rand_list_index = np.random.randint(N, size=N)
    rand_numbers = np.random.rand(N)
    for i in range(N):
        del_E, del_m = helix(S, rand_list_index[i], N, L, beta, 
                              rand_numbers[i], B, J_v, J_p)
        del_He_sweep += del_E
        del_m_sweep += del_m
    return del_He_sweep, del_m_sweep

I want to keep everything on numpy form as it is my understanding that this will help tremendously with computational speed. Basically I want to replace the for loop in random_sweep by just calling the function helix(S, rand_list_index, N, L, beta, rand_numbers, B, J_v, J_p), and using the numpy flexibility, it will automatically call the helix function N times (based on the arguments rand_list_index and rand_numbers instead of rand_list_index[i] and rand_numbers[i]). The reason I'm not able to this as of now is that I have to update the S-array at index i. Is there a way around this update?

I have been trying to implement this for quite a while now, so I would really appreciate any help!

rammelmueller
  • 1,092
  • 1
  • 14
  • 29
Sausage123
  • 11
  • 1
  • Why don't you code it in C or F, wrap it and call your function from python? It's not too challenging and a useful thing to learn how to do. But probably you just want to take advantage of vectorization in numpy - others will I'm sure help you here.. – jtlz2 Mar 24 '19 at 13:49
  • https://docs.scipy.org/doc/numpy/user/c-info.ufunc-tutorial.html – jtlz2 Mar 24 '19 at 13:50
  • Also Ising models were never much my thing but have you looked at how others have done it - https://rajeshrinet.github.io/blog/2014/ising-model? And: have you profiled your code to check this is what you need to optimize? – jtlz2 Mar 24 '19 at 13:53
  • 1
    Can you also post the equations? – jtlz2 Mar 24 '19 at 13:57
  • The equations used in the simulation are, for all practical purposes, the same as in the link you gave me. I only calculate the total energy and magnetization at the start of the simulation, afterwards I only calculate the change of energy a magnetization due to single flipped spin (del_E, and del_m). I will definetly try to use F or C if this can not be done in numpy, thank you for the tip! I just want to explore the possbilities of keeping all parts of the code in python as this slighly more convenient – Sausage123 Mar 24 '19 at 14:32
  • But this code is currently presuming that J_v = J_p = 1, and B =0 – Sausage123 Mar 24 '19 at 14:52
  • Also https://towardsdatascience.com/data-science-with-python-turn-your-conditional-loops-to-numpy-vectors-9484ff9c622e – jtlz2 Mar 24 '19 at 15:34
  • Could you format your code above ^^^^ so it's clearer where the if statements should be etc. - was going to try and code it up but I want to get it right. Or maybe github/gist? – jtlz2 Mar 25 '19 at 06:21
  • You could use Numba. For example: https://stackoverflow.com/a/45403017/4045774 – max9111 Mar 26 '19 at 14:00

0 Answers0