 Solution

INITIALIZE(N, sigma):
    for i = 1,...,N:
        r[i] = ((i-1)*sigma, 0, 0)
    return r

ENERGY(r):
    for each pair i < j:
        if norm(r[i]-r[j]) < sigma:
            return infinity
    U = 0
    for i = 1,...,N-1:
        U += k/2 * (norm(r[i+1]-r[i])-sigma)^2
    return U

TRIAL(r):
    choose i uniformly from 1,...,N
    choose delta uniformly from [-a,a]^3
    r_trial = copy(r)
    r_trial[i] += delta
    return r_trial

ACCEPT(U, U_trial, beta):
    if U_trial == infinity:
        return false
    return log(uniform(0,1)) < min(0, -beta*(U_trial-U))

MAIN:
    r = INITIALIZE(N, sigma)
    U = ENERGY(r)
    for sweep = 1,...,n_equilibrate+n_produce:
        repeat N times:
            r_trial = TRIAL(r)
            U_trial = ENERGY(r_trial)
            if ACCEPT(U, U_trial, 1/(k_B*T)):
                r = r_trial
                U = U_trial
        if sweep > n_equilibrate:
            d_ETE = norm(r[N]-r[1])
            store d_ETE and update its histogram
    estimate autocorrelation and block uncertainties