====== Monte Carlo simulations for the estimation of molecule pair interatcion ====== 
===== 
In this exercise you will perform a MC simulation for different coverages of "sumanene" mlecules
adsorbed on a Ag(111) substrate hypothesizing different posisble values for the nearest neighbor
energy in the molecule-molecule interaction =====


<note tip>
All files of this exercise can be downloaded from the wiki: {{exercise_3.zip|}} 
Download the exercise into your $HOME folder and unzip it. 

<code>
you@eulerX ~$ wget http://www.cp2k.org/_media/exercises:2017_ethz_mmm:exercise_3.zip
you@eulerX ~$ unzip exercises:2017_ethz_mmm:exercise_3.zip
you@eulerX ~$ cd exercise_3
</code>

</note>

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]

{{:exercises:2017_ethz_mmm:coverage.jpg?200|}}

<note tip>
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

{{:exercises:2017_ethz_mmm:hex_0.jpg?200|}}

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

{{:exercises:2017_ethz_mmm:plots_0.0_200.0.png?600|}}

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
</note>

<note warning>
**TASK 1**
Execute the program with parameters:

<code>

coverage 0.02
number of cycles 200
binding energy in eV 0.0
Temperature in K 200
</code>
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
To compute the concentration consider that at coverage 0.02, in the simulation, the total number of molecules is 50
</note>

<note warning>
**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 and 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 main MC steps in the MAIN part of the code

<code>
clusters_plot=[]
for i in range(nouter):

and also the section

#### DECIDE whether to accept or not the move
</code>
What do you think it would happen if you replace the condition
<code>
if np.random.random()<np.exp(-beta*deltae)
</code>
with the condition
<code>
if enew<eold
</code>
</note>
<note warning>
The most complex (and unefficient) sections of the code are the two functions
<code>
def allconnected(m,id,nx,ny)
</code>
which finds out all teh molecules that are connected to a given one
<code>
def neighbors(a,id,nx,ny)
</code>
which finds out all the molecules that are 1st neighbors to a given one.

The function "allconnected" is quite intuitive and inefficient. Can
you imagine roughly a more efficient function to perform the same task? 
</note>
<note warning>
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?**
</note>