CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace: • docs
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 [2016/08/23 20:32] – mwatkinsevents:2016_summer_school:gpw [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 4: Line 4:
 # heavily lifted from Marcella Iannuzzi's previous presentations, e.g # heavily lifted from Marcella Iannuzzi's previous presentations, e.g
 # https://www.cp2k.org/_media/events:2015_cecam_tutorial:iannuzzi_gpw_gapw_b.pdf # https://www.cp2k.org/_media/events:2015_cecam_tutorial:iannuzzi_gpw_gapw_b.pdf
 +#
 +# Matt Watkins
 +# School of Mathematics and Physics, University of Lincoln, UK
 +
 </code> </code>
 ===== GPW (GAPW) electronic structure calculations ===== ===== GPW (GAPW) electronic structure calculations =====
-Matt Watkins 
-School of Mathematics and Physics, University of Lincoln, UK 
  
 http://www.cp2k.org http://www.cp2k.org
Line 21: Line 23:
   * many observables available (spectroscopy)   * many observables available (spectroscopy)
   * cheaper (better scaling) than many wavefunction based approaches   * cheaper (better scaling) than many wavefunction based approaches
-  * *ab initio* molecular dynamics / Monte Carlo+  * //ab initio// molecular dynamics / Monte Carlo
  
  
-=== Hohenberg-Kohn Theorems ===+==== Hohenberg-Kohn Theorems ====
  
 **Theorem 1** **Theorem 1**
Line 45: Line 47:
   * $E_{xc}$ - non-classical Coulomb energy: electron correlation   * $E_{xc}$ - non-classical Coulomb energy: electron correlation
  
-{{https://upload.wikimedia.org/wikipedia/commons/d/dc/Kohn.jpg}} +==== Kohn-Sham: non-interacting electrons ====
-[[https://en.wikipedia.org/wiki/Walter_Kohn|Walter Kohn]] +
-{{http://physics.nyu.edu/~pch2/pierre.hohenberg.png}} +
-Pierre Hohenberg +
- +
-=== Kohn-Sham: non-interacting electrons ===+
  
 We map the problem of interacting electrons onto an equivalent one with non-interacting electrons with an altered external potential (xc). We map the problem of interacting electrons onto an equivalent one with non-interacting electrons with an altered external potential (xc).
  
-We expand the *density* as a sum of one electron orbital contributions+We expand the //density// as a sum of one electron orbital contributions
  
 === Electronic density === === Electronic density ===
Line 82: 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}}+==== KS equations ====
  
-[[https://en.wikipedia.org/wiki/Lu_Jeu_Sham|Lu Jeu Sham]]+Using the 2nd Hohenberg-Kohn theorem, we find the electron density that minimises the energy of the system. 
  
-=== KS equations ===+But, like in Hartree-Fock theory, we have to ensure that the electron orbitals are orthonormal to prevent the system imploding.
  
-Using the 2nd Hohenberg-Kohn theorem, we find the electron density that minimises the energy of the system. 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 96: 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 107: Line 102:
 $$ $$
  
-**Kohn-Sham functional**+=== Kohn-Sham functional ===
  
 $$ $$
Line 117: Line 112:
 $$ $$
  
-The exact functional form for the electron-electron repulsion is not known, but various levels of approximation are available (Jacob's La +The exact functional form for the electron-electron repulsion is not known, but various levels of approximation are available (Jacob's Ladder).  
-dder). 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.
  
-=== KS equations ===+==== KS equations ====
  
 $$ $$
Line 138: Line 134:
 - KS scheme is, in principle, exact (if we knew the exact $E_{XC}$) - KS scheme is, in principle, exact (if we knew the exact $E_{XC}$)
  
-=== 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).
  
-=== Basis Set Representation(s) ===+==== Basis Set Representation(s) ====
  
-to actually do calculations we need to expand the KS equations into a basis. In CP2K we use Gaussian functions, as the primary basis set:+to actually do calculations we need to expand the KS equations into a basis. In CP2K we use Gaussian functions as the primary basis set:
  
 the KS orbitals are then represented as the KS orbitals are then represented as
Line 201: Line 197:
 Where $E[\{\psi_i\};\{\mathbf{R}_I\}]$ indicates that the energy depends on the (occupied) KS orbitals (which then give a density) and parametrically on all the $I$ nuclear coordinates, $\mathbf{R}_I$. Where $E[\{\psi_i\};\{\mathbf{R}_I\}]$ indicates that the energy depends on the (occupied) KS orbitals (which then give a density) and parametrically on all the $I$ nuclear coordinates, $\mathbf{R}_I$.
  
-**Matrix formulation of the KS equations**+==== Matrix formulation of the KS equations ====
  
 using the variational principle we get using the variational principle we get
Line 209: Line 205:
 $$ $$
  
-=== Algorithmic challenges ===+==== Algorithmic challenges ====
  
   * construct the Kohn-Sham matrix   * construct the Kohn-Sham matrix
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.1471984330.txt.gz · Last modified: 2020/08/21 10:14 (external edit)