Python script for Metropolis Monte Carlo (global optimisation)

Question

I have written a python script to do global optimisation on interatomic potential (IP) using Metropolis Monte Carlo.

I have attached the part of the script which is the Metropolis MC part.

I have plotted the plot: IP energy (no Coulomb interaction in the energy) as a function of step number (from the first accepted order to the last).

I expect the plot shows the negative trend data because I expect it converges as the step increases. However, it is not.

Is the plot represent the Markov chain?

# Initialize variables
energy_previous = 0.0
energy_vec = []

while step_cnt < max_steps:
    # Calculate energy of proposed configuration
    IP_energy = calculate_energy(proposed_configuration)

    if step_cnt == 0:
        energy_previous = IP_energy
        step_cnt = 1
    else:
        delta = IP_energy - energy_previous
        energy_tol = np.random.uniform()

        # Accept or reject the new configuration based on the Metropolis-Hastings criterion
        if energy_tol < np.exp(-delta):
            energy_vec.append(IP_energy)    # store the accepted structure energy
            energy_previous = IP_energy

        # update the proposed configuration for the next iteration
        proposed_configuration = update_configuration(proposed_configuration)
        step_cnt += 1

Here is the output of
step_no, IP_energy, previous_energy, delta, energy_tol, np.exp(-delta), energy_tol < np.exp(-delta)

 1   -0.00179377    -0.00967768     0.00788391  0.08959475  0.99214709  True
 2   17.20824264    -0.00179377     17.21003641     0.47674332  0.00000003  False
 2   -0.02482590    -0.00179377     -0.02303213     0.94119261  1.02329942  True
 3   -0.02276726    -0.02482590     0.00205864  0.46212012  0.99794348  True
 4   -0.00245595    -0.02276726     0.02031131  0.03979358  0.97989358  True
 5   -0.01868081    -0.00245595     -0.01622486     0.76820130  1.01635720  True
 6   -0.00820708    -0.01868081     0.01047373  0.36041987  0.98958093  True
 7   -0.01322191    -0.00820708     -0.00501483     0.95599216  1.00502743  True
 8   -0.01875933    -0.01322191     -0.00553742     0.35053292  1.00555278  True
 9   3.28646223     -0.01875933     3.30522156  0.57247874  0.03669108  False
 9   -0.00338494    -0.01875933     0.01537439  0.44678120  0.98474319  True
10   0.00644749     -0.00338494     0.00983243  0.04497439  0.99021575  True
... so on

Answer 1

My original answer is below, but I don't think it solves your problem, because your problem is a bit deeper...

if a move isn't accepted, you need to keep your old configuration, and make a new move from it. If it is accepted, then you make your next move from your new configuration. I do not see any mechanism for doing this properly in your code.

The code appears to always update from the proposed configuration. You need to be able to revert back to the old configuration if the move is not accepted.

Old Answer - not wrong, but not quite right...

Your current code updates the configuration regardless of the Metropolis criteria passing or not - proposed_configuration = update_configuration(proposed_configuration) needs to be indented into that if statement above it.

You are also only appending energy if the move succeeds, but, you should be appending the old energy if your move fails. This will mess up your average energy if you only save the energy of the configuration when a move is accepted. A failed move does not negate the current configuration, a failed move means the move was to the current configuration rather than the proposed configuration.

Allen and Tildesley have lots of good Python code examples for Monte Carlo: https://github.com/Allen-Tildesley/examples/blob/master/python_examples/mc_nvt_lj.py

Answer 2

A properly-run Metropolis Monte Carlo algorithm does not perform global optimisation. It returns the stationary distribution of a random process.

In this case, you have specified an acceptance rule that higher energies should be accepted at an exponentially-decreasing rate. And guess what? Your plot indeed shows that higher energies occur, and exponentially less often as the energy increases. In other words you have sampled the Boltzmann distribution of whatever process it is you were modelling.

Global optimisation is a completely different problem. The only Monte Carlo approach I know of to global optimisation is simulated annealing, which is not what you have implemented.