Алгоритм Метрополис: цикл не повторяется в заданном диапазоне

Question

Я нахожусь в процессе написания некоторого кода в попытке использовать алгоритм Метрополиса для аппроксимации энергий основного состояния и других свойств квантово-механических систем.Итак, сначала я определил функцию delta_action, а также ряд параметров следующим образом:

import numpy.random as rnd;
import numpy as np;
import math;
#Definition of variables

EL = 10
N = 10              #number of lattice points
A = EL/N            #lattice constant
M = 1000000         #number of iterations
ME = 500            #burn in iterations
MA = 5              #every MA'th configuration is measured
BIN = 40            #number of bins for the wavefunction
INTERV = 2          #interval for binning [-INTERV,INTERV]
MASS = 1.0
MU = .5             #coefficient of q**2
LAMBDA = 0.0        #coefficient of q**4
DELTA = 0.5         #change of lattice variable paramater
massl = MASS/A
lambdal = A*LAMBDA
muleff = MASS/A + A*MU
reject = 0
translate = BIN/2.0
stretch = 0.5*BIN/INTERV
q = np.zeros(N)


def del_action(x,y,xs):
    return (y-x)*((y+x)*(muleff+lambdal*(y*y+x*x))-massl*xs)

Теперь у меня есть следующее:

def metropolis(q):
    for i in range(0,MA):
        for j in range(0,N-1):
            jp = ((j+1)%N)
            jm = ((j+N-1)%N)
            qnew = q[j] + DELTA*(1-2*rnd.uniform(0,1))
            dS = del_action(qnew,q[j],q[jp]+q[jm])
            if dS<0:
                q[j] = qnew
            else:
                rn = rnd.uniform(0,1)
                if rn <= math.exp(-dS):
                    q[j] = qnew
                return

Итак, у меня проблема в том, чтоцикл по j постоянно застревает на 0,1,2 и не достигает N-1.Добавив print(j,jp,jm,qnew,dS) чуть выше if dS < 0: и запустив, я получаю следующий результат:

0 1 9 -0.10299220871627934 -0.01591109258438647
1 2 0 0.4040577654073293 -0.2865088183872116
2 3 1 0.22295515394986587 0.015523260281817304
0 1 9 0.07350728912179683 0.07912210294530213
0 1 9 0.4117795868978783 -0.10955711120902165
1 2 0 0.06283878213693184 0.022387404949800852
0 1 9 0.5864083527452271 -0.2504950330034522
1 2 0 0.16584582581090546 0.048035953984589264
0 1 9 0.20274281147744633 0.3905258343186439
0 1 9 -0.1002558575204866 -0.0036709485216505355
1 2 0 -0.03319832019662894 0.01518148753902383
0 1 9 -0.42958241606820047 -0.25080163430181573
1 2 0 -0.2491635871862874 -0.04684623524276284
2 3 1 0.21831753974485746 0.00422520336014744
0 1 9 -0.533602476871921 -0.12436781520409396
1 2 0 -0.49950732570395673 -0.2022080030198727
2 3 1 0.4907960348437076 -0.425932303885753
3 4 2 0.24493839231060366 0.030222567685481976
0 1 9 -0.20701571624479698 0.19968166541608104

Где столбцы j, jp, jm, qnew, dS соответственно.Jp и jm относятся к соседним точкам каждого значения решетки, которые необходимы для функции delta_action.

Итак, мой вопрос: почему j застревает внутри и около 0, как вы можете видеть сверхув конечном итоге он оказывается равным нулю гораздо чаще, чем нет, после тысяч итераций может быть блок из 10 j = 0.

Любая помощь будет высоко ценится.Я также хотел бы добавить, что я совершенно новичок в Python и программировании в целом, я на последнем курсе степени математической науки, но полностью пренебрегал любым аспектом компьютерного программирования или кодирования, предлагаемым моим курсом на сегодняшний день, и яЯ плачу за это сейчас.Спасибо!