This is an old revision of the document!
Table of Contents
Nanostructures and adsorption on metallic surfaces
Problem: compute activation barrier for the last step (6 to 7) of the cyclodehydrogenation reaction CHP@Cu(111) → TBC
Introduction
The system is composed by a molecule (56 atoms) adsorbed on a Cu(111) slab (~2900 atoms). The molecule is treated within DFT (PM6[1] in the exercise) and the substrate with EAM[2] potentials. The substrate and the molecule are coupled through a par potential mimicking the van der Waals interaction and Pauli repulsion.
Exercise 1: NEB
We need MULTIPLE_FORCE_EVALS:
- In a first section we define how to partition the whole system in a DFT+EMPIRICAL part
- in a second section we define the force fields for the EMPIRICAL part
- and in the third section we define the parameters for the DFT part.
- The fourth and final section defines the target of the calculation(GEO_OPT, NEB, Molecular dynamics..)
First Section
Here we specify that our system consists of two fragments: atoms 1 to 56 and atoms 57 to 2936,
the whole system is contained in the cell specified in FORCE_EVAL/SUBSYS section.
We declare that the first FORCE_EVAL section following the present one will be dedicated to FRAGMENTS 1 and 2 (it will contain the empirical potentials)
and the second FORCE_EVAL section will be devoted to FRAGMENT 1 only (the DFT treatment of the molecule).
&MULTIPLE_FORCE_EVALS
FORCE_EVAL_ORDER 2 3
MULTIPLE_SUBSYS T
&END
&FORCE_EVAL
METHOD MIXED
&MIXED
MIXING_TYPE GENMIX
GROUP_PARTITION 1 1
&GENERIC
ERROR_LIMIT 1.0E-10
MIXING_FUNCTION E1+E2
VARIABLES E1 E2
&END
&MAPPING
&FORCE_EVAL_MIXED
&FRAGMENT 1
1 56
&END
&FRAGMENT 2
57 2936
&END
&END
&FORCE_EVAL 1
DEFINE_FRAGMENTS 1 2
&END
&FORCE_EVAL 2
DEFINE_FRAGMENTS 1
&END
&END
&END
&SUBSYS
&COLVAR
&DISTANCE
ATOMS 49 50
&END
&END
&CELL
ABC 40.764229 39.715716 70.
&END CELL
&TOPOLOGY
COORD_FILE_NAME ./s.xyz
COORDINATE XYZ
CONNECTIVITY OFF
&END TOPOLOGY
&END SUBSYS
&END
Second Section
In this section we define the empirical potentials for the whole system: EAM will be used for the Cu-Cu interactions (parametrized in the file CP2KTUTORIAL/Files/CU.pot), a “C6/R^6 + pauli repulsion” will be defined for C-Cu and H-Cu and, finally, NULL interactions will be defined for C-C, C-H and H-H since these will be obtained from DFT in the next section. The charge of every specie (zero in our example) as well as the interaction within each possible type-pair has to be defined.
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
&SPLINE
EPS_SPLINE 1.0E-6
EMAX_SPLINE 0.9
&END
&CHARGE
ATOM Cu
CHARGE 0.0
&END CHARGE
&CHARGE
ATOM H
CHARGE 0.0
&END CHARGE
&CHARGE
ATOM C
CHARGE 0.0
&END CHARGE
&NONBONDED
&GENPOT
atoms Cu C
FUNCTION A*exp(-av*r)+B*exp(-ac*r)-C/(r^6)
VARIABLES r
PARAMETERS A av B ac C
VALUES 4.13643 1.33747 115.82004 2.206825 75.40708524085266692113
RCUT 15
&END GENPOT
&GENPOT
atoms Cu H
FUNCTION A*exp(-av*r)+B*exp(-ac*r)-C/(r^6)
VARIABLES r
PARAMETERS A av B ac C
VALUES 0.878363 1.33747 24.594164 2.206825 21.32834305207502214234
RCUT 15
&END GENPOT
&LENNARD-JONES
atoms C H
EPSILON 0.0
SIGMA 3.166
RCUT 15
&END LENNARD-JONES
&LENNARD-JONES
atoms H H
EPSILON 0.0
SIGMA 3.166
RCUT 15
&END LENNARD-JONES
&LENNARD-JONES
atoms C C
EPSILON 0.0
SIGMA 3.166
RCUT 15
&END LENNARD-JONES
&EAM
atoms Cu Cu
PARM_FILE_NAME ../../Files/CU.pot
&END EAM
&END NONBONDED
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE none
&END EWALD
&END POISSON
&END
&SUBSYS
&CELL
ABC 40.764229 39.715716 70.
&END CELL
&TOPOLOGY
COORD_FILE_NAME ./s.xyz
COORDINATE XYZ
CONNECTIVITY OFF
&END TOPOLOGY
&END SUBSYS
&END
Third Section
In this section we specify the parameters for the DFT treatment of the molecule. Note that the simulation cell here is reduced with comparison to the simulation cell of the whole system.
&FORCE_EVAL
METHOD Quickstep
&DFT
&QS
METHOD PM6
EXTRAPOLATION USE_GUESS
&SE
&COULOMB
CUTOFF [angstrom] 20.0
&END
&EXCHANGE
CUTOFF [angstrom] 20.0
&END
&END
&END QS
&SCF
MAX_SCF 20
SCF_GUESS RESTART
EPS_SCF 1.0E-6
&OT
PRECONDITIONER FULL_SINGLE_INVERSE
MINIMIZER DIIS
N_DIIS 7
&END
&OUTER_SCF
MAX_SCF 10
EPS_SCF 1.0E-6
&END
&PRINT
&RESTART
&EACH
QS_SCF 0
&END
ADD_LAST NUMERIC
&END
&RESTART_HISTORY OFF
&END
&END
&END SCF
&POISSON
PERIODIC NONE
POISSON_SOLVER WAVELET
&WAVELET
SCF_TYPE 60
&END
&END POISSON
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
! &PRINT
! &MO_CUBES
! NHOMO 15
! NLUMO 15
! WRITE_CUBE T
! LOG_PRINT_KEY
! &END
! &END
&END DFT
&SUBSYS
&CELL
ABC 30 30 30
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME ./f.xyz
COORDINATE xyz
&CENTER_COORDINATES
&END
&END
# &KIND H
# BASIS_SET TZV2P-MOLOPT-GTH
# POTENTIAL GTH-BLYP-q1
# &END KIND
# &KIND C
# BASIS_SET TZV2P-MOLOPT-GTH
# POTENTIAL GTH-BLYP-q4
# &END KIND
&END SUBSYS
&END FORCE_EVAL
Fourth Section
In this section we define the target of the simulation as a geometry optimization this will be MD, or BAND in other examples.
&GLOBAL
PRINT_LEVEL LOW
PROJECT INI
RUN_TYPE GEO_OPT
WALLTIME 28000
&END GLOBAL
&MOTION
&GEO_OPT
MAX_FORCE 0.0001
MAX_ITER 500
OPTIMIZER LBFGS
&LBFGS
&END
&END
&END
!&EXT_RESTART
! RESTART_FILE_NAME constr_3cc_2h_41_dft-1.restart
!&END
Tasks
- First of all obtain the equilibrium geometry for the initial and the final states (in the directories INI and FIN you will find the input files (geo.inp) the output files (out) and the trajectory of the geometry optimization (INI-pos-1.xyz and FIN-pos-1.xyz). Please note that the geometry optimization of the final state is constrained with respect to the position of the final H2 molecule.
- Perform a NEB simulation to find a saddle point between ini_eq.xyz and fin_eq.xyz. A linear guess for the path is not a good idea in this case; a clever guess can be obtained from a series of constrained geometry optimizations where the H-H distance is forced to vary from its initial value (more than 3A)to the H2 equilibrium distance. Such a guess is given, in the directory EX2, with the files ini_eq.xyz, due.xyz,tre.xyz,…,tredici.xyz,fin_eq.xyz).
Check that the initial step of the NEB calculation is reasonable and determine the activation barrier at the and of the optimization. In the directory EX2 you will find the input for the calculation (neb.inp), the output (out) and the optimization trajectory of each image of the system (NEB-pos-Replica_nr_XX-1.xyz 14 images are used in the example)
The last section of the input described above is modified as follows:
&GLOBAL
PRINT_LEVEL LOW
PROJECT NEB
RUN_TYPE BAND
WALLTIME 28000
&END GLOBAL
&MOTION
&BAND
NPROC_REP 2
BAND_TYPE CI-NEB
NUMBER_OF_REPLICA 14
K_SPRING 0.05
&CONVERGENCE_CONTROL
MAX_FORCE 0.0030
RMS_FORCE 0.0050
MAX_DR 0.002
RMS_DR 0.005
&END
ROTATE_FRAMES F
ALIGN_FRAMES F
&CI_NEB
NSTEPS_IT 2
&END
&OPTIMIZE_BAND
OPT_TYPE DIIS
OPTIMIZE_END_POINTS F
&DIIS
MAX_STEPS 1000
&END
&END
#
&REPLICA
COORD_FILE_NAME ini_eq.xyz
&END
&REPLICA
COORD_FILE_NAME due.xyz
&END
.
.
.
.
.
&REPLICA
COORD_FILE_NAME fin_eq.xyz
&END
&PROGRAM_RUN_INFO
INITIAL_CONFIGURATION_INFO
&END
&CONVERGENCE_INFO
&END
&END
&END
We declare that we need 14 replicas and each replica has to be computed with 2 CPUS (if we run the calculation with 14 cpus, then every group of 2 cpus will work on two replicas). The initial guess for the reaction path is given through 14 geometry files.
SCF_GUESS RESTART in combination with EXTRAPOLATION USE_GUESS has to be declared in the FORCE_EVAL/DFT section to allow for an efficient guess of the wavefunctions at each step of the BAND optimization.