Monte Carlo simulations for the estimation of pair interactions

In this exercise you will perform a MC simulation for different coverages of “sumanene” mlecules adsorbed on a Ag(111) substrate hypothesizing different possible values for the nearest neighbor energy in the molecule-molecule interaction

Please download the program monte_carlo.py from this link

In an experiment performed at Empa, sumanene molecules were adsorbed on a Ag(111) surface. It was found that at very low coverage (0.02) 30% of the molecules were weakly bound into dimers. [ http://dx.doi.org/10.1021/ja504126z J. Am. Chem. Soc. 2014, 136, 13666−13671]

The python code mc.py mimics the Ag(111) substrate as a rigid hexagonal lattice onto which molecules (blue dots) are allowed to move by random discrete steps

While in execution, the code will show you snapshots of the positions of the molecules on the lattice on the left panel.

In the central panel average values for the number of isolated molecules, the number of dimers and the number of clusters is plotted. On the right panel the average energy of the system is plotted. Intermediate average values are also written in a output file data_de_T.out where de and T are the input values of the dimer energy and of the Temperature also a final snapshot of the graphic panels is saved as png image

TASK 1 Execute the program

python monte_carlo.py

with parameters:

coverage 0.02
number of cycles 200
binding energy in eV 0.0
Temperature in K 200

Estimate the dimer binding energy as DE=kT * ln(n0/nexp) where k is Boltzmann's constant T is the simulation (and experiment) temperature n0 is the concentration of dimers in the case of zero interaction nexp is the concentration of dimers found in the experiment

warning To compute the concentration consider that at coverage 0.02, in the simulation, the total number of molecules is 50

TASK2 repeat the simulation using as DE your estimate. What do you get?

TASK3 Repeat the simulation with coverage 0.1 and DE=-0.02 then DE=-0.1 Describe what you obtain. Now try coverage=0.1 T=400 DE=-0.1 Comment the result

TASK4 have a look at the pyhton code, identify the outer loop, mainly responsible of plotting and statistics

#### outer loop of the siulation
for i in range(nouter):

and then identify the inner loop where the Monte Carlo moves happen:

#### inner loop of 1000 steps

and

#### DECIDE whether to accept or not the move

What do you think it would happen if you replace the condition

if np.random.random()<np.exp(-beta*deltae)

with the condition

if enew<eold

TASK 5 The most complex (and unefficient) sections of the code are the two functions

def allconnected(m,id,nx,ny)

which finds out all teh molecules that are connected to a given one

def neighbors(a,id,nx,ny)

which finds out all the molecules that are 1st neighbors to a given one.

To understand these sections you have to refer to the graph of the hexagonal lattice with the convention adopted for the coordinate system. The function “allconnected” is quite intuitive and inefficient. Can you imagine roughly a more efficient function to perform the same task? Just think how you would do the task of identifying, among a set of particles, the ones that are alone, the ones that form an isolated dimer and thr ones that form clusters. Google can help you :) (you do not have to find a solution)

TASK 6 The code, at each step, moves a particle chosen randomly to a new site chosen randomly. Would it be correct to move all particles in a step? what would change?