This is an old revision of the document!
Table of Contents
Validation of a KCl QMMM model
(exercise by Matthew Watkins, University college, London)
In this exercise you will validate the mixed quamtum/classical model for a KCl slab. Hte : 10.1021%2Fja505936b.
- In the first part of the exercise you will consider the optimized configuration (already in the directory) and study the pure electronic adsorption energy, namely the difference between the total energy of the surface-molecule system and the energy of the molecule alone and surface alone in the same geometry as the surface-molecule system minimum structure. This will allow to show the binding pattern of the electronic density.
- In the second part, you will optimize the surface and the molecule separately; this will allow to compute the total adsorption energy.
1. Task: Familiarize yourself
The coordinates of the optimized configuration are provided to you as S_M.opt.xyz
(S stands for “Substrate”, M for “Molecule”, opt for “optimized”). Visualize the geometry with VMD and familiarize yourself with the system.
2. Task: Bond induced density differences
Compute the density difference induced by the adsorption bonding. For this you will have to run three separate energy calculations, using the *.ene.inp files.
- combined system (file
S_M.opt.xyz
) - lone acetylene (file
M.S_M.xyz
) - lone slab (file
S.S_M.xyz
)
In order to output the electronic densities as cube files, your input file has to contain the following snipped:
&DFT &PRINT &E_DENSITY_CUBE &END E_DENSITY_CUBE &END &END DFT
qsub run -v INP=prefix
. Check the run
file for the number of nodes.
To process the cube files we are going to use the cubecruncher tool. It is part of CP2K and is in your exercise directory.
you@eulerX ~$ ./cubecruncher.x -i S_M-ELECTRON_DENSITY-1_0.cube -subtract S-ELECTRON_DENSITY-1_0.cube -o tmp.cube you@eulerX ~$ ./cubecruncher.x -i tmp.cube -subtract M-ELECTRON_DENSITY-1_0.cube -o Delta_ads.cube
The generated cube file is not aligned with the simulation cell. Center the cube file with the cubecruncher.x tool:
you@eulerX ~$ ./cubecruncher.x -center geo -i Delta_ads.cube -o Delta_ads-centered.cube
You can visualize the resulting file delta_ads-centered.cube
with VMD. This has been covered in a previous exercise.
3. Task: Bonding energies
Compute the binding energy:
\[ E_\text{binding}=\sum E_\text{products} - \sum E_\text{reactants} \]
For this you will need the energy values of three systems:
- lone acetylene molecule (run geometry optimization, use energy of last step)
- lone slab (you can use the already geometry optimized coordinates from
S.opt.xyz
at the end of the exercise) - combined system adsorbed (can be reused from previous task)
Questions
- Sketch briefly the geometry of the molecule when adsorbed and in the gas phase.
- Report the system energy for the bonded system, lone slab, and lone molecule.
- Can you estimate the contribution due to the geometry relaxation?
- Briefly report the bond induced density difference on the system.
Required Files
/home/psd/Exercise_9
. Change the name of the xyz file accordingly in the input files.
- input.inp
@SET METHOD = QMMM # FIST all classical treatment # QS all quantum treatment &GLOBAL FLUSH_SHOULD_FLUSH PRINT_LEVEL LOW PROJECT KCl RUN_TYPE GEO_OPT &END GLOBAL &FORCE_EVAL METHOD $METHOD @include QS.inc @include MM.inc &QMMM #this defines the QS cell in the QMMM calc &CELL ABC 12.6 15.0 12.6 PERIODIC XZ &END CELL ECOUPL GAUSS # use GEEP method NOCOMPATIBILITY USE_GEEP_LIB 6 # use GEEP method &PERIODIC # apply periodic potential #in this case QM box = MM box in XZ so turn #off coupling/recoupling of the QM multipole &MULTIPOLE OFF &END &END PERIODIC #these are just the ionic radii of K Cl #but should be treated as parameters in general #fit to some physical property &MM_KIND K RADIUS 1.52 &END MM_KIND &MM_KIND Cl RADIUS 1.67 &END MM_KIND #define the model &QM_KIND K MM_INDEX 25..32 41..48 &END QM_KIND &MM_KIND Cl RADIUS 1.67 &END MM_KIND #define the model &QM_KIND K MM_INDEX 25..32 41..48 &END QM_KIND &QM_KIND Cl MM_INDEX 17..24 33..40 &END QM_KIND &END QMMM &SUBSYS #this defines the cell of the whole system #must be orthorhombic, I think &CELL ABC 12.6 100.0 12.6 &END CELL &TOPOLOGY COORD_FILE_NAME kcl.xyz COORD_FILE_FORMAT XYZ &GENERATE &ISOLATED_ATOMS #ignores bonds dihedrals etc in classical part LIST 1..48 &END &END &END &KIND K ELEMENT K BASIS_SET DZVP-MOLOPT-SR-GTH POTENTIAL GTH-PBE-q9 &END KIND &KIND Cl BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q7 &END &END SUBSYS &END FORCE_EVAL #should be able to use most motion sections #analytic stress tensor not available, I think @include motion.inc
and includes as separate files, using the @include macro, for the QS, MM and motion sections
- QS.inc
&DFT BASIS_SET_FILE_NAME BASIS_MOLOPT POTENTIAL_FILE_NAME GTH_POTENTIALS &MGRID COMMENSURATE CUTOFF 150 &END MGRID &QS EPS_DEFAULT 1.0E-12 &END QS &SCF EPS_SCF 1.0E-06 MAX_SCF 26 SCF_GUESS RESTART &OT MINIMIZER CG PRECONDITIONER FULL_SINGLE_INVERSE ENERGY_GAP 0.001 &END OT &OUTER_SCF EPS_SCF 1.0E-05 &END OUTER_SCF &END SCF &XC &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL &END XC &PRINT &MO_CUBES NLUMO 10 WRITE_CUBE T &END MO_CUBES &V_HARTREE_CUBE STRIDE 2 2 2 &END &END PRINT &END DFT
- MM.inc
&MM &FORCEFIELD &CHARGE ATOM K CHARGE 1.0 &END CHARGE &CHARGE ATOM Cl CHARGE -1.0 &END CHARGE &NONBONDED &WILLIAMS atoms K Cl A [eV] 4117.9 B [angstrom^-1] 3.2808 C [eV*angstrom^6] 0.0 RCUT [angstrom] 3.0 &END WILLIAMS &WILLIAMS atoms Cl Cl A [eV] 1227.2 B [angstrom^-1] 3.1114 C [eV*angstrom^6] 124.0 RCUT [angstrom] 3.0 &END WILLIAMS &WILLIAMS atoms K K A [eV] 3796.9 B [angstrom^-1] 3.84172 C [eV*angstrom^6] 124.0 RCUT [angstrom] 3.0 &END WILLIAMS &END NONBONDED &END FORCEFIELD &POISSON &EWALD EWALD_TYPE spme ALPHA .44 GMAX 40 &END EWALD &END POISSON &END MM
- motion.inc
&MOTION &GEO_OPT OPTIMIZER LBFGS &END &CONSTRAINT &FIXED_ATOMS LIST 1..16 EXCLUDE_MM .FALSE. EXCLUDE_QM .TRUE. &END FIXED_ATOMS &END CONSTRAINT &END MOTION
- kcl.xyz
48 Cl 0.00000 15.00000 0.00000 Cl 3.15000 15.00000 3.15000 Cl 0.00000 15.00000 6.30000 Cl 3.15000 15.00000 9.45000 Cl 6.30000 15.00000 0.00000 Cl 9.45000 15.00000 3.15000 Cl 6.30000 15.00000 6.30000 Cl 9.45000 15.00000 9.45000 K 3.15000 15.00000 0.00000 K 0.00000 15.00000 3.15000 K 3.15000 15.00000 6.30000 K 0.00000 15.00000 9.45000 K 9.45000 15.00000 0.00000 K 6.30000 15.00000 3.15000 K 9.45000 15.00000 6.30000 K 6.30000 15.00000 9.45000 Cl 3.15000 18.15000 0.00000 Cl 0.00000 18.15000 3.15000 Cl 3.15000 18.15000 6.30000 Cl 0.00000 18.15000 9.45000 Cl 9.45000 18.15000 0.00000 Cl 6.30000 18.15000 3.15000 Cl 9.45000 18.15000 6.30000 Cl 6.30000 18.15000 9.45000 K 0.00000 18.15000 0.00000 K 3.15000 18.15000 3.15000 K 0.00000 18.15000 6.30000 K 3.15000 18.15000 9.45000 K 6.30000 18.15000 0.00000 K 9.45000 18.15000 3.15000 K 6.30000 18.15000 6.30000 K 9.45000 18.15000 9.45000 Cl 0.00000 21.30000 0.00000 Cl 3.15000 21.30000 3.15000 Cl 0.00000 21.30000 6.30000 Cl 3.15000 21.30000 9.45000 Cl 6.30000 21.30000 0.00000 Cl 9.45000 21.30000 3.15000 Cl 6.30000 21.30000 6.30000 Cl 9.45000 21.30000 9.45000 K 3.15000 21.30000 0.00000 K 0.00000 21.30000 3.15000 K 3.15000 21.30000 6.30000 K 0.00000 21.30000 9.45000 K 9.45000 21.30000 0.00000 K 6.30000 21.30000 3.15000 K 9.45000 21.30000 6.30000 K 6.30000 21.30000 9.45000