MRChemSoft / mrchem

MultiResolution Chemistry
GNU Lesser General Public License v3.0
30 stars 21 forks source link

Add Hirshfeld charge analysis #494

Closed moritzgubler closed 2 months ago

moritzgubler commented 2 months ago

I implemented Hirshfeld charge partitioning. A method to compute partial charges of atoms in molecules which is implemented in many electronic structure code. The Hirshfeld charge of atom $i$ is defined as: $$q_i=Z_i-\int\frac{\rho^0_i(\mathbf{r})}{\sum_j\rho^0_j(\mathbf{r})}\rho(r)d\mathbf{r}$$ where $Z_i$ is the atomic number of element $i$, $\rho$ is the molecular density and $\rho_i^0$ is density of $i$ as an isolated atom. The spherically symmetric reference charge densities were obtained with an LDA calculation and are tabulated in share/hirshfeld/lda/{element_name}.density

To improve the numerical stability, the logarithm of the density is interpolated internally and the logsumexp trick was used to get rid of 0/0 issues in regions far away from all nuclei in the function in the integral above.

moritzgubler commented 2 months ago

Are the SAD input densities also spherically symmetric? If you think the radial interpolator is useful for other people too, we should probably put it in a different directory do you have some suggestions?

stigrj commented 2 months ago

Currently the SAD densities are not symmetric, because they are precomputed using unrestricted LDA in a tiny 3-21G basis, but I think it would be better if they were symmetric. Perhaps move it to a folder under utils?

moritzgubler commented 2 months ago

Yes, in that case I think it is better when they are symmetric. I can polish the interpolator a bit, implement some more interpolation modes and move it to utils.

moritzgubler commented 2 months ago

I moved specific parts of my implementation into more accessible locations:

I also added another extrapolation mode in the poynomial interpolator that is useful for interpolating densities directly: No extrapolation to the left and returning 0 if x > xmax.

The building blocks to replace the SAD densities by our new densities should now be in place and it should be quite easy to implement it in the future