MC simulation of polymer in Python

Suppose, a polymer has N monomers in its chain. I want to simulate its movement using the bead-spring model.

So, I wrote the following program in Python -

• GitHub : Polymer chain simulation with Monte Carlo method

Here is the main simulation part:

def run_simulation(steps_int):
    for i in range(steps_int):
        rand_index_int = random.randint(0, len(polymer_chain_vec)-1)  # choose a random index
        before_loc_pt = polymer_chain_vec[rand_index_int]  # obtain the bead location at the index
        before_pot_float = get_potential(rand_index_int)  # get the bead potential
        new_loc_pt = get_point_at_radius(before_loc_pt, max_atom_distance_float)  # get a random location
        polymer_chain_vec[rand_index_int] = new_loc_pt  # move the bead to the new point
        after_pot_float = get_potential(rand_index_int)  # get the bead potential again

        pot_diff_float = after_pot_float - before_pot_float

        if pot_diff_float < 0:  # e_before > e_after
            pass
        else:
            rand_float = random.random()
            if math.exp(-pot_diff_float / temperature_float) > rand_float:
                pass
            else:
                polymer_chain_vec[rand_index_int] = before_loc_pt  # restore the old location
                # print("C", end="")

I am using 100000 steps. The X-axis is steps. Y-axis is total energy.

While using LJ potential function, the energy is not fluctuating as expected.

I am especially concerned with the function that calculates the total energy of the polymer. Coz, the value of total energy is unexpectedly large.

Can anyone check the source code and tell me what I am doing incorrectly?

from MC simulation of polymer in Python