Differences

This shows you the differences between two versions of the page.

--- exercises:2015_pitt:gga [2015/02/26 12:56] – vondele
+++ exercises:2015_pitt:gga [2015/03/03 08:21] – [1. Task: Familiarize yourself with system and setup] vondele
@@ Line 23: / Line 23: @@
   * Use vmd to vizualize the geometries (provided below) named ''mode1.xyz'' and ''mode2.xyz''
   * To edit the input files provided below, use an editor such as ''vi'' or ''nano''. While ''nano'' is simple to use, ''vi'' can be [[http://www.cp2k.org/tools| configured to colour-code cp2k inputs]].
-  * You will need files named ''BASIS_MOLOPT'' ''GTH_POTENTIALS'' ''dftd3.dat'' that are provided as part of CP2K in a directory ''cp2k/data''
+  * You will need files named ''BASIS_MOLOPT'' ''GTH_POTENTIALS'' ''dftd3.dat'' that are provided as part of CP2K in a directory ''cp2k/data'', unless the code has been compiled with the proper flag ('' -D__DATA_DIR'') so that these are found automatically.
-  * Use a job script to submit jobs on the cluster, an example job submission script might look like (TODO: adjust for environment)
+  * Use a job script to submit jobs on the cluster, an example job submission script might look like
 <code - job>
-#!/bin/bash --login
+#PBS -N mode1
-#SBATCH --job-name=mode1
+#PBS -j oe
-#SBATCH --nodes=6
+#PBS -q dist_small
-#SBATCH --time=0:10:00
+#PBS -l nodes=4:ppn=16
-#SBATCH --account=ch5
+#PBS -l walltime=10:00
-aprun -n 144 -N 24 -d 1 ../../cp2k.popt -i mode1.inp -o mode1.out
+#PBS -A cp2k2015
-</code>
+cd $PBS_O_WORKDIR
-===== 2. Task: Bond induced density differences =====
+module purge
-Compute the density difference induced by the bonding for the first binding mode.
+module load cp2k/2.6
-For this you will have to run three separate energy calculations:
-  - combined system bound in the first mode (file ''mode1.xyz'')
-  - lone acetic acid molecule (just remove slab's coordinates from ''mode1.xyz'')
-  - lone TiO$_2$ slab (just remove the acid's coordinates from ''mode1.xyz'')
-In order to output the electronic densities as cube files, your input file has to contain the following snipped:
+prun cp2k.popt -i mode1.inp -o mode1.out
-<code>
-&DFT
-  &PRINT
-    &E_DENSITY_CUBE
-    &END E_DENSITY_CUBE
-  &END
-&END DFT
 </code>
-<note tip>
-The calculations involving the large TiO$_2$ slab should be run on 16 nodes with ''bsub  -n 16''.
-</note>
-To process the cube files we are going to use the cubecruncher tool. It is part of CP2K, but not installed on brutus.
+===== 2. Task: Binding induced density differences =====
-Therefore, a compiled binary of the tool is provided at ''/cluster/home03/matl/schuetto/cubecruncher.x''. Before invoking the cubecruncher, you have to load the cp2k module. To simplify the following steps you should create a symbolic link to it in your working directory:
-<code>
-you@brutusX ~$ ln -s /cluster/home03/matl/schuetto/cubecruncher.x .
-</code>
+We start with single point energy calculations on binding mode 1, to visualize the interaction between molecule and surface. The goal is to compute the binding induced density difference:
+\[ \rho_\text{induced}= \rho_\text{slab-dye-complex} - \rho_\text{dye} - \rho_\text{slab} \]
+First, we'll discuss in detail the structure and the choices made in the sample input file ''mode1.inp''.
+topics:
+  * Project name
+  * Runtype
+  * Gaussian Basis, pseudopotentials
+  * PW Cutoff
+  * thresholds
+  * SCF: OT
+  * XC and -D3 correction
+  * Unit cell choice
+Second, we run the cp2k input and store the output for analysis and discussion.
-Now subtract the densities of the lone systems from the bonded system:
 <code>
-you@brutusX ~$ ./cubecruncher.x -i mode1-ELECTRON_DENSITY-1_0.cube -subtract mode1_acid-ELECTRON_DENSITY-1_0.cube -o tmp.cube
+cp2k.popt -i mode1.inp -o mode1.out
-you@brutusX ~$ ./cubecruncher.x -i tmp.cube -subtract mode1_slab-ELECTRON_DENSITY-1_0.cube -o mode1_delta.cube
 </code>
+In addition to the output ''mode1.out'', other files are gerenated by CP2K named MODE1*
-The generated cube file is not aligned with the simulation cell. Center the cube file with the cubecruncher.x tool:
+topics:
+  * General overiew
+  * OT output
+  * Various grid quantities
+  * Density cube output
+  * Timing report
+Third, we compute the changes in density induced by the binding. For this you will have to run three separate energy calculations:
+  - combined system bound in the first mode (file ''mode1.xyz'')
+  - lone acetic acid molecule (just remove slab's coordinates from ''mode1.xyz''), name the file ''mode1_dye.xyz''
+  - lone TiO$_2$ slab (just remove the acid's coordinates from ''mode1.xyz''), name the file ''mode1_slab.xyz''
+Create the .xyz files (check with vmd that they contain the right subsystems), and create mode1_dye.inp, mode1_slab.inp by changing both ''COORD_FILE_NAME'' and ''PROJECT'' accordingly.
+After computing these input files, we analyze the results using a tool provided with cp2k ''cubecruncher.x'' that manipulates cube files (e.g. can compute difference).
 <code>
-you@brutusX ~$ cubecruncher.x -center geo -i mode1_delta.cube -o mode1_delta-centered.cube
+~$ cubecruncher.x -i MODE1-ELECTRON_DENSITY-1_0.cube -subtract MODE1_dye-ELECTRON_DENSITY-1_0.cube -o tmp.cube
+~$ cubecruncher.x -i tmp.cube -subtract MODE1_slab-ELECTRON_DENSITY-1_0.cube -center geo -o MODE1_delta.cube
 </code>
-You can visualize the resulting file ''mode1_delta-centered.cube'' with VMD. This has been covered in a [[mo_ethene| previous exercise]].
+You can visualize the resulting file ''MODE1_delta.cube'' with VMD.
 What you get should look similar to this:
-{{ dye_tio_bonding_density.png?300 |}}
+{{ exercises:2015_ethz_mmm:dye_tio_bonding_density.png?300 |}}
+===== 3. Task: relative stabilities  =====
+In order to compute the relative stability of mode1 and mode2, both configurations need to be geometry optimized.
+To do so, turn off the generation of cubes (''&E_DENSITY_CUBE OFF'') in mode1.inp, change to ''RUN_TYPE GEO_OPT'' and adjust the project name. Create and run a similar input file for mode2.
+input topics:
+  * BFGS vs LBFGS
+  * EPS_SCF, CUTOFF, MAX_DR, ..
+output topics:
+  * ''Informations at step''
+  * Trajectory ''MODE1_GEO-pos-1.xyz''
+Compare the final energies (''ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):''), and determine which mode is most stable. Does this agree with the values in table 1 of the manuscript cited ?
+===== 4. Task: ab initio molecular dynamics  =====
-===== 3. Task: Bonding energies  =====
+<note>This task is optional, and can be done near the end if time is available. Ab initio MD will be briefly discussed in a next exercise.</note>
-Compute the binding energy for both binding modes:
-\[ E_\text{binding}=\sum E_\text{products} - \sum E_\text{reactants} \]
+Perform a short ab initio molecular dynamics simulation of the system (~1000 steps, ~0.5ps) by changing to ''RUN_TYPE MD''. After a couple of hours the job should be finished. Now analyze the OH distance in VMD. A possible outcome is :
-For this you will need the energy values of four systems:
+{{exercises:2015_pitt:acetic_acid_mode1_md.png?600 |}}
-  - lone acetic acid molecule (run geometry optimization, use energy of last step)
-  - lone TiO$_2$ slab  (you can use the already geometry optimized coordinates from ''relaxed_slab.xyz'' at the end of the exercise)
-  - combined system bound in the first mode (can be reused from previous task)
-  - combined system bound in the second mode (file ''mode2.xyz'')
-<note important>
+What can you say about the hydrogen bond to the surface, relative acidity of the two oxygens ?
-You can not reuse the energy values for the lone sub-systems from the previous task. Since the unbound subsystems might relax into a different geometry, they have to be geometry optimized first. This has been covered in a
+Note that, in order to be statistically relevant, longer trajectories should be employed, and surface slab thickness will play an important role. Also compare to Fig. 7 of the paper referenced.
-[[geometry_optimization|previous exercise]].
+====== Required Files ======
-</note>
+(right) click on the filename to download to your local machine.
-===== Required Files =====
 <code - mode1.inp>
 &GLOBAL
@@ Line 121: / Line 145: @@
     &MGRID
        ! PW cutoff ... depends on the element (basis) too small cutoffs lead to the eggbox effect.
-       ! certain calculations (e.g. geometry optimization, vibrational frequencies, NPT and cell optimizations, need higher cutoffs)
+       ! certain calculations (e.g. geometry optimization, vibrational frequencies,
+       ! NPT and cell optimizations, need higher cutoffs)
        CUTOFF [Ry] 400
     &END
     &QS
-       METHOD GPW  ! use the GPW method (i.e. pseudopotential based calculations with the Gaussian and Plane Waves scheme).
+       ! use the GPW method (i.e. pseudopotential based calculations with the Gaussian and Plane Waves scheme).
-       EPS_DEFAULT 1.0E-10 ! default threshold for numerics ~ roughly numerical accuracy of the total energy per electron, sets reasonable values for all other thresholds
+       METHOD GPW
-       EXTRAPOLATION ASPC ! used for MD, the method used to generate the initial guess.
+       ! default threshold for numerics ~ roughly numerical accuracy of the total energy per electron,
+       ! sets reasonable values for all other thresholds.
+       EPS_DEFAULT 1.0E-10
+       ! used for MD, the method used to generate the initial guess.
+       EXTRAPOLATION ASPC
     &END
@@ Line 185: / Line 214: @@
     &END CELL
-    ! atom coordinates can be in the &COORD section, or provided as an external file.
+    ! atom coordinates can be in the &COORD section,
+    ! or provided as an external file.
     &TOPOLOGY
       COORD_FILE_NAME mode1.xyz
@@ Line 191: / Line 221: @@
     &END
-    ! MOLOPT basis sets are fairly costly, but in the 'DZVP-MOLOPT-SR-GTH' available for all elements
+    ! MOLOPT basis sets are fairly costly,
-    ! their contracted nature makes them suitable for condensed and gas phase systems alike.
+    ! but in the 'DZVP-MOLOPT-SR-GTH' available for all elements
+    ! their contracted nature makes them suitable
+    ! for condensed and gas phase systems alike.
     &KIND H
       BASIS_SET DZVP-MOLOPT-SR-GTH
@@ Line 225: / Line 257: @@
    TEMPERATURE [K] 300
    TIMESTEP [fs] 0.5
-   STEPS 100
+   STEPS 1000
    # GLE thermostat as generated at http://epfl-cosmo.github.io/gle4md
    # GLE provides an effective NVT sampling.