Processing math: 100%

CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2017_ethz_mmm:reaction_energy_2017

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:2017_ethz_mmm:reaction_energy_2017 [2017/04/27 09:05] – dpasseroneexercises:2017_ethz_mmm:reaction_energy_2017 [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 1: Line 1:
-====== Bond energy of ethylene ====== +====== Dehydration of ethanol ====== 
-<!-- +We will investigate today a very important chemical reaction, the production of ethene (ethylene) from ethanol. Ethanol is heated with an excess of concentrated sulphuric acid at a temperature of 170°C. The gases produced are passed through sodium hydroxide solution to remove the carbon dioxide and sulphur dioxide produced from side reactions. 
-===== Introduction ===== +Another way to favor this reaction is in presence of a catalyst surface, which makes the reaction exothermic even at room temperature, and the free energy barrier for the reaction is highly reduced, so to obtain reasonable reaction rates. 
-In this tutorial, we are going to show the reader how to perform a +The use of hematite (iron-based nanocrystal) as a catalyst is described computationally in the following paper by Lopes et al: [[doi>10.1039/c6ra08509a]]. We are interested in the gas phase reaction, in particular in the reactants and the products energy, to estimate the "energy of reaction".
-simple static self-consistent Kohn-Sham Density Functional Theory +
-energy and force calculation on a system using ''QUICKSTEP''.+
  
-We will use face centred cubic bulk Si, with 8 atoms in a cubic unit +The reaction is the following: 
-cell as an example. The example files are contained in {{:static_calculation.tgz|static_calculation.tgz}}. The calculations were carried out with +{{ :exercises:2017_ethz_mmm:screen_shot_2017-04-27_at_18.43.25.png?nolink&600 |}} 
-''CP2K'' version 2.4.+and this is the pictorial view in the gas phase: 
 +{{ :exercises:2017_ethz_mmm:screen_shot_2017-04-27_at_18.50.10.png?nolink&600 |}}
  
-The reader should note that the integration grid settings used in +In the supplementary material of the paper, which can be found [[http://www.rsc.org/suppdata/c6/ra/c6ra08509a/c6ra08509a1.pdf|here]], several calculation with different methods ranging from DFT to hybrid methods to high-level methods are described, together with the basis set dependence of the results.  
-the example calculations have already chosen to be sufficient for +
-the given accuracy. A [[converging_cutoff|separate tutorial]] is available to show how to +
-achieve this.+
  
  
-===== Input Files ===== +We will compare our results with the published ones.
-We first look at the input files required for this calculation. The +
-necessary input files are:+
  
-  * ''Si_bulk8.inp'': the main input file for the calculation, which defines the system and the job parameters +The input file structure is shown below:
-  * ''BASIS_SET'': file containing parameters for the basis sets that can be used by ''CP2K'' for this calculation +
-  * ''POTENTIAL'': file containing parameters for the pseudopotentials that can be used by ''CP2K'' for this calculation+
  
-A list of basis set and pseudopotential files may be found in +<code cp2k>
-''cp2k/data/'' that comes with a ''cp2k'' source release. These +
-should cover most of the commonly used elements. The user will, +
-however, need to produce their own main input file for a given +
-calculation.+
  
-Let us look at the main input: ''Si_bulk8.inp'' in detail. The name of +   &DFT 
-this file can be arbitrary, so is the file extension. The input file +      &POISSON                        ! Solver requested for non periodic calculations 
-is structured into ordered blocks and keywords, the order of which +         PERIODIC NONE 
-are unimportant. Each input block is referred to as a "section" in +         PSOLVER  WAVELET             ! Type of solver 
-this tutorial, and some sections are "subsections" of other +      &END POISSON 
-sections. Full definitions of the input file format and keywords is +      &QS                             ! Parameters needed to set up the Quickstep framework 
-available via the [[http://manual.cp2k.org/trunk/index.html|''CP2K'' input reference manual]].+         METHOD GAPW                  ! Method: gaussian and augmented plane waves 
 +      &END QS
  
-The input file is shown below:+# 
 +# Include the exchange and correlation information 
 +# 
 +@INCLUDE './pbe.inc'
  
-<code cp2k> +   &END DFT
-&GLOBAL +
-  PROJECT Si_bulk8 +
-  RUN_TYPE ENERGY_FORCE +
-  PRINT_LEVEL LOW +
-&END GLOBAL +
-&FORCE_EVAL +
-  METHOD Quickstep +
-  &SUBSYS +
-    &KIND Si +
-      ELEMENT   Si +
-      BASIS_SET DZVP-GTH-PADE +
-      POTENTIAL GTH-PADE-q4 +
-    &END KIND +
-    &CELL +
-      A     5.430697500    0.000000000    0.000000000 +
-      B     0.000000000    5.430697500    0.000000000 +
-      C     0.000000000    0.000000000    5.430697500 +
-    &END CELL +
-    &COORD +
-      Si    0.000000000    0.000000000    0.000000000 +
-      Si    0.000000000    2.715348700    2.715348700 +
-      Si    2.715348700    2.715348700    0.000000000 +
-      Si    2.715348700    0.000000000    2.715348700 +
-      Si    4.073023100    1.357674400    4.073023100 +
-      Si    1.357674400    1.357674400    1.357674400 +
-      Si    1.357674400    4.073023100    4.073023100 +
-      Si    4.073023100    4.073023100    1.357674400 +
-    &END COORD +
-  &END SUBSYS +
-  &DFT +
-    BASIS_SET_FILE_NAME  BASIS_SET +
-    POTENTIAL_FILE_NAME  GTH_POTENTIALS +
-    &QS +
-      EPS_DEFAULT 1.0E-10 +
-    &END QS +
-    &MGRID +
-      NGRIDS 4 +
-      CUTOFF 300 +
-      REL_CUTOFF 60 +
-    &END MGRID +
-    &XC +
-      &XC_FUNCTIONAL PADE +
-      &END XC_FUNCTIONAL +
-    &END XC +
-    &SCF +
-      SCF_GUESS ATOMIC +
-      EPS_SCF 1.0E-7 +
-      MAX_SCF 300 +
-      &DIAGONALIZATION  ON +
-        ALGORITHM STANDARD +
-      &END DIAGONALIZATION +
-      &MIXING  T +
-        METHOD BROYDEN_MIXING +
-        ALPHA 0.4 +
-        NBROYDEN 8 +
-      &END MIXING +
-    &END SCF +
-  &END DFT +
-  &PRINT +
-    &FORCES ON +
-    &END FORCES +
-  &END PRINT +
-&END FORCE_EVAL +
-</code>+
  
-The main sections in the input file are: 
  
-  * [[http://manual.cp2k.org/trunk/CP2K_INPUT/GLOBAL.html|''GLOBAL'']]: contains general options for the ''CP2K'' run, such as the name of the job, the type of run etc. +   &SUBSYS                            ! This section defines the system 
-  * [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL.html|''FORCE_EVAL'']]: contains all parameters associated with the evaluation of forces on atoms, this includes the initial atomic coordinates.+      &CELL                           ! Unit cell set up 
 +         PERIODIC NONE                ! Non periodic calculation 
 +         ABC 10 10 10                 ! Lengths of the cell vectors A, B, and C 
 +      &END CELL 
 +      &COORD                          ! This section specify all the atoms and their coordinates 
 +  H         2.5558925119        4.9201502245        3.9999999733 
 +  H         2.5617214834        3.0599837703        4.0000000267 
 +  C         3.3235178526        3.9916405319        4.0000000002 
 +  C         4.6493228423        3.9951429630        3.9999999995 
 +  H         5.4118318595        4.9261779106        4.0000000269 
 +  H         5.4167134504        3.0669045998        3.9999999735 
 +      &END COORD 
 +      &TOPOLOGY 
 +          &CENTER_COORDINATES 
 +          &END 
 +      &END TOPOLOGY
  
-We look at each section in detail. The ''GLOBAL'' section in +      &KIND H                                      ! potential and basis for H 
-''Si_bulk8.inp'' is presented below:+         &BASIS 
 +@INCLUDE './H_631Gdp.inc' 
 +         &END 
 +         POTENTIAL ALL 
 +         &POTENTIAL 
 +          1    0    0 
 +          0.20000000    0 
 +         &END 
 +      &END KIND 
 +      &KIND C                         ! potential and basis for C 
 +         &BASIS 
 +@INCLUDE './C_631Gdp.inc' 
 +         &END 
 +         POTENTIAL ALL 
 +         &POTENTIAL 
 +            4    2    0 
 +            0.34883045    0 
 +         &END 
 +      &END KIND 
 +   &END SUBSYS 
 +&END FORCE_EVAL                     ! This section defines method for calculating energy and forces
  
-<code cp2k> +&GLOBAL                             ! Section with general information regarding which kind of simulation to perform an parameters for the whole PROGRAM 
-&GLOBAL +   PRINT_LEVEL LOW                  ! Global print level 
-  PROJECT Si_bulk8 +   PROJECT c2h4                     ! Name of the project. This word will appear as part of a name of all ouput files (except main ouput file, specified with -o option) 
-  RUN_TYPE ENERGY_FORCE +   RUN_TYPE GEO_OPT                 ! Geometry optimization
-  PRINT_LEVEL LOW+
 &END GLOBAL &END GLOBAL
 +                                                                                                       
 </code> </code>
  
-We will be doing a static energy and force calculation, in this +The exercise is executed on the ''hypatia'' cluster. The instructions to connect are found [[exercises:2017_ethz_mmm:replica_2017|here]]. After creating a directory as usual, you can copy the files from the following location: 
-case, we must set [[http://manual.cp2k.org/trunk/CP2K_INPUT/GLOBAL.html#desc_RUN_TYPE|''RUNTYPE'']] to ''ENERGY_FORCE''. Keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/GLOBAL.html#desc_PROJECT_NAME|''PROJECT'']] is an +<note tip>~psd/Exercise_8/</note>
-alias for ''PROJECT_NAME'', which sets the root-name of the +
-calculation, in this case ''Si_bulk8''. Any output files automatically +
-generated by ''CP2K'' will have the name prefixed by +
-''Si_bulk8''. [[http://manual.cp2k.org/trunk/CP2K_INPUT/GLOBAL.html#desc_PRINT_LEVEL|''PRINT_LEVEL'']] controls the default verbosity of the main +
-output of ''CP2K'', in this example, it is set to "low". The +
-verbosity of the output can be fine-tuned by overriding this setting +
-in each individual subsection of the input.+
  
-We now explain the section ''FORCE_EVAL'' line-by-line.+Copy the files to the created directory in the ''/mnt/project/yourusername/'' path. 
 +<note tip> 
 +In the directory you will find the following files: 
 +  * ''h2o.inp'' for the geometry optimization of water 
 +  * ''ethanol.inp'' for the geometry optimization of ethanol 
 +  * ''ethylene.inp'' for the geometry optimization of ethylene (ethene) 
 +  * ''run'' to launch the jobs to the queue 
 +  * several ''*.inc'' file describing the level of theory (exchange and/or correlation functional) and the basis sets. 
 +</note>
  
-<code cp2k> +The command to launch the job is 
-METHOD Quickstep +<note important>> qsub run -v INP=file </note> 
-</code>+where file.inp has to be replaced by the relevant prefix of the input file (example: ''h2o''). 
 +In the ''*.xyz'' files you can look for the final energies by the following command: 
 +<note important>> grep i c2h6o-pos-1.xyz  </note>. 
 +This will list all the energies. Pick the last one for the optimized structure. You can also visualize it with vmd.
  
-The keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL.html#desc_METHOD|''METHOD'']] chooses the method for evaluating the forces on +You will run the calculation by changing in **ALL THREE INPUT FILES** the ''@INCLUDE'' lines, in the following way: 
-atoms to ''QUICKSTEP'', i.e. Density Functional Theory using the +At the line concerning the exchange and correlation potential, you prescribe the level of theory: 
-Gaussian and Planewaves (GPW) method.+<note important>   
 +  * hf.inc for Hartree-Fock 
 +  * pbe.inc for DFT/PBE, with gradient corrected local density 
 +  * b3lyp.inp for hybrid functional containing a percentage of exact exchange 
 +</note>
  
-<code cp2k> +When you are done with the three levels of theory, then you can redo the exercise with the larger basis set, and change EVERYWHERE in the input file the instances of ''H_631Gdp.inc'' to ''H_6311++G2d2p.inc'' and the same for the O_* and C_* instances. 
-&SUBSYS +Redo the calculations with the three levels ''pbe'', ''b3lyp'' and ''hf''. 
-  &KIND Si +  
-    ELEMENT   Si +<note warning> 
-    BASIS_SET DZVP-GTH-PADE +Assignments: 
-    POTENTIAL GTH-PADE-q4 +  - Compute the reaction energy for the dehydration reaction of ethanol 
-  &END KIND +  - Prepare a table with rows and columns: on the rows the level of theory, on the columns the basis set (3x2 table) 
-  &CELL +  - Compare the results with the published ones (note the conversion factors. You can use the tool at [[http://www.colby.edu/chemistry/PChem/Hartree.html|Energy converter]] 
-    A     5.430697500    0.000000000    0.000000000 +  - Comment on the dependence on the basis set and on the level of theory (hint: this also need the next theory lecture) 
-    B     0.000000000    5.430697500    0.000000000 +  - Is this information enough to determine the rates of reaction? Why? 
-    C     0.000000000    0.000000000    5.430697500 +</note>
-  &END CELL +
-  &COORD +
-    Si    0.000000000    0.000000000    0.000000000 +
-    Si    0.000000000    2.715348700    2.715348700 +
-    Si    2.715348700    2.715348700    0.000000000 +
-    Si    2.715348700    0.000000000    2.715348700 +
-    Si    4.073023100    1.357674400    4.073023100 +
-    Si    1.357674400    1.357674400    1.357674400 +
-    Si    1.357674400    4.073023100    4.073023100 +
-    Si    4.073023100    4.073023100    1.357674400 +
-  &END COORD +
-&END SUBSYS +
-</code>+
  
-The subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS.html|''SUBSYS'']] defines the simulation unit cell and the 
-initial coordinates of atoms in the calculation.  
  
-The subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/KIND.html|''KIND'']] gives definitions of elements in the +=== BONUS TRACK ===
-calculation. There must be one ''KIND'' subsection per element. In +
-this example, for Si, we have defined the basis set to be used: +
-''DZVP-GTH-PADE'' (double-ζ with polarisation basis optimised +
-for Geodecker-Teter-Hutter PADE LDA pseudopotential); and the +
-pseudopotential: ''GTH-PADE-q4'' (Geodecker-Teter-Hutter PADE LDA +
-pseudopotential with 4 valence electrons).+
  
-The basis set and pseudopotential names //must// correspond to an +<note tip>We may be interested in the visualisation of the electronic density. Copy the ''ethanol.inp'' into ''ethanol_dens.inp''.</note>
-existing entry in the corresponding basis set and pseudopotential +
-files defined by [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html#desc_BASIS_SET_FILE_NAME|''BASIS_SET_FILE_NAME'']] and [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html#desc_POTENTIAL_FILE_NAME|''POTENTIAL_FILE_NAME'']] +
-keywords in [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html|''DFT'']] subsection, in ''FORCE_EVAL'' section. The chosen +
-basis for Si corresponds to parameters: +
- +
-<code> +
-Si DZVP-GTH-PADE +
-  2 +
-  3  0  1  4  2  2 +
-        1.2032422345   0.3290350445   0.0000000000   0.0474539126   0.0000000000 +
-        0.4688409786  -0.2533118323   0.0000000000  -0.2594473573   0.0000000000 +
-        0.1679863234  -0.7870946277   0.0000000000  -0.5440929303   0.0000000000 +
-        0.0575619526  -0.1909898479   1.0000000000  -0.3624010364   1.0000000000 +
-  3  2  2  1  1 +
-        0.4500000000   1.0000000000 +
-</code> +
- +
-in file ''BASIS_SET''; and the chosen pseudopotential corresponds to +
-parameters: +
- +
-<code> +
-Si GTH-PADE-q4 GTH-LDA-q4 +
-    2    2 +
-     0.44000000    1    -7.33610297 +
-    2 +
-     0.42273813    2     5.90692831    -1.26189397 +
-                                        3.25819622 +
-     0.48427842    1     2.72701346 +
-</code> +
- +
-in file ''GTH_POTENTIALS''. +
- +
-The subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/CELL.html|''CELL'']] defines the simulation unit cell used in a +
-calculation.  In this example, we define the unit cell as cubic, +
-with lattice constant equal to 5.4306975 Angstroms. "Angstrom" is +
-the default unit for cell vectors. [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/CELL.html#desc_A|''A'']], [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/CELL.html#desc_B|''B'']] and [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/CELL.html#desc_C|''C'']] are the first, +
-second and third lattice (cell) vectors. There are many ways to +
-define the cell, see [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/CELL.html|''CP2K'' input reference manual]] for more details. +
- +
-The initial atomic coordinates are specified in the [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/COORD.html|''COORD'']] +
-subsection. The default input format for atomic coordinates in +
-''CP2K'' is: +
-<code> +
-<ATOM_KIND> X Y Z +
-</code> +
-where ''X'', ''Y'' and ''Z'' are Cartesian coordinates in Angstroms. This +
-can be changed by configuring keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/COORD.html#desc_SCALED|''SCALED'']] to ''.TRUE.'', in the +
-''COORD'' subsection, which makes the coordinate input ''X'' ''Y'' ''Z'' to +
-be fractional with respect to the lattice vectors. One can also +
-change the unit for the Cartesian coordinates by setting the keyword +
-[[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/COORD.html#desc_UNIT|''UNIT'']] with in the subsection. ''<ATOM_KIND>'' should be a label that +
-corresponds to the definition of the elements in the [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/KIND.html|''KIND'']] +
-subsections. +
- +
-After the ''SUBSYS'' section in the input file ''Si_bulk8.inp'' follows +
-the [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html|''DFT'']] subsection, which controls all aspects of the +
-self-consistent Kohn-Sham Density Functional Theory +
-calculation. This subsection is only relevant if and only if the +
-''METHOD'' keyword in ''FORCE_EVAL'' is set to ''QUICKSTEP''.+
  
 +Add the following sections:
 +**under &DFT**
 <code cp2k> <code cp2k>
-BASIS_SET_FILE_NAME  BASIS_SET +       &PRINT  
-POTENTIAL_FILE_NAME  GTH_POTENTIALS+          &E_DENSITY_CUBE 
 +          &END 
 +       &END       
 +       &SCF 
 +          SCF_GUESS RESTART 
 +       &END
 </code> </code>
 +This tells to read the old wavefunction and to print the cubefile of the density.
  
-As already mentioned above, the keywords [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html#desc_BASIS_SET_FILE_NAME|''BASIS_SET_FILE_NAME'']] and +At the end of the input file:
-[[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html#desc_POTENTIAL_FILE_NAME|''POTENTIAL_FILE_NAME'']] set the files that contains basis set and +
-pseudopotential parameters. +
 <code cp2k> <code cp2k>
-&QS +&EXT_RESTART 
-  EPS_DEFAULT 1.0E-10 +     RESTART_FILE_NAME ./c2h6o-1.restart 
-&END QS+&END
 </code> </code>
  
-The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/QS.html|''QS'']] subsection contains general control parameters used by +Then, change ''RUN_TYPE GEO_OPT'' to ''RUN_TYPE ENERGY'' to only run a single point calculation. It will generate a cubefile with the density which you may visualize with VMD.
-''QUICKSTEP''. [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/QS.html#desc_EPS_DEFAULT|''EPS_DEFAULT'']] sets the default value for all tolerances +
-used within ''QUICKSTEP''. The individual tolerances (''EPS_*'') can be +
-set, and they will override the ''EPS_DEFAULT'' value. +
- +
-<code cp2k> +
-&MGRID +
-  NGRIDS 4 +
-  CUTOFF 300 +
-  REL_CUTOFF 60 +
-&END MGRID +
-</code> +
- +
-The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/MGRID.html|''MGRID'']] subsection defines how the integration grid used in +
-''QUICKSTEP'' calculations should be setup. ''QUICKSTEP'' uses a +
-multi-grid method for representing Gaussian functions numerically on +
-the grid. Narrow and sharp Gaussians are mapped onto a finer grid +
-than wider and smoother Gaussians. In this case, we are telling the +
-code to set up 4 levels of multi-grids, with the planewave cutoff of +
-the finest grid set to be 300 Ry, and with the grid spacing +
-underneath any Gaussian functions to be finer than the equivalent +
-planewave cutoff of 60 Ry. The users should read the tutorial +
-"[[converging_cutoff|Converging the CUTOFF and REL_CUTOFF]]" for details on how these +
-parameters affect the grid constructed, and how to define a +
-sufficient grid for their calculation. In this example, the grid +
-defined has already been found to be sufficient for the energy and +
-force calculation. +
- +
-The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/XC.html|''XC'']] subsection follows: +
- +
-<code cp2k> +
-&XC +
-  &XC_FUNCTIONAL PADE +
-  &END XC_FUNCTIONAL +
-&END XC +
-</code> +
- +
-This defines which exchange-correlation density functional we want +
-to use. In this we choose PADE LDA functional, which is consistent +
-with the basis set and pseudopotential we have chosen. +
- +
-<code cp2k> +
-&SCF +
-  SCF_GUESS ATOMIC +
-  EPS_SCF 1.0E-7 +
-  MAX_SCF 300 +
-  &DIAGONALIZATION +
-    ALGORITHM STANDARD +
-  &END DIAGONALIZATION +
-  &MIXING +
-    METHOD BROYDEN_MIXING +
-    ALPHA 0.4 +
-    NBROYDEN 8 +
-  &END MIXING +
-&END SCF +
-</code> +
- +
-The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF.html|''SCF'']] subsection defines all the settings related to methods +
-used to find a self-consistent solution of the Kohn-Sham DFT +
-formalism. +
- +
-[[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF.html#desc_SCF_GUESS|''SCF_GUESS'']] sets how the initial trial electron density function +
-ρ(→r) is to be generated. In this example (''ATOMIC''), the +
-initial density is to be generated using overlapping of atomic +
-charge densities. A good starting point for the electron density in +
-the self-consistency loop is important in obtaining a convergent +
-result quickly.  [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF.html#desc_EPS_SCF|''EPS_SCF'']] sets the tolerance of the charge density +
-residual. This overrides the ''EPS_DEFAULT'' value set in ''QS'' +
-subsection. [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF.html#desc_MAX_SCF|''MAX_SCF'']] sets the maximum number of self-consistency +
-loops ''QUICKSTEP'' is allowed to perform for each ground-state energy +
-calculation. +
- +
-<code cp2k> +
-&DIAGONALIZATION  ON +
-  ALGORITHM STANDARD +
-&END DIAGONALIZATION +
-</code> +
- +
-The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/DIAGONALIZATION.html|''DIAGONALIZATION '']] subsection tells the code to use the +
-traditional diagonalisation method for finding the ground state +
-Kohn-Sham energy and electron density. The subsection heading also +
-takes an argument, and in this case is set to "''ON''", which +
-equivalent to "''.TRUE.''" or "''T''", and indicates that the +
-diagonalisation method is turned on. One can also omit the value of +
-the subsection heading, which defaults to "''.TRUE.''". The +
-alternative to diagonalisation is to use the Orbital Transform (OT) +
-method, in which case, the user should either delete the +
-''DIAGONALIZATION'' block or change "''ON''" to "''OFF''" (or +
-"''.FALSE.''"), and add the [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/DIAGONALIZATION/OT.html|''OT'']] subsection instead. The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/DIAGONALIZATION.html#desc_ALGORITHM|''ALGORITHM'']] +
-keyword sets the algorithm to use for diagonalisation of the +
-Kohn-Sham Hamiltonian. "''STANDARD''" means the standard +
-LAPACK/SCALAPACK subroutines are to be used for diagonalisation. +
- +
-<code cp2k> +
-&MIXING  T +
-  METHOD BROYDEN_MIXING +
-  ALPHA 0.4 +
-  NBROYDEN 8 +
-&END MIXING +
-</code> +
- +
-The [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/MIXING.html|''MIXING'']] subsection contains all the parameters associated with +
-charge mixing in a self-consistency calculation. The subsection also +
-admits a value, which can be either ''.TRUE.'' (''T'') or ''.FALSE.'' +
-(''F''), which switches charge mixing on or off. The default is +
-''.TRUE.''. Note that this subsection //only applies to the traditional +
-diagonalisation method//. The OT method uses a different approach for +
-charge mixing, and is explained in other tutorials. The keyword +
-[[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/MIXING.html#desc_ALPHA|''ALPHA'']] sets the mixing parameter; in this example 0.4 of the output +
-density will be mixed with 0.6 of the input density to form the new +
-input density in the next SCF iteration. The keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/MIXING.html#desc_METHOD|''METHOD'']] sets +
-the mixing method; in this case, we will use Broyden mixing. The +
-keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/MIXING.html#desc_NBUFFER|''NBROYDEN'']] is an alias to the parameter ''NBUFFER'', and it +
-sets the number of histories to be used in the Broyden mixing +
-algorithm. +
- +
-The final [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PRINT.html|''PRINT'']] subsection in ''FORCE_EVAL'' section: +
- +
-<code cp2k> +
-&PRINT +
-  &FORCES ON +
-  &END FORCES +
-&END PRINT +
-</code> +
- +
-tells ''CP2K'', in this case, to print out atomic forces in the main +
-output of the calculation. +
- +
- +
-===== Running the Calculation ===== +
-To run the calculation, the reader needs to have a working ''CP2K'' +
-package compiled, and with the path to its binaries in the system +
-''PATH''. The files ''Si_bulk8.inp'', ''BASIS_SET'' and ''GTH_POTENTIAL'' +
-should be in the same working directory. In this example, we will +
-use the MPI version of ''CP2K''. Type command: +
-<code> +
-mpirun -n 2 cp2k.popt -o Si_bulk8.out Si_bulk8.inp & +
-</code> +
-in the working directory to run ''CP2K'' in parallel with 2 MPI +
-processes in the background. The ''-o'' option redirects the ''CP2K'' +
-output to file ''Si_bulk8.out''. Note that the ''-o'' option //appends// +
-output of successive runs to ''Si_bulk8.out'', so if the reader wants +
-to start completely from afresh, they must delete ''Si_bulk8.out'' +
-before running a new calculation. +
- +
- +
-===== Obtaining the Results ===== +
-After the job has finished, we should obtain the following files: +
- +
-  * ''Si_bulk8.out'' +
-  * ''Si_bulk8-RESTART.wfn'' +
-  * ''Si_bulk8-RESTART.wfn.bak-1'' +
- +
-The file ''Si_bulk8.out'' contains the main output of the +
-job. ''Si_bulk8-RESTART.wfn'' is the final Kohn-Sham wavefunctions +
-from the calculation. ''Si_bulk8-RESTART.wfn.bak-<n>'' records the +
-Kohn-Sham wavefunctions obtained from the ''<n>''-th previous SCF +
-step; in this case, ''Si_bulk8-RESTART.wfn.bak-1'' contains the +
-wavefunctions obtained from the last SCF step. Should the reader +
-want to start a new calculation from the previous calculated +
-wavefunctions, they need to change the ''SCF_GUESS'' keyword in ''SCF'' +
-subsection of ''FORCE_EVAL'' section to ''RESTART'': +
-<code cp2k> +
-SCF_GUESS RESTART +
-</code> +
-provided that the new calculation shares the same ''PROJECT_NAME'' as +
-the one that generated the wavefunctions. Otherwise, we would need +
-to rename the restart wavefunction files to correspond to the +
-project name of the new calculation. +
- +
-We now look at the main ''CP2K'' output in more detail. +
- +
-<code> +
-Number of electrons:                                                         32 +
-Number of occupied orbitals:                                                 16 +
-Number of molecular orbitals:                                                16 +
- +
-Number of orbital functions:                                                104 +
-Number of independent orbital functions:                                    104 +
- +
-Extrapolation method: initial_guess +
- +
- +
-SCF WAVEFUNCTION OPTIMIZATION +
- +
- Step     Update method      Time    Convergence         Total energy    Change +
- ------------------------------------------------------------------------------ +
-    1 NoMix/Diag. 0.40E+00    0.6     0.75558724       -32.2320848878 -3.22E+01 +
-    2 Broy./Diag. 0.40E+00    1.1     0.05667976       -31.1418135481  1.09E+00 +
-    3 Broy./Diag. 0.40E+00    1.1     0.09691469       -31.1974003416 -5.56E-02 +
-    4 Broy./Diag. 0.40E+00    1.1     0.00245608       -31.3378474040 -1.40E-01 +
-    5 Broy./Diag. 0.40E+00    1.1     0.00235460       -31.3009654398  3.69E-02 +
-    6 Broy./Diag. 0.40E+00    1.1     0.00007565       -31.2972158934  3.75E-03 +
-    7 Broy./Diag. 0.40E+00    1.1     0.00009004       -31.2977293749 -5.13E-04 +
-    8 Broy./Diag. 0.40E+00    1.1     0.00000186       -31.2978454163 -1.16E-04 +
-    9 Broy./Diag. 0.40E+00    1.1     0.00000252       -31.2978835492 -3.81E-05 +
-   10 Broy./Diag. 0.40E+00    1.1     5.6405E-09       -31.2978852054 -1.66E-06 +
- +
- *** SCF run converged in    10 steps *** +
-</code> +
- +
-The above shows a typical output from a self-consistent Kohn-Sham +
-ground state calculation. It states that we are using +
-diagonalisation method with Broyden charge mixing, and it took 10 +
-Broyden mixing steps (each containing a diagonalisation process to +
-obtain the wavefunctions) to reach the required tolerance for +
-self-consistency. +
- +
-<code> +
- Electronic density on regular grids:        -31.9999999889        0.0000000111 +
- Core density on regular grids:               31.9999999939       -0.0000000061 +
- Total charge density on r-space grids:        0.0000000051 +
- Total charge density g-space grids:           0.0000000051 +
- +
- Overlap energy of the core charge distribution:               0.00000000005320 +
- Self energy of the core charge distribution:                -82.06393942512820 +
- Core Hamiltonian energy:                                     18.06858429706010 +
- Hartree energy:                                              42.41172824581682 +
- Exchange-correlation energy:                                 -9.71425832315952 +
- +
- Total energy:                                               -31.29788520535761 +
- +
-ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):              -31.297885372811002 +
- +
- +
-ATOMIC FORCES in [a.u.] +
- +
-# Atom   Kind   Element          X              Y              Z +
-     1      1      Si          0.00000000     0.00000000     0.00000000 +
-     2      1      Si          0.00000000     0.00000001     0.00000001 +
-     3      1      Si          0.00000001     0.00000001     0.00000000 +
-     4      1      Si          0.00000001     0.00000000     0.00000001 +
-     5      1      Si         -0.00000001    -0.00000001    -0.00000001 +
-     6      1      Si         -0.00000001    -0.00000001    -0.00000001 +
-     7      1      Si         -0.00000001    -0.00000001    -0.00000001 +
-     8      1      Si         -0.00000001    -0.00000001    -0.00000001 +
-SUM OF ATOMIC FORCES          -0.00000000    -0.00000000    -0.00000000     0.00000000 +
-</code> +
- +
-The above shows the results on final energies and forces.  One +
-should always check if the total number of electrons calculated from +
-the final electron density, in this case 31.9999999939, is correct. +
- +
-The results show that the force on the atoms are almost zero. This +
-means the system is more or less relaxed, and its geometry is close +
-to its optimal at ground state. +
- +
- +
-===== Adding Smearing ===== +
-In the example so far, we have not used any smearing on electron +
-occupation. This is fine for system with a large band gap. However, +
-for metals or systems with a small gap, this may cause the +
-calculation to be unstable, and the self-consistency loop may never +
-converge, due to the discontinuity in the electron occupation +
-function. +
- +
-To add smearing, we need to add the subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF/SMEAR.html|''SMEAR'']] inside +
-subsection ''SCF'': +
- +
-<code cp2k> +
-&SMEAR ON +
-  METHOD FERMI_DIRAC +
-  ELECTRONIC_TEMPERATURE [K] 300 +
-&END SMEAR +
-</code> +
- +
-This tells ''QUICKSTEP'' to use a smearing function for the electron +
-occupation. In this example, we use the Fermi-Dirac smearing +
-function, with the electron temperature being set to 300 K. Note +
-that, in ''CP2K'', one can explicitly define the unit of a given input +
-value by using a unit enclosed in square bracket, such as "''[K]''", +
-before the value. +
- +
-This is not all, since smearing may lead to occupation of molecular +
-orbitals in the conduction band, we must tell ''CP2K'' to include +
-extra, empty, molecular orbitals into the calculation, which +
-otherwise would be omitted (for reducing computational cost). To do +
-this, we need to add the keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/SCF.html#desc_ADDED_MOS|''ADDED_MOS'']] in the ''SCF'' +
-subsection: +
- +
-<code cp2k> +
-ADDED_MOS 10 +
-</code> +
- +
-In this example, we have asked ''CP2K'' to not to omit 10 of the +
-lowest empty molecular orbitals in the calculation. It should be +
-noted that given a chosen basis set, there is a maximum number of +
-molecular orbitals, i.e. the number of eigenvectors of the +
-Hamiltonian, one can generate. In theory, the maximum should be the +
-rank of the Hamiltonian generated in the calculation. +
- +
-In the output of the calculation using smearing, we first notice +
-that: +
- +
-<code> +
-Number of electrons:                                                         32 +
-Number of occupied orbitals:                                                 16 +
-Number of molecular orbitals:                                                26 +
- +
-Number of orbital functions:                                                104 +
-Number of independent orbital functions:                                    104 +
-</code> +
- +
-unlike in the previous case with no smearing, now 26 molecular +
-orbitals have been used during the calculation. There are a total of +
-104 basis functions ("atom centred orbitals" spanned by Cartesian +
-Gaussians) used in the calculation for the given basis set, so 26 is +
-well within the limit of the calculation. +
- +
-<code> +
- Electronic density on regular grids:        -31.9999999889        0.0000000111 +
- Core density on regular grids:               31.9999999939       -0.0000000061 +
- Total charge density on r-space grids:        0.0000000051 +
- Total charge density g-space grids:           0.0000000051 +
- +
- Overlap energy of the core charge distribution:               0.00000000005320 +
- Self energy of the core charge distribution:                -82.06393942512820 +
- Core Hamiltonian energy:                                     18.06842027191411 +
- Hartree energy:                                              42.41184371469986 +
- Exchange-correlation energy:                                 -9.71419454555998 +
- Electronic entropic energy:                                  -0.00001687947145 +
- Fermi energy:                                                 0.20867150262130 +
- +
- Total energy:                                               -31.29788686349247 +
- +
-ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):              -31.297887031736590 +
-</code> +
- +
-In the final energy section of the output, we notice that there is +
-an extra entropy (TS) term: +
- +
-<code> +
-Electronic entropic energy:                                  -0.00001687947145 +
-</code> +
- +
-This should be small for the calculation to be a reliable +
-approximation to the zero electron temperature result. The final +
-free energy is the sum of the total DFT energy and the entropic +
-energy. The total DFT energy is given by: +
- +
-<code> +
-Total energy:                                               -31.29788686349247 +
-</code> +
- +
-and the final free energy extrapolated for TS→0 is given by: +
- +
-<code> +
-ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):              -31.297887031736590 +
-</code>+
  
-This is the energy to be quoted as the final result. --!> 
exercises/2017_ethz_mmm/reaction_energy_2017.1493283959.txt.gz · Last modified: 2020/08/21 10:15 (external edit)