CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2017_uzh_cp2k-tutorial:hybrid

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_uzh_cp2k-tutorial:hybrid [2017/07/11 19:16] – gtocciexercises:2017_uzh_cp2k-tutorial:hybrid [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 3: Line 3:
 The purpose of this section is to explain how to perform hybrid functional calculations (or Hartree-Fock exchange, HFX) with CP2K in condensed phase systems. It is based on the developments described in [[doi>10.1021/ct900494g]] and [[doi>10.1063/1.2931945]], and its efficient extension (ADMM) described in [[doi>10.1021/ct1002225]]. The purpose of this section is to explain how to perform hybrid functional calculations (or Hartree-Fock exchange, HFX) with CP2K in condensed phase systems. It is based on the developments described in [[doi>10.1021/ct900494g]] and [[doi>10.1063/1.2931945]], and its efficient extension (ADMM) described in [[doi>10.1021/ct1002225]].
  
-Hartree-Fock exchange in CP2K is based on four center electron repulsion integrals (ERI), these are computed with an external library ([[http://sourceforge.net/projects/libint/|libint]]). To do these exercises, CP2K must be linked to this library.+Hartree-Fock exchange in CP2K is based on four center electron repulsion integrals (ERI), these are computed with an external library ([[http://sourceforge.net/projects/libint/|libint]]).
  
-This approach has a computational cost that depends strongly on the nature of the basis, unless combined with ADMM (see below), do not use MOLOPT basis sets with HFX. We use basis sets from ''HFX_BASIS'', which are suitable.+This approach has a computational cost that depends strongly on the nature of the basis. Unless combined with ADMM (see below), do not use MOLOPT basis sets with HFX. We use basis sets from ''HFX_BASIS'', which are suitable.
  
 This tutorial is also based on the ''GPW'' method and the calculations will be run using the GTH pseudopotentials. This tutorial is also based on the ''GPW'' method and the calculations will be run using the GTH pseudopotentials.
  
 <note> <note>
-Tutorial re-adapted from [[https://www.cp2k.org/exercises:2016_summer_school:hfx]] +Tutorial re-adapted from the [[https://www.cp2k.org/exercises:2016_summer_school:hfx| Tutorial on HFX and ADMM 2016 at KCL]] 
-For more info see also [[https://www.cecam.org/upload/talk/presentation_5766.pdf]]+For more info see also the slides from Joost VandeVondele [[pdf>https://www.cecam.org/upload/talk/presentation_5766.pdf| presentation 2011]], Sanliang Ling [[pdf>https://www.cp2k.org/_media/events:2016_user_meeting:cp2k-uk-2016-ling.pdf | presentation 2016]] and 
 +Matt Watkins [[https://mattatlincoln.github.io/cp2k_workshop_2017 | presentation 2017]]
 </note> </note>
  
 ===== Truncated Coulomb operator ===== ===== Truncated Coulomb operator =====
  
-To enable HFX in the condensed phase (described at the Gamma point only), CP2K employs a truncated Coulomb operator for the exchange part. The physical picture is that we do not want to have 'self-exchange interactions' of an electron with its image in neighboring unit cells. As a rule of thumb, the maximum range (truncation radius) is L/2 where L is the smallest edge of the unit cell. The convergence of the exchange energy is exponential wrt. this radius. Typically, 5-6A provides good results, but this depends on the nature of the system, i.e. the band gap or the range of the maximally localized Wannier orbitals.+To enable HFX in the condensed phase CP2K employs a truncated Coulomb operator for the exchange part. The physical picture is that we do not want to have 'self-exchange interactions' of an electron with its image in neighboring unit cells. As a rule of thumb, the maximum range (truncation radius) is L/2 where L is the smallest edge of the unit cell. The convergence of the exchange energy is exponential wrt. this radius. Typically, 5-6A provides good results, but this depends on the nature of the system, i.e. the band gap or the range of the maximally localized Wannier orbitals.
  
 ==== 1st task : GGA restart wfn ==== ==== 1st task : GGA restart wfn ====
Line 30: Line 31:
    ! various runtypes (energy, geo_opt, etc.) available.    ! various runtypes (energy, geo_opt, etc.) available.
    RUN_TYPE ENERGY    RUN_TYPE ENERGY
-   ! limit the runs to 5min +   ! amount of information printed to output
-   WALLTIME 1800 +
-   ! reduce the amount of IO+
    IOLEVEL  MEDIUM    IOLEVEL  MEDIUM
 &END GLOBAL &END GLOBAL
Line 73: Line 72:
          &E_DENSITY_CUBE OFF          &E_DENSITY_CUBE OFF
          &END E_DENSITY_CUBE          &END E_DENSITY_CUBE
-         ! compute eigenvalues and homo-lumo gap each 10nd MD step+         ! compute eigenvalues and homo-lumo gap 
          &MO_CUBES          &MO_CUBES
             NLUMO 4             NLUMO 4
Line 87: Line 86:
          SCF_GUESS ATOMIC ! can be used to RESTART an interrupted calculation          SCF_GUESS ATOMIC ! can be used to RESTART an interrupted calculation
          MAX_SCF 30          MAX_SCF 30
-         EPS_SCF 1.0E-6 ! accuracy of the SCF procedure typically 1.0E-6 - 1.0E-7+         EPS_SCF 1.0E-6 ! accuracy of the SCF procedure, for OT typically 1.0E-6 - 1.0E-7, for diagonalization may have to be smaller
          &OT          &OT
             ! an accurate preconditioner suitable also for larger systems             ! an accurate preconditioner suitable also for larger systems
Line 152: Line 151:
    &END SUBSYS    &END SUBSYS
 &END FORCE_EVAL &END FORCE_EVAL
- 
- 
  
 </code> </code>
  
 Run the input and rename the generated wfn file (''WATER-RESTART.wfn'') to ''WATER-RESTART-GGA.wfn''. Run the input and rename the generated wfn file (''WATER-RESTART.wfn'') to ''WATER-RESTART-GGA.wfn''.
-Also make a note of the ''HOMO - LUMO gap [eV]''+Also make a note of the ''HOMO - LUMO gap [eV]'' (the experimental water band-gap is of about 6.9 eV).
  
 ==== 2nd task: PBE0-D3 water ==== ==== 2nd task: PBE0-D3 water ====
Line 168: Line 165:
   * ''WFN_RESTART_FILE_NAME WATER-RESTART-GGA.wfn'' in the ''&DFT'' section   * ''WFN_RESTART_FILE_NAME WATER-RESTART-GGA.wfn'' in the ''&DFT'' section
  
-And replace the previous ''&XC...&END XC'' section with the following one:+And replace the ''&XC...'' section with the following one:
  
 <code> <code>
Line 185: Line 182:
         ! 25 % HFX exchange         ! 25 % HFX exchange
         FRACTION 0.25         FRACTION 0.25
-        &SCREENING +        ! Important to improve scaling from O(N^4) to O(N) 
-          ! important parameter to get stable HFX calcs+        &SCREENING         
 +          ! important parameter to get stable HFX calcs (contributions to hfx smaller than EPS_SCHWARZ are not considered)
           EPS_SCHWARZ 1.0E-6           EPS_SCHWARZ 1.0E-6
           ! needs a good (GGA) initial guess            ! needs a good (GGA) initial guess 
 +          ! screening on the product between maximum of density matrix elements and ERI
           SCREEN_ON_INITIAL_P TRUE           SCREEN_ON_INITIAL_P TRUE
         &END         &END
Line 194: Line 193:
           ! for condensed phase systems           ! for condensed phase systems
           POTENTIAL_TYPE TRUNCATED           POTENTIAL_TYPE TRUNCATED
-          ! should be less than halve the cell +          ! should be less than half the cell 
-          CUTOFF_RADIUS 6.0+          CUTOFF_RADIUS  6.0
           ! data file needed with the truncated operator           ! data file needed with the truncated operator
           T_C_G_DATA ./t_c_g.dat           T_C_G_DATA ./t_c_g.dat
Line 219: Line 218:
  
 Topics: Topics:
-  * ''EPS_SCHWARZ'', ''EPS_PGF_ORB'', ''SCREEN_ON_INITIAL_P'' : parameters to guarantee stable SCF.+  * ''EPS_PGF_ORB'' controls the sparse pattern of the overlap matrix, any contribution of the density matrix smaller than ''EPS_PGF_ORB'' is treated as zero; 
 +  * ''EPS_SCHWARZ'' parameter for Schwarz screening, where contributions to HFX smaller than ''EPS_SCHWARZ'' treated as zero. Typical values range between $10^{-6} and 10^{-9}$; 
 +  * ''SCREEN_ON_INITIAL_P'' uses the density matrix initially available to SCF and starts to use screening on that; 
 +  * ''EPS_PGF_ORB'', ''EPS_SCHWARZ'', ''SCREEN_ON_INITIAL_P'' and ''EPS_FILTER_MATRIX'' are parameters to guarantee stable SCF.
   * Fraction of exchange (''SCALE_X'', ''FRACTION'').    * Fraction of exchange (''SCALE_X'', ''FRACTION''). 
  
-Run this input u 
 Have a look at the output, in the section where CP2K is performing the SCF loop, using the OT method. Have a look at the output, in the section where CP2K is performing the SCF loop, using the OT method.
 The first iteration should look something like this: The first iteration should look something like this:
Line 267: Line 268:
 You can see from the output the appearance of a ''WARNING'' that is related to the You can see from the output the appearance of a ''WARNING'' that is related to the
 tags ''EPS_PGF_ORB'' and ''EPS_FILTER_MATRIX''. It is recommended to check and understand every warning in tags ''EPS_PGF_ORB'' and ''EPS_FILTER_MATRIX''. It is recommended to check and understand every warning in
-the output and if it is OK to discard it, or not. This particular warning is related to the stability of the SCF cycle. For more information see [[https://www.cp2k.org/faq:hfx_eps_warning]].+the output and if it is OK to discard it or not. This particular warning is related to the stability of the SCF cycle. For more information see [[https://www.cp2k.org/faq:hfx_eps_warning]].
  
 +**Questions:**
   * Look at the output where the HOMO-LUMO gap has been printed out. How does this compare to the GGA result?    * Look at the output where the HOMO-LUMO gap has been printed out. How does this compare to the GGA result? 
-  * Most of the time in the SCF cycle is spent in the first step. Why do you think is that? 
   * Adjust the fraction of exchange (modify the input in two places!) to 20% and/or 30%, how does this influence the gap ?   * Adjust the fraction of exchange (modify the input in two places!) to 20% and/or 30%, how does this influence the gap ?
-  * Optional task n. 1. You can check if the SCF cycle is stable by changing the values of ''EPS_PGF_ORB'' and ''EPS_FILTER_MATRIX'', re-run the calculation and see how this affects the ''ENERGY''. +  * The most expensive part of the whole SCF cycle is represented by the first step, while the other steps are much faster. Why is that? 
-  * Optional task n. 2. The most expensive part in hybrid calculations is the calculation of the electron repulsion integrals (ERI). To run efficiently CP2K tries to store the ERI in-core and avoid to calculate them at each SCF step. Especially for large systems that can be run on large HCP machines it is important to maximize the in-core operation and fit the calculations of the ERI into memory. To see the detrimental effect of not having enough memory on the time required for the calculation of exact exchange, modify the tag of ''MAX_MEMORY'' to 40 and rerun the calculation. How do the timings compare with the previous calculations?+  * __Optional__ You can check if the SCF cycle is stable by decreasing the values of ''EPS_PGF_ORB'', ''EPS_FILTER_MATRIX'' as well as ''EPS_SCHWARZ''. Re-run the calculation and see how this affects the ''ENERGY''. 
 +  * __Optional__ CP2K tries to store the ERI in-core and avoid to calculate them at each SCF step. Especially for large systems that can be run on large HCP machines it is important to run in-core operation and fit the calculations of the ERI into memory. To see the effect of not having enough memory on the time for the ''SCF'' cycle, modify the tag of ''MAX_MEMORY'' to 40 and rerun the calculation. How do the timings compare with those where ''MAX_MEMORY'' was larger? 
 + 
 +<note> 
 +It is recommended to run HFX calculations in-core, either by setting larger ''MAX_MEMORY'' or using more MPI processes. 
 +</note>
  
 ===== Truncated Coulomb operator with long range correction ===== ===== Truncated Coulomb operator with long range correction =====
  
-Like in the HSE functional, the difference between the operator used for exchange and 1/r, can be accounted for by a special GGA exchange functional. Also for the truncated coulomb operator this is possible, and allows for xc functionals that embed very short range exchange operators only. This can be used to speedup the calculation, while retaining the benefits of HFX. The functional employed in this way smoothly goes from PBE to PBE0 as the range goes from 0 to Infinity.+In HSE and other screened hybrid functionals, the 1/r exchange interaction potential can be separated into a long range and a short range part, where only the short range part is computed using exact exchange and the long-range part is computed using GGA exchange. This is what is done also when using a truncated Coulomb operator. This can be used to speedup the calculation, while retaining the benefits of HFX. The functional employed in this way smoothly goes from PBE to PBE0 as the ''CUTOFF_RADIUS'' goes from 0 to Infinity.
  
 ==== 3rd task ==== ==== 3rd task ====
Line 294: Line 300:
 Rerun the single point energy calculation and note the band gap.  Rerun the single point energy calculation and note the band gap. 
   * Is such a short range sufficient to have a sizable effect on the band gap ?   * Is such a short range sufficient to have a sizable effect on the band gap ?
-  * is ''HFX_MEM_INFO| Number of cart. primitive ERI's calculated'' very different for calculations with 2.5 and 6.0A truncation radius ? +  * is ''HFX_MEM_INFO| Number of cart. primitive ERI's calculated'' very different for calculations with 2.5 and 6.0A truncation radius? Remember to use the same cutoff radius under the ''&PBE_HOLE_T_C_LR'' and ''&INTERACTION_POTENTIAL'' sections.
  
 ===== Auxiliary Density Matrix Methods (ADMM) ===== ===== Auxiliary Density Matrix Methods (ADMM) =====
Line 303: Line 308:
 ==== 4rd task : introduce ADMM ==== ==== 4rd task : introduce ADMM ====
  
-Make the following changes: +Make the following changes to the input: 
-  * insert and additional line ''BASIS_SET_FILE_NAME BASIS_ADMM'' (and copy that file from cp2k/data as needed).+  * insert an additional line ''BASIS_SET_FILE_NAME BASIS_ADMM''.
   * insert for each ''&KIND'' a line ''AUX_FIT_BASIS_SET cFIT3''   * insert for each ''&KIND'' a line ''AUX_FIT_BASIS_SET cFIT3''
   * insert a secion ''&AUXILIARY_DENSITY_MATRIX_METHOD''   * insert a secion ''&AUXILIARY_DENSITY_MATRIX_METHOD''
Line 321: Line 326:
  
 <note> In this tutorial we combine ADMM using a very small basis set (cFIT3) with a small primary basis (DZVP-GTH), so gains are small at best (and the results not very accurate). ADMM is most useful with good quality primary basis sets, such as e.g. MOLOPTs</note> <note> In this tutorial we combine ADMM using a very small basis set (cFIT3) with a small primary basis (DZVP-GTH), so gains are small at best (and the results not very accurate). ADMM is most useful with good quality primary basis sets, such as e.g. MOLOPTs</note>
- 
-<note> For CP2K 2.7, a [[https://sourceforge.net/p/cp2k/code/HEAD/tree/trunk/cp2k/data/BASIS_ADMM_MOLOPT|new library of ADMM basis sets]] is available.</note> 
  
 Run the input, what's the ''HOMO - LUMO gap'' ? Run the input, what's the ''HOMO - LUMO gap'' ?
Line 336: Line 339:
 adapt the admm input for water to reflect the ionized state: adapt the admm input for water to reflect the ionized state:
 <code> <code>
-    ! Charge and multiplicity+    ! Spin polarization, charge and multiplicity
     LSD     LSD
     CHARGE 1     CHARGE 1
Line 360: Line 363:
   RUN_TYPE ENERGY   RUN_TYPE ENERGY
   ! limit the runs to 30min   ! limit the runs to 30min
-  WALLTIME 1800+  WALLTIME 1800 
   ! reduce the amount of IO   ! reduce the amount of IO
   IOLEVEL  MEDIUM    IOLEVEL  MEDIUM 
Line 376: Line 379:
  
     ! Charge and multiplicity     ! Charge and multiplicity
-    CHARGE 0 +    CHARGE 1 
-    MULTIPLICITY 1+    MULTIPLICITY 2
  
     &MGRID     &MGRID
exercises/2017_uzh_cp2k-tutorial/hybrid.1499800605.txt.gz · Last modified: 2020/08/21 10:15 (external edit)