Monte Carlo with Metropolis algorithm extremely slow in Python

Question

I'm trying to implement a simple Monte Carlo in Python (to which I'm fairly new). Coming from C I'm probably following the wrongest path since my code is far too slow for what I'm asking: I have a potential hard sphere-like (see V_pot(r) in the code) for 60 3d particles and periodic boundary conditions (PBC), so I defined the following functions

import timeit
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from numpy import inf
#
L, kb, d, eps, DIM = 100, 1, 1, 1, 3
r_c, T = L/2, eps/(.5*kb)
beta = 1/(kb*T)
#
def dist(A, B):
    d = A - B
    d -= L*np.around(d/L)
    return np.sqrt(np.sum(d**2))

#
def V_pot(r):
    V = -eps*(d**6/r**6 - d**6/r_c**6)
    if r > r_c:
        V = 0
    elif r < d:
        V = inf

    return V

#
def ener(config):
    V_jk_val, j = 0, N
    #
    while (j > 0):
        j -= 1
        i = 0
        while (i < j):
            V_jk_val += V_pot(dist(config[j,:], config[i,:]))
            i += 1

    #
    return V_jk_val

#
def acc(en_n, en_o):
    d_en = en_n-en_o
    if (d_en <= 0):
        acc_val = 1
    else:
        acc_val = np.exp(-beta*(d_en))

    return acc_val

#

then, starting from the configuration (where every line of the array represents the coordinates of a 3D particle)

config = np.array([[16.24155657, 57.41672173, 94.39565792],
       [76.38121764, 55.88334066,  5.72255163],
       [38.41393783, 58.09432145,  6.26448054],
       [86.44286438, 61.37100899, 91.97737383],
       [37.7315366 , 44.52697269, 23.86320444],
       [ 0.59231801, 39.20183376, 89.63974115],
       [38.00998141,  3.84363202, 52.74021401],
       [99.53480756, 69.97688928, 21.43528924],
       [49.62030291, 93.60889503, 15.73723259],
       [54.49195524,  0.6431965 , 25.37401196],
       [33.82527814, 25.37776021, 67.4320553 ],
       [64.61952893, 46.8407798 ,  4.93960443],
       [60.47322732, 16.48140136, 33.26481306],
       [19.71667792, 46.56999616, 35.61044526],
       [ 5.33252557,  4.44393836, 60.55759256],
       [44.95897856,  7.81728046, 10.26000715],
       [86.5548395 , 49.74079452,  4.80480133],
       [52.47965686, 42.831448  , 22.03890639],
       [ 2.88752006, 59.84605062, 22.75760029],
       [ 9.49231045, 42.08653603, 40.63380097],
       [13.90093641, 74.40377984, 32.62917915],
       [97.44839233, 90.47695772, 91.60794836],
       [51.29501624, 27.03796277, 57.09525454],
       [10.30180295, 21.977336  , 69.54173272],
       [59.61327648, 14.29582325, 11.70942289],
       [89.52722796, 26.87758644, 76.34934637],
       [82.03736088, 78.5665713 , 23.23587395],
       [79.77571695, 66.140968  , 53.6784269 ],
       [82.86070472, 40.82189833, 51.48739072],
       [99.05647523, 98.63386809,  6.33888993],
       [31.02997123, 66.99709163, 95.88332332],
       [97.71654767, 59.24793618,  5.20183793],
       [ 6.79964473, 45.01258652, 48.69477807],
       [93.34977049, 55.20537774, 82.35693526],
       [17.35577815, 20.45936211, 29.27981422],
       [55.51942207, 52.22875901,  3.6616131 ],
       [61.45612224, 36.50170405, 62.89796773],
       [23.55822368,  7.09069623, 37.38274914],
       [39.57082799, 58.95457592, 48.0304924 ],
       [93.94997617, 64.34383203, 77.63346308],
       [17.47989107, 90.01113402, 81.00648645],
       [86.79068539, 66.35768515, 56.64402907],
       [98.71924121, 38.33749023, 73.4715132 ],
       [ 0.42356139, 78.32172925, 15.19883322],
       [77.75572529,  2.60088767, 56.4683935 ],
       [49.76486142,  3.01800153, 93.48019286],
       [42.54483899,  4.27174457,  4.38942325],
       [66.75777178, 41.1220603 , 19.64484167],
       [19.69520773, 41.09230171,  2.51986091],
       [73.20493772, 73.16590392, 99.19174281],
       [94.16756184, 72.77653334, 10.32128552],
       [29.95281655, 27.58596604, 85.12791195],
       [ 2.44803886, 32.82333962, 41.6654683 ],
       [23.9665915 , 49.94906612, 37.42701059],
       [30.40282934, 39.63854309, 47.16572743],
       [56.04809276, 30.19705527, 29.15729635],
       [ 2.50566522, 70.37965564, 16.78016719],
       [28.39713572,  4.04948368, 27.72615789],
       [26.11873563, 41.49557167, 14.38703697],
       [81.91731981, 12.10514972, 12.03083427]])

