Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision		 Previous revision
the_density_functional_program [2023/06/23 15:25] – [Implementing new functionals] hesselmannthe_density_functional_program [2025/01/09 09:54] (current) – fix layout doll
Line 307: Line 307:
 |                           | ''NEESE_INT_ACC''                         | 4.67          | Neese parameter for number of radial points for Neese type grid when ''NEESE_INDEX'' < 3  |  | ✔ | |                           | ''NEESE_INT_ACC''                         | 4.67          | Neese parameter for number of radial points for Neese type grid when ''NEESE_INDEX'' < 3  |  | ✔ |
 | other                     | ''WEIGHT_CUT|WCUT''                    | 1d-20         | threshold value to discard grid points with small integration weights | ✔ | ✔ |  | other                     | ''WEIGHT_CUT|WCUT''                    | 1d-20         | threshold value to discard grid points with small integration weights | ✔ | ✔ | 
 +|                           | ''ORIENT''                             | 2             | orientation of atomic grids | ✔ | |  
  
  
Line 469: Line 470:
  
 Grid points very close to the nucleus can have very small grid weights which can drop below machine precsion (~10d-16) when the radial grid is very large. These can be discarded with the option ''WEIGHT_CUT''=//thr//, i.e., grid points with weights smaller than //thr// will then not be used in the numerical integration anymore. While this can make the integration a little bit more efficient, a truncation of the grid with this option can prevent to achieve high accuracies in the numerical integration below, e.g., micro-hartree accuracies for energies. Therefore by default ''WEIGHT_CUT'' is set to a very small number that is well below machine precision (''WEIGHT_CUT''=1d-20) so that effectively all grid points close to the nucleus are taken into account for the integration.  Grid points very close to the nucleus can have very small grid weights which can drop below machine precsion (~10d-16) when the radial grid is very large. These can be discarded with the option ''WEIGHT_CUT''=//thr//, i.e., grid points with weights smaller than //thr// will then not be used in the numerical integration anymore. While this can make the integration a little bit more efficient, a truncation of the grid with this option can prevent to achieve high accuracies in the numerical integration below, e.g., micro-hartree accuracies for energies. Therefore by default ''WEIGHT_CUT'' is set to a very small number that is well below machine precision (''WEIGHT_CUT''=1d-20) so that effectively all grid points close to the nucleus are taken into account for the integration. 
-==== Grid caching (GRIDSAVE, NOGRIDSAVE) ==== 
- 
-''NOGRIDSAVE'' 
- 
-disables the disk caching of the grid, i.e, forces the recalculation of the grid each time it is needed. 
- 
-''GRIDSAVE'' 
- 
-forces the use of a grid cache where possible. 
  
 ==== Grid symmetry (GRIDSYM,NOGRIDSYM) ==== ==== Grid symmetry (GRIDSYM,NOGRIDSYM) ====
Line 495: Line 487:
 controls printing of the grid, which by default is not done. At present, the only possible value for //key// is ''GRID'', and //value// should be specified as an integer. ''GRID=0'' causes the total number of integration points to be evaluated and reported; ''GRID=1'' additionally shows the number of points on each atom; ''GRID=2'' causes the complete set of grid points and weights to be printed. controls printing of the grid, which by default is not done. At present, the only possible value for //key// is ''GRID'', and //value// should be specified as an integer. ''GRID=0'' causes the total number of integration points to be evaluated and reported; ''GRID=1'' additionally shows the number of points on each atom; ''GRID=2'' causes the complete set of grid points and weights to be printed.
  
 +==== Orientation of atomic grids (ORIENT) ====
 +
 +The orientation of the atomic grids can be controled with the option ''ORIENT'' that can take the values 0, 1 and 2. ''ORIENT''=0 means that the grids are not oriented with respect to the atomic coordinates. This, however, can lead to non-invariant KS energies with respect to rotations of the molecule which is particularly significant for small grids. The setting ''ORIENT''=1 uses the method of B. G. Johnson //et al.// [[https://doi.org/10.1016/0009-2614(94)00199-5|Chem. Phys. Lett.]] **220**, 377 (1994) which constructs rotationally invariant DFT grids. The method of Johnson //et al.// uses, however, the whole molecule to determine the orientation on each atom which might not be ideal for very large molecules or molecular clusters in various applications (determination of conformer energy differences, potential energy surfaces, or geometry optimisations). This can be resolved by the setting ''ORIENT''=2 which emphasises the local environment in the orientation of the atomic grids. This is the current default in Molpro (version 2023.1 and higher).
 ===== Adaptive Molpro Grid (AMG) ==== ===== Adaptive Molpro Grid (AMG) ====
 By default Molpro generates a quadrature grid which prunes the angular grid using an adaptive scheme which takes the molecular environment into account. Several pruning methods and functions can be set. The pruning method used can be controled via the option By default Molpro generates a quadrature grid which prunes the angular grid using an adaptive scheme which takes the molecular environment into account. Several pruning methods and functions can be set. The pruning method used can be controled via the option
