CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2015_cecam_tutorial:neb

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
exercises:2015_cecam_tutorial:neb [2015/08/20 09:37] – tmuellerexercises:2015_cecam_tutorial:neb [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 2: Line 2:
  
 Problem: compute activation barrier for the last step (6 to 7) of the cyclodehydrogenation reaction CHP@Cu(111) -> TBC Problem: compute activation barrier for the last step (6 to 7) of the cyclodehydrogenation reaction CHP@Cu(111) -> TBC
 +
 +{{ :exercises:2015_cecam_tutorial:exercise.jpg?direct&500 |}}
 +
 +  * Original author: Carlo Pignedoli
 +  * Complete source and output files: [[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB.tar.xz|NEB.tar.xz]]
  
 ===== Introduction ===== ===== Introduction =====
  
 The system is composed by a molecule (56 atoms) adsorbed on a Cu(111) slab (~2900 atoms). 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 molecule is treated within DFT (PM6[(Stewart2007>J.J.P. Stewart, J. Mol. Model. 13, 1173 (2007))] in the exercise) 
-The substrate and the molecule are coupled through a par potential mimicking the van der Waals interaction and Pauli repulsion.+and the substrate with EAM[(Daw1983>M.S. Daw, M.I. Baskes, Phys. Rev. Lett. 50, 1285 (1983))] potentials. 
 +The substrate and the molecule are coupled through a par potential mimicking the van der Waals interaction and Pauli repulsion.[(Levi2007>A.C. Levi, P.Calvini, Surf. Sci. 601, 1494 (2007))]
  
 ===== Exercise 1: NEB ===== ===== Exercise 1: NEB =====
Line 17: Line 23:
   * and in the third section we define the parameters for the DFT 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..)   * The fourth and final section defines the target of the calculation(GEO_OPT, NEB, Molecular dynamics..)
 +
 +The files for this part of the exercise are found in the **''EX2''** folder.
  
 ==== First Section ==== ==== First Section ====
Line 79: Line 87:
  
 In this section we define the empirical potentials for the whole system: 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), +EAM will be used for the Cu-Cu interactions (parametrized in the file ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/Files/CU.pot|Files/CU.pot]]''), 
-a “C6/R^6 + pauli repulsion” will be defined for C-Cu and H-Cu and, finally,+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. 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. The charge of every specie (zero in our example) as well as the interaction within each possible type-pair has to be defined.
Line 280: Line 288:
 ==== Tasks ==== ==== 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. +  - 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 (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/INI/geo.inp|INI/geo.inp]]''/''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/FIN/geo.inp|FIN/geo.inp]]'') the output files (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/INI/geo.inp|INI/out]]''/''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/FIN/geo.inp|FIN/out]]'') and the trajectory of the geometry optimization (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/INI/INI-pos-1.xyz|INI-pos-1.xyz]]'' and ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/FIN/FIN-pos-1.xyz|FIN-pos-1.xyz]]''). 
