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 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 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 embedded atom model (EAM), a computationally simple, yet effective model for metals.

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.

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$.

Here, we will use the so-called phase coexistence technique developed by 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.

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

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 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:

###################################
## 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.

Use geo.in to perform a geometry optimization of cyclohexaphenylene adsorbed on the Cu(111) surface.

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