Line 970: Line 965:
  
  
 +**Important note regarding some functionals from LibXC that also are available in Molpro's internal functional library:**
 +A number of functionals which are both available in LibXC and Molpro may give slightly different results in total energies (in the micro hartree range for small molecules) due to truncation of significant digits of underlying parameters. An example for this is the PBE functional [[https://doi.org/10.1103/PhysRevLett.77.3865]] in which case the implementation in Molpro uses parameters taken only from the PBE paper some of which are truncated compared to the original source, see references in the PBE paper. Compared to this, LibXC uses a more accurate representation of the parameters in terms of significant digits, resulting to differences in the micro hartree range. Such differences are unlikely to have any effect on energy differences, such as reaction or atomisation energies, but should be kept in mind when comparing total energies calculated with the Molpro and LibXC functionals, respectively.
  
 +Functionals known to yield slighty different results according to this:
  
- +^ Molpro ^ LibXC ^ 
- +| PBEX    | GGA_X_PBE | 
- +| PBEXREV  | GGA_X_PBE_R |  
- +| PBEC     | GGA_C_PBE | 
 +| PBESOLC  | GGA_C_PBE_SOL | 
 +| PW91C    | GGA_C_PW91 | 
 +| HCTH93   | GGA_XC_HCTH_93 | 
 +| HCTH120  | GGA_XC_HCTH_120 | 
 +| M06LX    | MGGA_X_M06_L | 
 +| B95      | MGGA_C_BC95 | 
 +| TPSSC    | MGGA_C_TPSS | 
 +| M06C     | MGGA_C_M06 | 
 +| M06HFC   | MGGA_C_M06_HF | 
 +| M062XC   | MGGA_C_M06_2X | 
 +| M06LC    | MGGA_C_M06_L | 
 +| M05C     | MGGA_C_M05 | 
 +| M052XC   | MGGA_C_M05_2X | 
 +| M06X     | HYB_MGGA_X_M06 | 
 +| M062XX   | HYB_MGGA_X_M06_2X |
  
 ==== Implementing new functionals ==== ==== Implementing new functionals ====
Line 1036: Line 1048:
 It has been observed by Mardirossian and Head-Gordon that several density functionals from the Minnesota group, including M06, M11, M06-L, M06-HF and M11-L, exhibit a very slow convergence of intermolecular interaction energies with respect to the basis set, see [[https://pubs.acs.org/doi/10.1021/ct400660j|JCTC 9 (2013) 4453]]. Apparently, this problem can not be cured by larger integration grids and originates from a contribution to the exchange enhancement factor whose shape can strongly depend on the chosen basis set. Based on these observations it is recommended to use basis sets which were originally employed to develop these functionals. It has been observed by Mardirossian and Head-Gordon that several density functionals from the Minnesota group, including M06, M11, M06-L, M06-HF and M11-L, exhibit a very slow convergence of intermolecular interaction energies with respect to the basis set, see [[https://pubs.acs.org/doi/10.1021/ct400660j|JCTC 9 (2013) 4453]]. Apparently, this problem can not be cured by larger integration grids and originates from a contribution to the exchange enhancement factor whose shape can strongly depend on the chosen basis set. Based on these observations it is recommended to use basis sets which were originally employed to develop these functionals.
  
 +==== Note on MK00 and MK00B functionals === 
 +The MK00 and MK00B functionals (termed MGGA_X_MK00 and MGGA_X_MK00B in LibXC) [[https://doi.org/10.1063/1.481298|J. Chem. Phys. 112, 7002–7007 (2000)]] should better not be used in general calculations with Molpro due to the problem that the denominators τ−υ/4 in the functional can be zero or even turn negative, producing unphysical results.
  
  
Line 1293: Line 1306:
  
 The total energy will then be calculated as EDFT-D2=EDFT+Edisp(D2) if x=d2 (see Ref. [2]), The total energy will then be calculated as EDFT-D2=EDFT+Edisp(D2) if x=d2 (see Ref. [2]),
-EDFT-D3=EDFT+Edisp(D3) if x=d3 (see Ref. [3]), and EDFT-D4=EDFT+Edisp(D4) if x=d4 (see Ref. [4]).+EDFT-D3=EDFT+Edisp(D3) if x=d3 (see Ref. [3]), and 
 +EDFT-D4=EDFT+Edisp(D4) if x=d4 (see Ref. [4]).
  
 Currently the default dispersion correction added to the DFT energy is the D4 dispersion correction developed by Grimme //et al.//, see Ref. [4].  Currently the default dispersion correction added to the DFT energy is the D4 dispersion correction developed by Grimme //et al.//, see Ref. [4]. 
Line 1383: Line 1397:
  
 Gradient contributions from the D3 and D4 dispersion correction are automatically computed in DFT geometry optimisations. Gradient contributions from the D3 and D4 dispersion correction are automatically computed in DFT geometry optimisations.
 +
 +Note that not all functionals implemented in Molpro (and Libxc) are known to the D4 program (and vice versa). And in some cases the functionals supported by D4 might have a different identifier than used in Molpro. To see which functionals can be used with D4, see the documentation at [[https://github.com/dftd4/dftd4]]. Most functionals supported by D4 can also be given by their Libxc names, for example
 +<code> ks,gga_x_b88,gga_c_p86,disp=d4 </code>
 +using the ''B88'' exchange and ''P86'' correlation functional names as defined in the Libxc library. For compound functionals like B+P86 it is important to //insert the exchange functional first and the correlation functional second//!
 +
 +
  
 References:\\ References:\\
Line 1446: Line 1466:
 where ''orbital_index'' is the index of the orbital in the order used for constructing the density matrix, and ''occupation_number'' is the occupation number of this orbital between 0 and 1. These commands can be repeated in case of several fractionally occupied orbitals. When specifying the charge and the wave function, all fractionally occupied orbitals must be considered as occupied orbitals. where ''orbital_index'' is the index of the orbital in the order used for constructing the density matrix, and ''occupation_number'' is the occupation number of this orbital between 0 and 1. These commands can be repeated in case of several fractionally occupied orbitals. When specifying the charge and the wave function, all fractionally occupied orbitals must be considered as occupied orbitals.
  
-For DFT calculations, the implementation has only been done for the option ''matrix=0'' in the ''uks'' command.+For DFT calculations, the implementation has only been done for the option ''matrix=0'' in the ''uks'' command. Furthermore, the option ''old'' must be added to the ''uks'' command in order to redirect Molpro to the old unrestricted SCF programn (the setting of fractional occupations is not available in the new SCF code).
  
 Example of a PBE calculation on the fractional C cation with 5.3 electrons: Example of a PBE calculation on the fractional C cation with 5.3 electrons:
Line 1454: Line 1474:
 geom={C} geom={C}
 fracocca,4,0.3 fracocca,4,0.3
-{uks,pbe,matrix=0;wf,6,0,2}+{uks,pbe,matrix=0,old; wf,6,0,2}
 </code> </code>
 Example of a RSH calculation on the fractional CO anion with 14.8 electrons: Example of a RSH calculation on the fractional CO anion with 14.8 electrons:
Line 1469: Line 1489:
 fracocca,8,0.8 fracocca,8,0.8
 {int;erf,0.5} {int;erf,0.5}
-{uks,ecerfpbe,exerfpbe,matrix=0;rangehybrid;wf,15,0,1}+{uks,ecerfpbe,exerfpbe,matrix=0,old; rangehybrid;wf,15,0,1}
 </code> </code>
 Example of a HF calculation on the H atom with 0.5 alpha electron and 0.5 beta electron: Example of a HF calculation on the H atom with 0.5 alpha electron and 0.5 beta electron:
Line 1479: Line 1499:
 fracocca,1,.5 fracocca,1,.5
 fracoccb,1,.5 fracoccb,1,.5
-{uhf;wf,2,0,0}+{uhf,old; wf,2,0,0}
 </code> </code>
-Subsequent RPA correlation calculations (see [[the_density_functional_program#random-phase_approximation|here]]) will be automatically done with fractional orbital occupation numbers.+Subsequent RPA correlation calculations (see [[kohn-sham_random-phase_approximation#random-phase_approximation_rpatddft_program|here]]) will be automatically done with fractional orbital occupation numbers.
  
 References\\ References\\