High-order differentiation of the density with respect to the position

[RaphaelRobidas]

Hello,

In the documentation (https://chemtools.org/sci_doc_dftbased_it.html), the steric potential and steric charge density are presented. They are not implemented directly in the code, but I believe I have managed to calculate the steric potential using density, its gradient and the Laplacian. However, higher order differentiations are necessary to compute the steric charge density, and I have not found how to do that in chemtools (non-numerically). Is it possible or am I missing the point here?

Thanks!

[PaulWAyers]

The workflow for this would use gbasis which can evaluate arbitrary derivatives of the electron density and one-electron reduced density matrix. I do not think that a front end linking those capabilities into ChemTools is developed yet. (@FarnazH would know for sure.) But it would be a nice project.

I'd be happy to try to help you work it out if you are interested.

[RaphaelRobidas]

It seems like the number of grid points isn't the only important factor: the grid must start very close to the nuclei, most likely because of the wild fluctuations at short distances.

Here are some examples of grid settings and the corresponding integral of the steric charge for neon:

exp:0.001:5:600:26
Integral: -8.642083283836156

exp:0.00001:5:600:26
Integral: -0.0002992803651711061

exp:0.0000001:5:600:26
Integral: -0.0002606104362368989

exp:0.0000000001:5:600:26
Integral: -0.00022561087383682669

exp:0.0000000001:8:600:26
Integral: 9.304317069567754e-11

We do get an integral of around 0 as expected. I don't have much more insight on why that must be the case however, but it would be interesting to have an intuition for it indeed.

[PaulWAyers]

So the integral is going to zero well. That's good. The ill-conditioning is not unexpected, and occurs for exactly the reason you suspect. So I think that's a good indication that the steric charge is implemented correctly.

I'm still a bit befuddled by the integration, but no reason to look a gift horse in the mouth.

[RaphaelRobidas]

I managed to get an integral of close to 0 (0.054) with ethanol as well, but it required absolutely aug-cc-pVQZ and a pretty fine grid (exp:0.0000000001:8:1000:86).

It also seems that random rotations of the grid throw off the integral by quite a lot. Otherwise, I think we have it working alright.

Did you manage to install it on ComputeCanada/the Digital Research Alliance's clusters? Conda is not allowed and I did not have any luck installing chemtools without it on my local machine. The calculations necessary for the ethanol molecule took a couple hours on my workstation, so calculations at a larger scale would be better done on a cluster.

[leila-pujal]

Hi @RaphaelRobidas, looks you and @PaulWAyers have made a lot of progress on this issue, thanks a lot. Regarding installing chemtools on Compute Canada (aka Digital Research Alliance) you are right, you are not allowed to use conda but you can use python virtual enviroments. You get more or less the same functionality you have with conda except you have to use pip to install packages instead of using conda-forge. I added below some instructions to install chemtools along with horton using a virtual enviroment. You should be able to get the same working version of chemtools as the one you have installed using conda on your local machine. I tested the instructions they worked for me but let me know if you have any problem. I got some expected test fails but you shouldn't worry too much about them, you can double-check everything is working fine by reproducing some of your results.

# make a Python 3 environment on Compute Canada
module load python/3.9 libxc/3.0.0 libint/2.6.0 git-lfs/2.11.0
virtualenv --no-download envs/env_denspart_py39
# Next line activates the enviroment and is the analog as conda activate
source ~/envs/env_denspart_py39/bin/activate
pip install numpy cython scipy h5py matplotlib nose --no-index
# if you already don't have ${HOME}/modules directory, create one
mkdir ~/modules

# install HORTON
cd ~/modules
curl -kLO https://github.com/theochem/horton/releases/download/2.3.0/horton-2.3.0.tar.gz
tar -xvf horton-2.3.0.tar.gz && mv horton horton-2.3.0
cd horton-2.3.0
./setup.py install
# You need to export this variable so horton knows where the data is. I am not sure if the line below is a permanent thing or is only for a particular terminal session. If you realise you need to export this variable every time just add this to your .bashrc file.
export HORTONDATA=${HOME}/modules/horton-2.3.0/data
cd ~/
# run tests (some tests fail, but most should succeed)
cd ~; nosetests -v horton

# install ChemTools
cd ~/modules
git clone https://github.com/theochem/chemtools.git
cd chemtools
pip install -e .
# run test (all tests should pass successfully except 1)
nosetests -v chemtools

[RaphaelRobidas]

Thanks @leila-pujal, your script worked after some modifications. I had to install qc-iodata, grid and gbasis from source and downgrade numpy to precisely version 1.23.0. With this, chemtools runs on the cluster.

@PaulWAyers I compared the partitioned steric charges of LiH obtained from my calculations to those in the reference paper. They are quite off, even though the integral converges to relatively close to zero (with aug-cc-pV5Z-quality SCF and 1200 radial/170 spherical grid points). I noticed that the small density amount added to avoid the long range abnormal behavior actually influences the partitioning (and steric charge profile) a lot, although not the overall integral.

Here is an example with neon. I compared the profile of the steric charge with different $\delta \rho$ values:

With a value of $10^{-7}$ or above, the integral is essentially zero. However, the profile changes considerably. These changes are probably very different for different atoms, resulting in extremely variable partial (partitioned) steric charges.

[PaulWAyers]

I thought some about this, and the unit tests we can use. Obviously one of them is using finite differences.

But there are others, which are more global.

1. Test for the Steric Potential: Because the Weizsacker kinetic-energy function is homogeneous of degree one with respect to density scaling,

$$ T_{\text{Weizsacker}}\left[\rho \right] = \int \rho(\mathbf{r}) v_{\text{steric}} (\mathbf{r}) d \mathbf{r} $$

2. Test for the Steric force: Because the Weizsacker kinetic energy is homogenous of degree two with respect to coordinate scaling of the density,

$$ T_{\text{Weizsacker}}\left[\rho \right] =\tfrac{1}{2} \int \rho(\mathbf{r}) \mathbf{r} \cdot \nabla v_{\text{steric}} (\mathbf{r}) d \mathbf{r} $$

3. Test for the Steric charge: And from the (Green's function for the Laplace equation),

$$ v_{\text{steric}} = \int \frac{q_{\text{steric}} (\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d \mathbf{r}' $$

[PaulWAyers]

P.S. For the Steric Charge, I think I'd compute the Steric Hessian (the full second derivative matrix of the Steric Potential) and then take the trace of that matrix. The reason is that the Steric Potential is relatively new, and topological eigenvalue analysis of it, similar to what is done for the density and other scalar fields in QTAIM/QCT, would be interesting.

[PaulWAyers]

P.S. I asked Shubin Liu about looking at the steric potential Hessian and he thought it was interesting.

[PaulWAyers]

I heard back from Shubin and Tian Lu. (Thanks to both of them.) The issue is that these plots were generated with multiwfn, and multiwfn regularizes the functions in a way described in the below figure. If they don't do that, their plots (also below) look like what you're getting.

In ChemTools we are more selective about how we avoid divide-by-zero errors. Generally we identify divide-by-zero errors and then either trim the function (leaving it undefined where there are divide-by-zero-errors) or shift to an asymptotic form when the denominator of expressions becomes too small. That's something we should keep in mind in this implementation too. Fortunately, because the asymptotic decay here is so well-defined, shifting to the "right" asymptotic expression is not difficult. And it won't really matter until you are far from all the nuclei. It is important (at least potentially) for integrated properties, however.