CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
events:2016_summer_school:gpw

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
events:2016_summer_school:gpw [2018/05/29 21:32] – [Hohenberg-Kohn Theorems] mwatkinsevents:2016_summer_school:gpw [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 79: Line 79:
 \Psi = \frac{1}{\sqrt{N!}} \text{det} [\psi_1 \psi_2 \psi_3 \cdots \psi_N] \Psi = \frac{1}{\sqrt{N!}} \text{det} [\psi_1 \psi_2 \psi_3 \cdots \psi_N]
 $$ $$
-{{https://physics.ucsd.edu/images/people/lsham/lsham1_144x216.jpg}} 
- 
-[[https://en.wikipedia.org/wiki/Lu_Jeu_Sham|Lu Jeu Sham]] 
- 
 ==== KS equations ==== ==== KS equations ====
  
Line 89: Line 85:
 But, like in Hartree-Fock theory, we have to ensure that the electron orbitals are orthonormal to prevent the system imploding. But, like in Hartree-Fock theory, we have to ensure that the electron orbitals are orthonormal to prevent the system imploding.
  
-**Orthogonality constraint**+=== Orthogonality constraint ===
  
 $$ $$
Line 95: Line 91:
 $$ $$
  
-**Variational search in the space of the orbitals**+=== Variational search in the space of the orbitals ===
  
 We correct the non-interacting electron model by adding in an (in principle unknown) XC potential that accounts for **all** quantum mechanical many-body interactions (electron-electron repulsion) We correct the non-interacting electron model by adding in an (in principle unknown) XC potential that accounts for **all** quantum mechanical many-body interactions (electron-electron repulsion)
  
-**Classical election-electron repulsion**+=== Classical election-electron repulsion ===
  
 $$ $$
Line 106: Line 102:
 $$ $$
  
-**Kohn-Sham functional**+=== Kohn-Sham functional ===
  
 $$ $$
Line 118: Line 114:
 The exact functional form for the electron-electron repulsion is not known, but various levels of approximation are available (Jacob's Ladder).  The exact functional form for the electron-electron repulsion is not known, but various levels of approximation are available (Jacob's Ladder). 
  
-The existence of this functional is guarenteed by the 1st Hohenberg-Kohn Theroem.+The existence of this functional is guaranteed by the 1st Hohenberg-Kohn Theroem.
  
 This maps mathematically onto the familiar Hartree-Fock model of electronic structure. This maps mathematically onto the familiar Hartree-Fock model of electronic structure.
Line 140: Line 136:
 ==== Self-consistency ==== ==== Self-consistency ====
  
-''generate an initial guess, from a superposition of atomic densities (typical PW code) or atomic block diagonal density matrix (CP2K)''+//generate an initial guess, from a superposition of atomic densities (typical PW code) or atomic block diagonal density matrix (CP2K)//
  
-  - generate a starting density $\Rightarrow n^{init}$ +  - generate a starting density $\Rightarrow n^{0}$ 
-  - generate the KS potential   $\Rightarrow V_{KS}^{init}$+  - generate the KS potential   $\Rightarrow V_{KS}^{0}$
   - solve the KS equations  $\Rightarrow \epsilon, \psi$   - solve the KS equations  $\Rightarrow \epsilon, \psi$
  
- +//then repeat until convergence//
-''then repeat until convergence''+
  
   - calculate the new density $\Rightarrow n^{1}$    - calculate the new density $\Rightarrow n^{1}$ 
Line 154: Line 149:
  
  
-''calculate properties from final density and orbitals''+//calculate properties from final density and orbitals//
  
 Stopping criteria in CP2K are: Stopping criteria in CP2K are:
 +
   * the largest MO derivative ($\frac{\partial E}{\partial C_{ \alpha i}} $) is smaller than `EPS_SCF` (OT)   * the largest MO derivative ($\frac{\partial E}{\partial C_{ \alpha i}} $) is smaller than `EPS_SCF` (OT)
   * the largest change in the density matrix ($P$) is smaller than `EPS_SCF` (traditional diagonalization).   * the largest change in the density matrix ($P$) is smaller than `EPS_SCF` (traditional diagonalization).
Line 244: Line 240:
 This leads to linear scaling KS construction for Gaussian Type Orbitals (GTO) This leads to linear scaling KS construction for Gaussian Type Orbitals (GTO)
  
-  * Guassian basis sets (many matrix elements can be done analytically)  +==== Guassian basis sets (many matrix elements can be done analytically) ==== 
 we go a bit further than implied above - to be more accurate, we *contract* several Gaussians to form approximate atomic orbitals we go a bit further than implied above - to be more accurate, we *contract* several Gaussians to form approximate atomic orbitals
 $$ $$
Line 254: Line 251:
 $$ $$
 and $m_x + m_y + m_z = l$, the angular momentum quantum number of the functions. and $m_x + m_y + m_z = l$, the angular momentum quantum number of the functions.
 +
   * Pseudo potentials   * Pseudo potentials
   * Plane waves auxiliary basis for Coulomb integrals   * Plane waves auxiliary basis for Coulomb integrals
Line 279: Line 277:
 $$ $$
  
-=== Data files ===+==== Data files ====
  
 Parameter files distributed with CP2K in $CP2K_ROOT/cp2k/data Parameter files distributed with CP2K in $CP2K_ROOT/cp2k/data
Line 298: Line 296:
   * non-local dispersion functionals `vdW_kernel_table.dat, rVV10_kernel_table.dat`   * non-local dispersion functionals `vdW_kernel_table.dat, rVV10_kernel_table.dat`
  
-=== Basis set libraries ===+==== Basis set libraries ====
  
 There are two main types of basis sets supplied with CP2K There are two main types of basis sets supplied with CP2K
Line 327: Line 325:
 </code> </code>
  
-Here there are gaussians with five different exponents ( in Bohr$^{-2}$). From these 5 sets of functions are built, two $s$ functions (2nd and 3rd columns), 2 sets of $p$ functions (4th and fifth columns) and one set of $d$  functions (last column).+Here there are Gaussians with five different exponents ( in Bohr$^{-2}$). From these 5 sets of functions are built, two $s$ functions (2nd and 3rd columns), 2 sets of $p$ functions (4th and fifth columns) and one set of $d$  functions (last column).
  
 Note  that the contraction coefficients are not varied during calculation. For the nitrogen basis above we have $2 + 2 \times 3 + 1 \times 5 = 13$ variables to optimize for each nitrogen atom in the system. Note  that the contraction coefficients are not varied during calculation. For the nitrogen basis above we have $2 + 2 \times 3 + 1 \times 5 = 13$ variables to optimize for each nitrogen atom in the system.
Line 334: Line 332:
  
  
-=== GTH pseudopotentials ===+==== GTH pseudopotentials ====
  
 Accurate and transferable with few parameters.   Accurate and transferable with few parameters.  
Line 342: Line 340:
   * Norm-conserving, separable, dual-space   * Norm-conserving, separable, dual-space
   * Local PP : short-range and long-range terms   * Local PP : short-range and long-range terms
 +
 $$ $$
 V_{loc}^{PP}(r) = \sum_{i=1}^4 C_i^{PP} (\sqrt{2} \alpha^{PP} r)^{2i-2}e^{-(\alpha^{PP}r)^2} - \frac{Z_{ion}}{r}erf(\alpha^{PP}r) V_{loc}^{PP}(r) = \sum_{i=1}^4 C_i^{PP} (\sqrt{2} \alpha^{PP} r)^{2i-2}e^{-(\alpha^{PP}r)^2} - \frac{Z_{ion}}{r}erf(\alpha^{PP}r)
 $$ $$
 +
 first term in sum is short ranged, and analytically calculated in Gaussian basis. Second term is long ranged, and merged into the electrostatic calculation (see later) first term in sum is short ranged, and analytically calculated in Gaussian basis. Second term is long ranged, and merged into the electrostatic calculation (see later)
  
   * non-local PP with Gaussian type operators   * non-local PP with Gaussian type operators
 +
 $$ $$
 V_{nl}^{PP}(\mathbf{r},\mathbf{r'}) = \sum_{lm} \sum_{ij} \big{<} \mathbf{r} \mid p_i^{lm} \big{>}h^l_{ij} V_{nl}^{PP}(\mathbf{r},\mathbf{r'}) = \sum_{lm} \sum_{ij} \big{<} \mathbf{r} \mid p_i^{lm} \big{>}h^l_{ij}
 \big{<} p_j^{lm} \mid \mathbf{r'} \big{>}  \big{<} p_j^{lm} \mid \mathbf{r'} \big{>} 
 $$ $$
 +
 $$ $$
 \big{<} \mathbf{r} \mid p_i^{lm} \big{>} = N_i^l Y^{lm}(\hat{r})r^{l+2i-2}e^{-\frac{1}{2}(\frac{r}{r_l})^2} \big{<} \mathbf{r} \mid p_i^{lm} \big{>} = N_i^l Y^{lm}(\hat{r})r^{l+2i-2}e^{-\frac{1}{2}(\frac{r}{r_l})^2}
 $$ $$
  
-You can scan through potentials available at [[http://cp2k.web.psi.ch/potentials/|Matthias Krack's website]]+You can scan through potentials available at [[https://www.cp2k.org/static/potentials/|here]]
  
 Original papers:   Original papers:  
-[[Goedeker, Teter, Hutter, PRB 54 (1996), 1703]http://journals.aps.org/prb/abstract/10.1103/PhysRevB.54.1703]   +[[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.54.1703| Goedeker, Teter, Hutter, PRB 54 (1996), 1703]] and   
-[[Hartwigsen, Goedeker, Hutter, PRB 58 (1998) 3641]http://journals.aps.org/prb/abstract/10.1103/PhysRevB.58.3641]+[[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.58.3641| Hartwigsen, Goedeker, Hutter, PRB 58 (1998) 3641]]
  
-=== Electrostatics ===+==== Electrostatics ====
  
   * long-range term: Non-local Hartree Potential     * long-range term: Non-local Hartree Potential  
 +
 Poisson equation solved using the auxiliary plane-wave basis Poisson equation solved using the auxiliary plane-wave basis
 $$ $$
 E_H[n_{tot}] = \frac{1}{2} \int_r \text{d}\mathbf{r} \int_{r'} \text{d}\mathbf{r'} \frac{n(\mathbf{r})n(\mathbf{r'})}{\mid \mathbf{r} - \mathbf{r'}} E_H[n_{tot}] = \frac{1}{2} \int_r \text{d}\mathbf{r} \int_{r'} \text{d}\mathbf{r'} \frac{n(\mathbf{r})n(\mathbf{r'})}{\mid \mathbf{r} - \mathbf{r'}}
 $$ $$
 +
 where $n_{tot}$ includes the nuclear charge as well as the electronic.   where $n_{tot}$ includes the nuclear charge as well as the electronic.  
 (The nuclear charge density is (of course) represented as a Gaussian distribution with parameter $R_I^c$ chosen to cancel a similar term from the local part of the pseudopotential) (The nuclear charge density is (of course) represented as a Gaussian distribution with parameter $R_I^c$ chosen to cancel a similar term from the local part of the pseudopotential)
Line 374: Line 378:
   * FFT (scaling as $N\text{log}N$) gives    * FFT (scaling as $N\text{log}N$) gives 
 $$\tilde{n}(\mathbf{r}) = \frac{1}{\Omega} \sum_G \tilde{n}(\mathbf{G})e^{i\mathbf{G \cdot r}})$$ $$\tilde{n}(\mathbf{r}) = \frac{1}{\Omega} \sum_G \tilde{n}(\mathbf{G})e^{i\mathbf{G \cdot r}})$$
 +
 In the $G$ space representation the Poisson equation is diagonal and the Hartree energy is easily evaluated In the $G$ space representation the Poisson equation is diagonal and the Hartree energy is easily evaluated
 +
 $$ $$
 E_H[n_{tot}] = 2 \pi \Omega \sum_G  \frac{\tilde{n}(\mathbf{G}) ^* \tilde{n}(\mathbf{G}) }{\mathbf{G}^2} E_H[n_{tot}] = 2 \pi \Omega \sum_G  \frac{\tilde{n}(\mathbf{G}) ^* \tilde{n}(\mathbf{G}) }{\mathbf{G}^2}
Line 382: Line 388:
 {{http://ars.els-cdn.com/content/image/1-s2.0-S0010465505000615-gr001.gif}} {{http://ars.els-cdn.com/content/image/1-s2.0-S0010465505000615-gr001.gif}}
  
-=== Real space integration ===+==== Real space integration ====
  
 Finite cutoff and simulation box define a realspace grid Finite cutoff and simulation box define a realspace grid
Line 392: Line 398:
 \sum_{\alpha \beta} P_{\alpha \beta} \bar{\phi}_{\alpha \beta} (\mathbf{r}) = n(\mathbf{R}) \sum_{\alpha \beta} P_{\alpha \beta} \bar{\phi}_{\alpha \beta} (\mathbf{r}) = n(\mathbf{R})
 $$ $$
-where $n(\mathbf{R})$ is the density at grid points in the cell, and $\bar{\phi}_{\alpha \beta}$ are the products of two basis functions +where $n(\mathbf{R})$ is the density at grid points in the cell, and $\bar{\phi}_{\alpha \beta}$ are the products of two basis functions.
-{{materials/collocate.bmp}}+
  
   * numerical approximation of the gradient $n(\mathbf{R}) \rightarrow \nabla n(\mathbf{R})$   * numerical approximation of the gradient $n(\mathbf{R}) \rightarrow \nabla n(\mathbf{R})$
Line 420: Line 425:
 We see that 11.9999999997 electrons are mapped onto the grids, along with 11.9999999994 nuclear charge. You should aim for at a very minimum accuracy of 10−8. If not, increase the cutoff. We see that 11.9999999997 electrons are mapped onto the grids, along with 11.9999999994 nuclear charge. You should aim for at a very minimum accuracy of 10−8. If not, increase the cutoff.
  
-=== Energy ripples ===+==== Energy ripples ====
  
 Low density regions can cause unphysical behaviour of $XC$ terms (such as $\frac{\mid \nabla n \mid ^2}{n^{\alpha}}$) Low density regions can cause unphysical behaviour of $XC$ terms (such as $\frac{\mid \nabla n \mid ^2}{n^{\alpha}}$)
Line 426: Line 431:
 {{http://ars.els-cdn.com/content/image/1-s2.0-S0010465505000615-gr002.gif}} {{http://ars.els-cdn.com/content/image/1-s2.0-S0010465505000615-gr002.gif}}
  
-GTH pseudos have small density at the core - graph of density and $v_{XC}$ through a water molecule. These spikes can cause ripples in the energy as atoms move relative to the grid. +GTH pseudos have small density at the core - graph of density and $v_{XC}$ through a water molecule. These spikes can cause ripples in the energy as atoms move relative to the grid. These can be very problematic when trying to calculate vibrational frequencies.
  
-There are smoothing routines `&XC_GRID / XC_DERIV`, but probably best to stick with the defaults ... whatever you do don't change settings between simulations you want to compare.+ 
 +There are smoothing routines `&XC_GRID / XC_DERIV`, see the exercise [[events:2018_summer_school:converging_cutoff]].
  
 {{http://ars.els-cdn.com/content/image/1-s2.0-S0010465505000615-gr003.gif}} {{http://ars.els-cdn.com/content/image/1-s2.0-S0010465505000615-gr003.gif}}
  
-avoid with higher cutoff, or GAPW methodology.+Avoid ripples with higher a cutoff, or GAPW methodology. 
  
-These can be very problematic when trying to calculate vibrational frequencies.+Whatever you do don't change settings between simulations you want to compare.
  
-=== Multigrids ===+ 
 +==== Multigrids ====
  
 When we want to put (collocate) a Gaussian type function onto the realspace grid, we can gain efficiency by using multiple grids with differing cutoff / spacing. When we want to put (collocate) a Gaussian type function onto the realspace grid, we can gain efficiency by using multiple grids with differing cutoff / spacing.
Line 446: Line 453:
 specified in input as: specified in input as:
  
 +<code>
   &MGRID   &MGRID
     CUTOFF 400     CUTOFF 400
Line 452: Line 459:
     NGRIDS 5     NGRIDS 5
   &END MGRID   &END MGRID
 +</code>
  
 you can see in the output you can see in the output
Line 469: Line 477:
 For this system (formaldeyhde with an aug-TZV2P-GTH basis) that 45436 density matrix elements ($\bar{\phi}_{\alpha \beta}$) were mapped onto the grids. To be efficient, all grids should be used. For this system (formaldeyhde with an aug-TZV2P-GTH basis) that 45436 density matrix elements ($\bar{\phi}_{\alpha \beta}$) were mapped onto the grids. To be efficient, all grids should be used.
  
-To fully converge calculations both ''CUTOFF'' and ''REL_CUTOFF'' need to be increased together.+To fully converge calculations both ''CUTOFF'' and ''REL_CUTOFF'' need to be increased together otherwise the highest grid may not be used. 
 + 
 +==== Timings ====
  
 At the end of the run you'll see timings - these can be very useful for understanding performance. At the end of the run you'll see timings - these can be very useful for understanding performance.
Line 521: Line 531:
  -------------------------------------------------------------------------------  -------------------------------------------------------------------------------
 </code> </code>
 +
 +==== Warnings ====
  
 Also check if you get an output like Also check if you get an output like
events/2016_summer_school/gpw.1527629564.txt.gz · Last modified: 2020/08/21 10:14 (external edit)