====== QM/MM: Cyclohexaphenylene on Cu(111) ======

Systems containing thousand atoms and more (e.g. polymers, reactions in solution, molecules on surfaces) are hard to describe fully //ab initio//. In many applications, however, the main attention focusses on a small subset //A// of the large system (active sites, the dissolved/adsorbed molecules). While the interaction of part //A// with the remaining part //B// cannot be neglected, one can often resort to treating //B// by less computationally expensive molecular mechanics (MM) models that interact with a quantum mechanical (QM) description of //A//. The usefulness of this hybrid QM/MM approach was recognized by the [[http://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/press.pdf|2013 Nobel Prize in Chemistry]], awarded to A. Warshel, M. Levitt and A. Karplus for "the development of multiscale models for complex chemical systems".

In this exercise, we are studying the adsorption of cyclohexaphenylene on Cu(111), following the work by [[doi>10.1140/epjb/e2010-00038-1|Pignedoli et al.]].
We will start by examining the molecular mechanics model for the copper substrate, calculating the surface energy and the melting point. We then go on to add the cyclohexaphenylene molecule, which is described quantum mechanically, and study its interaction with the Cu(111) surface.

The Cu(111) substrate will be described by the [[http://en.wikipedia.org/wiki/Embedded_atom_model|embedded atom model]] (EAM), a computationally simple, yet effective model for metals.

<note>**TASK 1**

Read through the first two pages of the paper by [[doi>10.1103/PhysRevB.33.7983|Foiles et al.]] describing the embedded atom model and give a summary of the EAM in your own words. (2P)
</note>

Cutting a solid into pieces breaks bonds and thus costs energy. The //surface energy// is the energy required to create a unit of surface area.

One way to simulate a surface within periodic boundary conditions is the so-called //slab geometry//.
A slab has an upper and a lower surface and the total energy $E(N)$ of the slab contains the energetic contributions
from both surfaces. If the slab is too thin, upper and lower surface may interact, causing their energetic contributions to deviate from the surface energy of a semi-infinite slab. One therefore increases the number of layers $N$ until the energy $\triangle E(N) = E(N) − E(N − 1)$ for adding one layer to the slab converges. At this point, the energy for adding one layer has reached its bulk value $\triangle E(\infty)$ and the surface energy $\sigma$ can be extracted via the formula
$$\sigma =\lim_{N\rightarrow\infty} \frac{1}{2} \frac{1}{A} (E(N)−N \triangle E(N))$$
where $A$ is the surface area.

We have prepared a set of Cu(111) slabs with thicknesses ranging from 1 to 9 layers. Use the provided script ''./run-slabs.sh'' to optimize their geometry with CP2K.
<note tip>
Have a quick look inside ''run-slabs.sh'' to see what it does.
You may need to adapt the name of the CP2K executable.
</note>

<note>**TASK 2**

  - Visualize an optimized slab with vmd, draw the simulation box and create a replica along x,y and z to familiarize yourself with the geometry. Why is the //slab geometry// used in atomistic simulations? What are the relevant convergence parameters? (2P)
  - The initial geometries were cut from a perfect bulk Cu crystal. Why do they still need to be optimized?
  - **BONUS** Do //all// of them need to be optimized?
  - How many layers are required to achieve a bulk environment in the center of the slab? //Hint:// Study the energy $\triangle E(N)$ as a function of $N$. You may want to use the script ''analyze-slabs.sh''.
  - Calculate the surface energy $\sigma$ of Cu(111).
  - Compare with table V of the work by [[doi>10.1103/PhysRevB.33.7983|Foiles et al.]] //Note:// $1\,\text{erg} = 10^{-7}\,\text{J}$.
</note>

Another important parameter of a potential for a metal is its melting temperature $T_m$.
The naive way of determining $T_m$ would be to set up a bulk crystal, perform MD simulations at different temperatures and note at which temperature the system melts and the crystal structure is destroyed.
The problem with this approach is that the lack of any surfaces makes it very hard for the liquid phase to nucleate in a small cell and within a reasonable simulation time. In such a simulation, melting will typically be observed only at temperatures significantly above $T_m$.

{{ folded.png?700 |Half-molten Cu slab.}}

Here, we will use the so-called //phase coexistence technique// developed by [[doi>10.1007/978-94-011-1729-6|Ercolessi et al.]].
We work with a bulk cell with length $l_z$ significantly larger than $l_x$ and $l_y$ and periodic boundary conditions (shown above with z-axis along horizontal direction). In a first step, we have molten the upper half of the crystal, keeping the atoms in the lower half artificially fixed.

Starting from the structure of the half-molten cell, perform an MD simulation at constant energy with an initial temperature $T_0$ somewhere near where you expect the melting temperature.

<note>**TASK 3**

  - Before your start, look inside ''equilibrate.in'' and define the initial temperature. What type of [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/MD.html#desc_ENSEMBLE|ensemble]] are we using? Why aren't we using the NVE ensemble? (2P)
  - The simulation time in ''equilibrate.in'' is 1 ns, which may be longer than required. While the simulation is running, monitor the temperature and the atomic structure (see hint below) to find out when the system reaches equilibrium.
  - Explain, why the average value of the temperature $\langle T \rangle$ must converge to $T_m$ in equilibrium. Calculate $\langle T \rangle$ for the equilibrated system and compare against the [[https://sites.google.com/site/eampotentials/Cu|literature value]] calculated for this potential. //Note:// An error of 2-3% is still expected due to finite size effects. (2P)
  - What happens if you choose $T_0$ much lower (or larger) than $T_m$? Describe qualitatively the expected effects in the atomic structure as well in the temperature as a function of simulation time.
</note>
<note tip>
Since the cell may change during the simulation, the trajectory is saved in the binary [[http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/dcdplugin.html|DCD format]].
It contains the cell, but it does not contain information on the atomic types.
When visualizing it with VMD, one therefore also needs to provide an ''.xyz'' file containing a single frame (and tell VMD the cell for this frame):
<code bash>
vmd -f half-molten.xyz run-CU-EQUIL-pos-1.dcd
pbc set {10.2247 17.7098 100.1818} -first 0 -last 0
pbc wrap -all
</code>
</note>

Chemically inert surfaces can influence the reaction pathways of adsorbed molecules by van der Waals interactions, making the molecule more planar. This is one of the motivations for studying the adsorption of cyclohexaphenylene (CHP), a molecular precursor for [[http://en.wikipedia.org/wiki/Graphene|graphene]], on Cu(111).

We have prepared an input file ''geo.in'' for performing a geometry optimization of cyclohexaphenylene on Cu(111) in the QM/MM approach.
Since this file is a bit more complicated than usual, let's have a look at its overall structure:
<code bash>
###################################
## CHP (PM6) on Cu(111) (EAM)   ###
###################################

@SET LIBDIR ../cp2klib

&GLOBAL
  PROJECT run-GEO
  RUN_TYPE GEO_OPT          # we want to perform a geometry optimization
  PRINT_LEVEL LOW
&END GLOBAL

&MOTION                     # Describes what geometry optimization algorithm to use
  ...                       # and which quantities to plot while doing so.
    OPTIMIZER LBFGS         # The algorithm is completely independent of *how*
  ...                       # the forces are actually computed.
&END MOTION                 

&MULTIPLE_FORCE_EVALS       # notify CP2K that this input file contains multiple
  MULTIPLE_SUBSYS .TRUE.    # FORCE_EVAL sections
  FORCE_EVAL_ORDER 2 3
&END

########################
## How to mix QM & MM ##
########################

&FORCE_EVAL                 # The first FORCE_EVAL section only has an organizational role.
  ...                       # It defines that the total energy of the combined QM/MM system
    MIXING_FUNCTION E1+E2   # should be obtained as a sum of the QM energy (E1, from FORCE_EVAL 2)
  ...                       # and the MM energy (E2, from FORCE_EVAL 3).
END FORCE_EVAL              

#######################
## QM - CHP with PM6 ##
#######################

&FORCE_EVAL                 # The second FORCE_EVAL is responsible for the quantum-mechanical
  ...                       # description of the CHP molecule
END FORCE_EVAL  

###############################################
## MM - Cu(111) (EAM) + interaction with CHP ##
###############################################

&FORCE_EVAL                 # The third FORCE_EVAL section describes the MM model
  ...                       # for the Cu(111) slab *and* the MM model for the interaction of the
END FORCE_EVAL              # Cu(111) slab with the adsorbed molecule.
</code>
Use ''geo.in'' to perform a geometry optimization of cyclohexaphenylene adsorbed on the Cu(111) surface.


<note>**TASK 4**

  - Read through section 2.1 of the work by [[doi>10.1140/epjb/e2010-00038-1|Pignedoli et al.]] and report the methods used for the different kinds of interactions. //Note:// In ''geo.in'', we are using a less computationally demanding method to treat the QM part. (2P)
  - Look inside ''geo.in'' and describe its overall structure. How is the coupling of QM and MM represented in the input file?
  - Run the geometry optimization, visualize the optimized geometry with VMD and render a nice image for the report.
  - What is the adsorption height of the molecule?
  - In order to calculate the //binding energy// of CHP adsorbed on Cu(111), perform geometry optimizations of the isolated molecule as well as the clean slab. //Hint:// Start from ''geo.in'' and remove the sections that you don't need. Make sure that the ''.xyz'' files you use in the input only contain the atoms you need for the respective calculation.
  - Calculate the //interaction energy// of CHP adsorbed on Cu(111). //Hint:// How is the interaction energy defined? ''RUNTYPE ENERGY'' will calculate the total energy without optimizing the geometry.
</note>
<note tip>
In order to calculate the adsorption height, use the [[http://www.ks.uiuc.edu/Research/vmd/current/ug/node116.html|Tcl text interface]] of VMD to calculate the geometric mean of the molecule and surface layer of the slab. Type either directly on the terminal or open the slightly more convenient Extensions -> Tk Console.
<code>
set mol [atomselect top "element C or element H"]  # select molecule
measure center $mol                                # calculate geometric mean of coordinates in selection
set lay [atomselect top "element Cu and z > 12"]   # select surface layer of slab
</code>
</note>