exercises:2015_uzh_molsim:chp_cu111
Differences
This shows you the differences between two versions of the page.
— | exercises:2015_uzh_molsim:chp_cu111 [2020/08/21 10:15] (current) – created - external edit 127.0.0.1 | ||
---|---|---|---|
Line 1: | Line 1: | ||
+ | ====== 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, | ||
+ | |||
+ | In this exercise, we are studying the adsorption of cyclohexaphenylene on Cu(111), following the work by [[doi> | ||
+ | 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, | ||
+ | |||
+ | The Cu(111) substrate will be described by the [[http:// | ||
+ | |||
+ | < | ||
+ | |||
+ | Read through the first two pages of the paper by [[doi> | ||
+ | </ | ||
+ | |||
+ | 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 '' | ||
+ | <note tip> | ||
+ | Have a quick look inside '' | ||
+ | You may need to adapt the name of the CP2K executable. | ||
+ | </ | ||
+ | |||
+ | < | ||
+ | |||
+ | - 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? | ||
+ | - 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 '' | ||
+ | - Calculate the surface energy $\sigma$ of Cu(111). | ||
+ | - Compare with table V of the work by [[doi> | ||
+ | </ | ||
+ | |||
+ | 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? | ||
+ | |||
+ | Here, we will use the so-called //phase coexistence technique// developed by [[doi> | ||
+ | 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. | ||
+ | |||
+ | < | ||
+ | |||
+ | - Before your start, look inside '' | ||
+ | - The simulation time in '' | ||
+ | - 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:// | ||
+ | - 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 tip> | ||
+ | Since the cell may change during the simulation, the trajectory is saved in the binary [[http:// | ||
+ | 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 '' | ||
+ | <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 | ||
+ | </ | ||
+ | </ | ||
+ | |||
+ | Chemically inert surfaces can influence the reaction pathways of adsorbed molecules by van der Waals interactions, | ||
+ | |||
+ | We have prepared an input file '' | ||
+ | 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 | ||
+ | PRINT_LEVEL LOW | ||
+ | &END GLOBAL | ||
+ | |||
+ | & | ||
+ | ... # 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_SUBSYS .TRUE. | ||
+ | FORCE_EVAL_ORDER 2 3 | ||
+ | &END | ||
+ | |||
+ | ######################## | ||
+ | ## How to mix QM & MM ## | ||
+ | ######################## | ||
+ | |||
+ | & | ||
+ | ... # 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 ## | ||
+ | ####################### | ||
+ | |||
+ | & | ||
+ | ... # description of the CHP molecule | ||
+ | END FORCE_EVAL | ||
+ | |||
+ | ############################################### | ||
+ | ## MM - Cu(111) (EAM) + interaction with CHP ## | ||
+ | ############################################### | ||
+ | |||
+ | & | ||
+ | ... # for the Cu(111) slab *and* the MM model for the interaction of the | ||
+ | END FORCE_EVAL | ||
+ | </ | ||
+ | Use '' | ||
+ | |||
+ | |||
+ | < | ||
+ | |||
+ | - Read through section 2.1 of the work by [[doi> | ||
+ | - Look inside '' | ||
+ | - Run the geometry optimization, | ||
+ | - 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 '' | ||
+ | - Calculate the // | ||
+ | </ | ||
+ | <note tip> | ||
+ | In order to calculate the adsorption height, use the [[http:// | ||
+ | < | ||
+ | set mol [atomselect top " | ||
+ | measure center $mol # calculate geometric mean of coordinates in selection | ||
+ | set lay [atomselect top " | ||
+ | </ | ||
+ | </ |