# Python script for Metropolis Monte Carlo (global optimisation)

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

• Is this for a monatomic LJ system? Was the initial configuration generated at random, or on a lattice? What is your y-axis in the plot? You are using reduced units? There is a lot of information you are leaving out... Commented Feb 22, 2023 at 17:32

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

• I edited the question after I gained some knowledge from your reply. Could you please give me more advice on the edited question? I appreciate any advice. Commented Feb 23, 2023 at 16:45
• @DGKang sorry we prefer for users to ask new questions in new posts rather than adding more questions to a post ad infinitum. Commented Feb 23, 2023 at 17:09
• In short: No, your plot is not Markov-Chain. You are definitely not maintaining Detailed Balance. If you want a global optimization as your main goal, that is an unsolved problem. We usually use MD with a force minimization. Some sort of simulated annealing would be your best bet. Markov-chain MC is useful for calculating average properties, not max/mins. That said, the configurations at equilibrium will be not far from what you are after. Commented Feb 23, 2023 at 19:40

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.

• this is a good point. I purposefully ignored the global optimisation comment, assuming it was poor wording, but, if that was the goal, the script makes a bit more sense... still isn't right though. A Markov chain converges to equilibrium, not a min or max. Commented Feb 23, 2023 at 16:43
• I edited the question after I gained some knowledge from your reply. Could you please give me more advice on the edited question? I appreciate any advice. Commented Feb 23, 2023 at 16:45
• @DGKang sorry we prefer for users to ask new questions in new posts rather than adding more questions to a post ad infinitum. Commented Feb 23, 2023 at 17:09