I make the 5000 time steps of the simulation with the following code

N = 60
TIME_MC = 5000
DELTA_LIST = [d]
#d/6, d/3, d, 2*d, 3*d 
np.random.seed(19680801)
en_mc_delta = np.zeros((TIME_MC, len(DELTA_LIST)))
start = timeit.default_timer()
config_tmp = config
#
for iD, Delta in enumerate(DELTA_LIST):
    t=0

    while (t < TIME_MC):
        for k in range(N):
            RND = np.random.rand()
            config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)
            en_o, en_n = ener(config), ener(config_tmp)
            ACC = acc(en_n, en_o)
            if (RND < ACC):
                config[k,:] = config_tmp[k,:]
                en_o = en_n


        en_mc_delta[t][iD] = en_o
        t += 1


stop = timeit.default_timer()
print('Time: ', stop-start)

following the rule of the Metropolis algorithm for the acceptance of the proposed move extracted with config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5).

I made some attempts to check where the code get stuck and I found that the function ener (also because of the function dist) is extremely slow: it takes something like ~0.02s to calculate the energy of a configuration, which means something around ~6000s to run the complete simulation (60 particles, 5000 proposed moves).

The outer for it's just to calculate the results for different values of Delta.

Running this code with TIME_MC=60 can make you an idea of how much slow is this code (~218s) which takes just some seconds if implemented in C. I read some other question about how to speed up Python codes but I can't understand how to do it here.

EDIT:

I'm now almost sure that the problem is in the function dist, since just to calculate PBC distance between two 3D vectors it takes around ~0.0012s which gives crazy long times when you calculate it 5000*60 times.

Answer 1

Note that this is a partial answer continued from comments on the original question.

Here's an example of how "unrolling" numpy's function can improve performance when replaced with a more direct calculation of the distance. Note that this was not verified to be equivalent, especially concerning the rounding. The principle still applies, I think.

import random
import time
import numpy as np

L = 100
inv_L = 0.01
vec_length = 10
repetitions = 100000

def dist_np(A, B):
    d = A - B
    d -= L*np.around(d/L)
    return np.sqrt(np.sum(d**2))

def dist_direct(A, B):
    sum = 0
    for i in range(0, len(A)):
        diff = (A[0,i] - B[0,i])
        diff -= L * int(diff * inv_L)
        sum += diff * diff
    return np.sqrt(sum)

vec1 = np.zeros((1,vec_length))
vec2 = np.zeros((1,vec_length))

for i in range(0, vec_length):
    vec1[0,i] = random.random()
    vec2[0,i] = random.random()

print("with numpy method:")
start = time.time()
for i in range(0, repetitions):
    dist_np(vec1, vec2)
print("done in {}".format(time.time() - start))

print("with direct method:")
start = time.time()
for i in range(0, repetitions):
    dist_direct(vec1, vec2)
print("done in {}".format(time.time() - start))

Output:

with numpy method:
done in 6.332799911499023
with direct method:
done in 1.0938000679016113

Play around with the average vector length and the repetitions to see where the sweet spot is. I expect the performance gain is not constant when varying these meta-parameters.