I have been assigned to use the metropolis algorithm in order to solve numerically an integral. What am I doing wrong?

Question

I have been trying to solve the following for a couple of days now, using the Metropolis algorithm:

This is taken from: https://arxiv.org/pdf/astro-ph/0505561.

I basically approximate A with the Metropolis algorithm, by randomly placing a random particle somewhere inside the box. The following sum (m is the number of iterations for the Monte Carlo) is the one I use in my code:

Am I using the Metropolis algorithm correctly, or am I doing something wrong?

Any help would be greatly appreciated, since I have been stuck on it for a couple of days now :)

My code is the following :

@jit(nopython=True)
def potential(coord_array, A, N):
    energy_sum = 0
    for l in range(N-1):
        for j in range(l+1, N):
            dif = np.sqrt((coord_array[l, 0] - coord_array[j, 0])**2 +
                          (coord_array[l, 1] - coord_array[j, 1])**2 +
                          (coord_array[l, 2] - coord_array[j, 2])**2)
            if dif >= A:
                energy = -1/dif
            else:
                energy = 1/A
            energy_sum += energy
    return energy_sum / N

@jit(nopython=True)
def potential_move(coord_array, random_particle, A, N):
    energy_sum = 0
    for l in range(N):
        if l == random_particle:
            continue
        dif = np.sqrt((coord_array[l, 0] - coord_array[random_particle, 0])**2 +
                      (coord_array[l, 1] - coord_array[random_particle, 1])**2 +
                      (coord_array[l, 2] - coord_array[random_particle, 2])**2)
        if dif >= A:
            energy = -1 / dif
        else:
            energy = 1 / A
        energy_sum += energy
    return energy_sum / N



@jit(nopython=True)
def func(eta, N, A, coord_array):
    return np.exp(-eta * potential(coord_array, A, N))

@jit(nopython=True)
def func_move(eta, N, A, coord_array, final_coord, random_particle):
    return np.exp(-eta * (potential_move(final_coord, random_particle, A, N)-potential_move(coord_array, random_particle, A, N)))


'''
@jit(nopython=True)
def canonical_weight(eta, N, A, initial_coord, final_coord, random_particle):
    initial_func = func_move(eta, N, A, initial_coord, random_particle)
    final_func = func_move(eta, N, A, final_coord, random_particle)
    return final_func/initial_func
'''

@jit(nopython=True)
def canonical_ensemble(coord_array, A, N, eta, iterations, delta, m):
    final_coords = np.zeros((N, 3, iterations))
    acceptance_count = 0

    for i in range(iterations):
        random_particle = random.randint(0, N-1)
        if i == 0:
            coordinates = np.copy(coord_array)
        else:
            coordinates = np.copy(final_coords[:, :, i-1])

        for j in range(3):
            coordinates[random_particle, j] = np.random.random()#+= delta * (2 * np.random.random() - 1)#= np.random.random()#

        if i == 0:
            stat_weight = func_move(eta, N, A, coord_array, coordinates, random_particle)
        else:
            stat_weight = func_move(eta, N, A, final_coords[:, :, i-1], coordinates, random_particle)

        if stat_weight >=1: #(potential_move(coordinates, random_particle, A, N)-potential_move(coord_array, random_particle, A, N))<=0:
            final_coords[:, :, i] = coordinates
            acceptance_count += 1
        else:
            Omega = np.random.random()
            if min(1, stat_weight) >= Omega:
                final_coords[:, :, i] = coordinates
                acceptance_count += 1
            else:
                final_coords[:, :, i] = final_coords[:, :, i-1]

    integral_sum = 0
    for i in range(m,iterations):
        distance = np.sqrt((final_coords[0,0,i] - final_coords[1,0,i])**2 +
                           (final_coords[0,1,i] - final_coords[1,1,i])**2 +
                           (final_coords[0,2,i] - final_coords[1,2,i])**2)
        if distance >= A:
            integral_sum += (1 / distance)
        else:
            integral_sum += -(1/A)
    Phi =1 - eta * (N-1) * integral_sum / (2* 3 * N * (iterations-m))

    acceptance_rate = acceptance_count / iterations * 100
    return Phi, acceptance_rate, final_coords[:,:,-1]

N = 500
m= 5000
A = 10**(-4)-10**(-8)
initial_coords = np.random.rand(N,3)
#print(initial_coords)
delta = 0.1
iterations = 200000

eta_arr = np.arange(1.3,1.6,0.01)
PVNT = np.zeros((len(eta_arr)))
acceptance_rates = np.zeros((len(eta_arr)))

for i in range(len(eta_arr)):
    PVNT[i]= canonical_ensemble(initial_coords, A, N, eta_arr[i], iterations, delta, m)[0]
    acceptance_rates[i] = canonical_ensemble(initial_coords, A, N, eta_arr[i], iterations, delta, m)[1]
    print(f"eta: {eta_arr[i]}, Acceptance rate: {acceptance_rates[i]:.2f}%")

plt.plot(eta_arr, PVNT, marker='o')
plt.xlabel('Eta')
plt.ylabel('PV/NT')
plt.title(r'PV/NT vs $\eta$')
plt.ylim(0,1)
plt.xlim(1.3,1.6)
plt.show()

plt.plot(eta_arr, acceptance_rates, marker='o', color='red')
plt.xlabel('Eta')
plt.ylabel('Acceptance Rate (%)')
plt.title(r'Acceptance Rate vs $\eta$')
plt.show()

The expectation is Figure 1 in the paper, but i got nothing like it