-  - 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)+    - 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 ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/INI/ini_eq.xyz|ini_eq.xyz]]'' and ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/FIN/fin_eq.xyz|fin_eq.xzy]]''. 
 +    - 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 3 Å) to the H2 equilibrium distance. Such a guess is given, in the directory ''EX2'', with the files ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/ini_eq.xyz|ini_eq.xyz]]'', ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/due.xyz|due.xyz]]'',''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/tre.xyz|tre.xyz]]'',…,''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/tredici.xyz|tredici.xyz]]'',''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/fin_eq.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 (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/neb.inp|neb.inp]]''), the output (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX2/out|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: The last section of the input described above is modified as follows:
Line 347: Line 357:
 <note important>''[[inp>FORCE_EVAL/DFT/SCF#SCF_GUESS|SCF_GUESS RESTART]]''  in combination with  ''[[inp>FORCE_EVAL/DFT/QS#EXTRAPOLATION|EXTRAPOLATION USE_GUESS]]'' has to be declared in the ''[[inp>FORCE_EVAL/DFT]]'' section to allow for an efficient guess of the wavefunctions at each step of  the BAND optimization.</note> <note important>''[[inp>FORCE_EVAL/DFT/SCF#SCF_GUESS|SCF_GUESS RESTART]]''  in combination with  ''[[inp>FORCE_EVAL/DFT/QS#EXTRAPOLATION|EXTRAPOLATION USE_GUESS]]'' has to be declared in the ''[[inp>FORCE_EVAL/DFT]]'' section to allow for an efficient guess of the wavefunctions at each step of  the BAND optimization.</note>
  
 +{{ :exercises:2015_cecam_tutorial:distance_energy.jpg?direct&400 |}}
 ===== Exercise 2: Metadynamics ===== ===== Exercise 2: Metadynamics =====
  
 +The files for this part of the exercise are found in the **''EX1''** folder.
 +
 +==== Part 1 ====
 +
 +Obtain two restart files for independent configurations of the system at the geometry of the initial state thermalized at 450 K.
 +In the directory ''EQUIL'' the input file ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/EQUIL/inp|inp]]'' is used to run 3000 MD steps of  MD in the  NVT ensamble with a CSVR thermostat (a longer simulation should be used) two restart files obtained during the simulation (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/EQUIL/EQUIL-1_2000.restart|EQUIL-1_2000.restart]]'' and ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/EQUIL/EQUIL-1_3000.restart|EQUIL-1_3000.restart]]'') will be used to run a Metadynamics simulation with two WALKERS.
 +The trajectory of the MD should be also used to define the shape of the gaussians added in the metadynamics run.
 +Trajectories are printed in dcd format (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/EQUIL/EQUIL-pos-1.dcd|EQUIL-pos-1.dcd]]'' use, e.g. ''vmd [[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/EQUIL/s.xyz|s.xyz]] [[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/EQUIL/EQUIL-pos-1.dcd|EQUIL-pos-1.dcd]]'' to visualize the trajectory).
 +
 +The ''[[inp>MOTION]]'' section will look like:
 +
 +<code>
 +&MOTION
 +  &MD
 +    ENSEMBLE NVT
 +    STEPS    3000
 +    TIMESTEP 0.5
 +    TEMPERATURE 450.0
 +    &THERMOSTAT
 +       TYPE CSVR
 +       REGION MASSIVE
 +       &CSVR
 +          TIMECON [fs] 200.0
 +       &END
 +    &END
 +  &END MD
 +  &PRINT
 +    &TRAJECTORY
 +      FORMAT DCD
 +      &EACH
 +        MD 100
 +      &END
 +    &END
 +    &RESTART
 +      &EACH
 +        MD 1000
 +      &END
 +      FILENAME EQUIL
 +    &END
 +    &RESTART_HISTORY
 +     &EACH
 +      MD 1000
 +     &END
 +    &END
 +  &END
 +&END
 +</code>
 +
 +
 +==== Part 2 ====
 +
 +The two restart files are copied in the directory ''Files'' (''r1'' and ''r2'' respectively).
 +In the directory also the geometry for the whole system in the initial configuration as well as the geometry of the molecular fragment are included (files ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/Files/s|s]]'' and ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/Files/f|f]]'' respectively).
 +
 +We perform a metadynamics simulation with two collective variables: the sum of C-H distances (CV1) and the  H-H distance (CV2).
 +Walls have to be introduced to limit the range of variability of the two CVs.
 +
 +{{ :exercises:2015_cecam_tutorial:distances.jpg?direct&400 |}}
 +
 +To define the collective variables introduce a ''[[inp>FORCE_EVAL/SUBSYS/COLVAR|COLVAR]]'' card in the ''[[inp>FORCE_EVAL/SUBSYS|SUBSYS]]'' section of the first ''[[inp>FORCE_EVAL]]'':
 +
 +<code>
 +  &SUBSYS
 +    &COLVAR
 +       &COMBINE_COLVAR
 +          &COLVAR
 +             &DISTANCE
 +                ATOMS 4 50
 +             &END
 +          &END COLVAR
 +          &COLVAR
 +             &DISTANCE
 +                ATOMS 7 49
 +             &END
 +          &END COLVAR
 +          FUNCTION CV1+CV2
 +          VARIABLES CV1 CV2
 +          ERROR_LIMIT 1.0E-8
 +       &END
 +    &END
 +    &COLVAR
 +     &DISTANCE
 +      ATOMS 49 50
 +     &END
 +    &END
 +</code>
 +
 +The ''[[inp>MOTION]]'' section has to be modified as follows to activate metadynamics on the two CVs:
 +
 +<code>
 +&MOTION
 +  &MD
 +    ENSEMBLE NVT
 +    STEPS    50000
 +    TIMESTEP 0.5
 +    TEMPERATURE 450.0
 +    &THERMOSTAT
 +       TYPE CSVR
 +       REGION MASSIVE
 +       &CSVR
 +          TIMECON [fs] 200.0
 +       &END
 +    &END
 +  &END MD
 +  &FREE_ENERGY
 +    &METADYN
 +      DO_HILLS
 +      NT_HILLS 50
 +      WW 2.0e-3
 +      &METAVAR
 +        SCALE 0.1
 +        COLVAR 1
 +        &WALL
 +          TYPE QUADRATIC
 +          POSITION [angstrom] 5
 +          &QUADRATIC
 +             DIRECTION WALL_PLUS
 +             K [kcalmol] 40.0
 +          &END
 +        &END
 +        &WALL
 +          TYPE QUADRATIC
 +          POSITION [angstrom] 1.90
 +          &QUADRATIC
 +             DIRECTION WALL_MINUS
 +             K [kcalmol] 40.0
 +          &END
 +        &END
 +      &END METAVAR
 +      &METAVAR
 +        SCALE 0.1
 +        COLVAR 2
 +        &WALL
 +          TYPE QUADRATIC
 +          POSITION [angstrom] 3.80
 +          &QUADRATIC
 +             DIRECTION WALL_PLUS
 +             K [kcalmol] 40.0
 +          &END
 +        &END
 +        &WALL
 +          TYPE QUADRATIC
 +          POSITION [angstrom] 0.75
 +          &QUADRATIC
 +             DIRECTION WALL_MINUS
 +             K [kcalmol] 40.0
 +          &END
 +        &END
 +      &END METAVAR
 +#     &METAVAR
 +#       SCALE 0.1
 +#       COLVAR 3
 +#       &WALL
 +#         TYPE QUADRATIC
 +#         POSITION [angstrom] 4.00
 +#         &QUADRATIC
 +#            DIRECTION WALL_PLUS
 +#            K [kcalmol] 40.0
 +#         &END
 +#       &END
 +#       &WALL
 +#         TYPE QUADRATIC
 +#         POSITION [angstrom] 1.30
 +#         &QUADRATIC
 +#            DIRECTION WALL_MINUS
 +#            K [kcalmol] 40.0
 +#         &END
 +#       &END
 +#     &END METAVAR
 +      &PRINT
 +        &COLVAR
 +         COMMON_ITERATION_LEVELS 4
 +        &END
 +      &END
 +      &MULTIPLE_WALKERS
 +        NUMBER_OF_WALKERS 2
 +        WALKER_ID 1
 +        WALKER_COMM_FREQUENCY 4
 +        &WALKERS_FILE_NAME
 +           ../WALK_DATA_FILES/WALKER_1.data
 +           ../WALK_DATA_FILES/WALKER_2.data
 +        &END
 +      &END
 +    &END METADYN
 +  &END
 +  &PRINT
 +    &TRAJECTORY
 +      FORMAT DCD
 +      &EACH
 +        MD 100
 +      &END
 +    &END
 +    &RESTART
 +      &EACH
 +        MD 100
 +      &END
 +      FILENAME META
 +    &END
 +    &VELOCITIES ON
 +      &EACH
 +        MD 5000
 +      &END
 +    &END
 +    &RESTART_HISTORY OFF
 +    &END
 +  &END
 +&END
 +&EXT_RESTART
 + RESTART_FILE_NAME ../../Files/r1
 +&END
 +</code>
 +
 +The section specific to the multiple walkers could be removed to run the simulation with a single walker:
 +
 +<code>
 +      &MULTIPLE_WALKERS
 +        NUMBER_OF_WALKERS 2
 +        WALKER_ID 1
 +        WALKER_COMM_FREQUENCY 4
 +        &WALKERS_FILE_NAME
 +           ../WALK_DATA_FILES/WALKER_1.data
 +           ../WALK_DATA_FILES/WALKER_2.data
 +        &END
 +      &END
 +</code>
 +
 +Here we declare that we are using two walkers, that the input belongs to the walker with ID 1, that the data relative to the different walkers are stored in the files contained in the directory ''WALK_DATA_FILES'' and that the history of the //gaussian hills deposited// by each walker should be exchanged every 4 hills deposited.
 +
 +We need two independent directories (''WALK1'' and ''WALK2'' in the example) belonging to each WALKER, each directory will contain an input in which the WALKER ID is modified (1 and 2).
 +
 +The calculation will run as a FARMING: we need a common input that specifies that we are running two cp2k problems:
 +
 +The input file for the calculation (''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/inp|inp]]'' in the directory ''EX1'') will be:
 +
 +<code>
 +&GLOBAL
 +  PROGRAM FARMING
 +  RUN_TYPE NONE
 +&END GLOBAL
 +
 +&FARMING
 +  NGROUP 2
 +  &JOB
 +    DIRECTORY WALK1
 +    INPUT_FILE_NAME inp
 +  &END JOB
 +  &JOB
 +    DIRECTORY WALK2
 +    INPUT_FILE_NAME inp
 +  &END JOB
 +&END FARMING
 +</code>
 +
 +We declare that the cpus used for the calculation belongs to two groups of calculations,
 +one job has to run in the directory ''WALK1'' according to the prescriptions contained in the input file ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/WALK1/inp|WALK1/inp]]''
 +and a second job has to run in the directory ''WALK2'' (using ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/WALK2/inp|WALK2/inp]]'').
 +
 +If we run the calculation with 8 cpus, 4 will be devoted to the JOB in directory ''WALK1'' and 4 to ''WALK2''.
 +Each group of 10 cpus will be then partitioned according to partition instruction contained in each ''inp'' file (GROUP_PARTITION 1 3).
 +
 +During execution you can monitor the evolution of the CVs of each walker from the files 
 +''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/WALK1/META-COLVAR.metadynLog|META-COLVAR.metadynLog]]'' and ''[[http://cp2k.org/static/exercises/2015_cecam_tutorial/NEB/EX1/WALK2/META-COLVAR.metadynLog|META-COLVAR.metadynLog]]''
 +
 +e.g.
 +<code>
 +plot  'WALK1/META-COLVAR.metadynLog'   u 2:3 pt 7.0 ps 1.0 t "",\
 +      'WALK2/META-COLVAR.metadynLog'   u 2:3 pt 7.0 ps 1.0 t ""
 +</code>
 +
 +You can reconstruct the free energy profile from a restart file with the command:
 +
 +<code>
 +fes.popt -cp2k -mathlab -file META-META-1.restart -ndim 2 -ndw 1 2 -out fes_12.dat
 +</code>
 +
 +{{ :exercises:2015_cecam_tutorial:fes.jpg?direct&400 |}}
 +===== References =====
  
 +~~REFNOTES~~
exercises/2015_cecam_tutorial/neb.1440063476.txt.gz · Last modified: 2020/08/21 10:14 (external edit)