The python code is here.

I tried to do a DFT calculation for H2 Restricted Hartree-Fock with Local Density Approximation. The optimization is a scf-loop with Lagrange-multiplier only. I use the same atomic orbitals as the STO3G ab initio computation counterpart. I get the following Energy plot : enter image description here

My guess is that the exchange(-correlation) energy/operator computation is not good since I am using Riemann integrals (uniform cubes formed with a uniform point grid) to compute it.

How are computed exchange-correlation integrals usually? Via Gauss quadrature? Analytic integration if possible?

If you have any questions, dont hesitate.


1 Answer 1


The exchange-correlation energy and potential of a molecular calculation are usually integrated on a grid that is dense near nuclei (to integrate the nuclear cusps accurately) and sparse far from nuclei (to save computational cost). Since the grid spacing is non-uniform, the grid points near the nuclei have to receive smaller weights, and the points far from the nuclei receive larger weights.

Firstly, consider an atomic calculation. It is very natural to perform the integration in polar coordinates, using the position of the atomic nucleus as the origin. In many (maybe most) quantum chemistry programs, the radial integration is done with Gaussian quadrature, after a nonlinear coordinate transformation that makes the radial grid very dense near $r=0$ ($r$ is the distance of the grid point to the nucleus) and sparse when $r$ is large. The angular (i.e. $\theta$ and $\phi$) integration is usually done with Lebedev quadrature, which is an efficient quadrature that exactly integrates spherical harmonics up to a high degree, just as the Gaussian quadrature exactly integrates polynomials up to a high degree. Integrating under polar coordinates not only utilizes and keeps all relevant spherical symmetries, but also allows us to separate out the radial degree of freedom, and therefore apply well-known techniques to put more grid points near the nuclei while still being able to assign suitable weights to the grid points.

Doing the integrations for molecules is more complicated, since it is no longer possible to separate out the radial degrees of freedom. The usual strategy is to partition the integrand as a sum of atomic contributions, using some smoothened analog of the well-known Voronoi partitioning, and integrate each contribution using the atomic grid of the corresponding atom (which is constructed according to the previous paragraph). One early example is the Becke partitioning, which is later improved by the Stratmann-Scuseria-Frisch (SSF) partitioning, that has better numerical stability and costs much less in calculating the weights of the molecular grid points.

There are a few more techniques for reducing the number of grid points for a given integration accuracy, such as using small Lebedev grids for very small and very large $r$, and large Lebedev grids for intermediate $r$ (where an accurate angular integration is more important). These kind of grid optimizations can go up to quite elaborate levels, as exemplified by the newest ORCA grids. Some earlier protocols for optimizing the molecular grids can be found in the reference list of this paper as well.

  • $\begingroup$ Thank you for this detailed response. This gives a general recipe and was clear. I'll read and try the methods you mentioned. In this case, which is the simplest method to implement? Since, I have the intention to have a test case that works. I expected H2 to be the simplest molecule to start learning about DFT without any particular mesh techniques. $\endgroup$
    – mle
    May 6, 2023 at 19:09
  • 1
    $\begingroup$ For diatomic molecules, especially symmetric diatomic molecules, special axially symmetric grids are available. Once you want to calculate non-linear molecules, however, the simplest way would be Gaussian radial quadrature + Lebedev angular quadrature + Becke partitioning. The SSF partitioning is only slightly more complex than the Becke partitioning (just a few more lines) and is also worth trying. The really complicated technical stuff is how to efficiently prune the grids, and this can be ignored at the initial stage of program development. $\endgroup$
    – wzkchem5
    May 6, 2023 at 20:03

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .