Repository: ifilot/pyqint (language: python)

Issue:
Possible bugs about GTO gradient calculation and Bragg-Slater radius
Hi Ivo,

I recently tried to follow up on your codes to learn DFT, they are very helpful for beginners, thank you! 

I found there are some codes I am confused and they might be bugs:
1. [Here](https://github.com/ifilot/pyqint/blob/88821cbec18f82d53fd0e5d7a051e81e745aca36/pyqint/cgf.cpp#L100-L108), you use $(x-x_0)^{l-1}$ represent the gradient for the angular term, shouldn't it be $l\times (x-x_0)^{l-1}$?
2.  [Here](https://gitlab.tue.nl/ifilot/pydft/-/blob/master/pydft/atomicgrid.py?ref_type=heads#L438) in your `pydft` package, it seems that you get Bragg-Slater radius by directly passing the atom index, but since your Bragg-Slater radius list is started from hydrogen instead of atom placeholder, so the true way should be `bragg_slater.BSRADII[self.__aidx-1]`.


Discussion / hints:

Fix patch:
diff --git a/COMPILATION.md b/COMPILATION.md
index 21360dd..a11d95c 100644
--- a/COMPILATION.md
+++ b/COMPILATION.md
@@ -11,7 +11,7 @@ pipx run cibuildwheel --only cp310-manylinux_x86_64
 To install the `whl` file
 
 ```bash
-pip3 install wheelhouse/pyqint-0.17.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
+pip3 install wheelhouse/pyqint-0.17.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
 ```
 
 and to locally test
diff --git a/meta.yaml b/meta.yaml
index 6c6dac8..ec52621 100644
--- a/meta.yaml
+++ b/meta.yaml
@@ -1,6 +1,6 @@
 package:
   name: "pyqint"
-  version: "0.17.2"
+  version: "0.17.3"
 
 source:
   path: .
diff --git a/pyqint/_version.py b/pyqint/_version.py
index 6109b9e..81b4fea 100644
--- a/pyqint/_version.py
+++ b/pyqint/_version.py
@@ -1,2 +1,2 @@
-__version__ = "0.17.2"
+__version__ = "0.17.3"
 
diff --git a/pyqint/cgf.cpp b/pyqint/cgf.cpp
index a95fa81..a5b8bd5 100644
--- a/pyqint/cgf.cpp
+++ b/pyqint/cgf.cpp
@@ -84,6 +84,8 @@ const double GTO::get_amp(const Vec3& r) const {
  * @return gradient
  */
 Vec3 GTO::get_grad(const Vec3& r) const {
+    // calculate exponential term and its product with the cartesian terms
+    // for x,y,z components
     const double ex = std::exp(-this->alpha * std::pow(r[0]-this->position[0],2));
     const double fx = std::pow(r[0] - this->position[0], this->l) * ex;
 
@@ -93,20 +95,24 @@ Vec3 GTO::get_grad(const Vec3& r) const {
     const double ez = std::exp(-this->alpha * std::pow(r[2]-this->position[2],2));
     const double fz = std::pow(r[2] - this->position[2], this->n) * ez;
 
+    // calculate first derivative of the exponential term
     double gx = -2.0 * this->alpha * (r[0]-this->position[0]) * fx;
     double gy = -2.0 * this->alpha * (r[1]-this->position[1]) * fy;
     double gz = -2.0 * this->alpha * (r[2]-this->position[2]) * fz;
 
+    // if there is a Cartesian component (l,m,n > 0), apply the product rule
+    // and add the contribution of this term
     if(this->l > 0) {
-        gx += std::pow(r[0] - this->position[0], this->l-1) * ex;
+        gx += this->l * std::pow(r[0] - this->position[0], this->l-1) * ex;
     }
     if(this->m > 0) {
-        gy += std::pow(r[1] - this->position[1], this->m-1) * ey;
+        gy += this->m * std::pow(r[1] - this->position[1], this->m-1) * ey;
     }
     if(this->n > 0) {
-        gz += std::pow(r[2] - this->position[2], this->n-1) * ez;
+        gz += this->n * std::pow(r[2] - this->position[2], this->n-1) * ez;
     }
 
+    // return vector with derivative towards x,y and z
     return Vec3(this->norm * gx * fy * fz,
                 this->norm * fx * gy * fz,
                 this->norm * fx * fy * gz);


Test patch:
diff --git a/tests/test_grad.py b/tests/test_grad.py
index c9312ef..774eba8 100644
--- a/tests/test_grad.py
+++ b/tests/test_grad.py
@@ -1,6 +1,7 @@
 import unittest
 from pyqint import cgf
 import numpy as np
+import itertools
 
 class TestGrad(unittest.TestCase):
     """
@@ -14,7 +15,8 @@ def test_cgf_grad(self):
         """
         h = 1e-4 # set step size for finite difference
 
-        pp = [[0,0,0], [1,0,0], [0,1,0], [0,0,1], [1,1,0], [0,1,1], [1,0,1], [1,1,1]]
+        pp = list(itertools.product(range(4), repeat=3))
+
         for exp in pp:
             for p in pp:
                 l,m,n = exp
@@ -32,7 +34,7 @@ def test_grad_density(self):
         """
         h = 1e-4 # set step size for finite difference
 
-        pp = [[0,0,0], [1,0,0], [0,1,0], [0,0,1], [1,1,0], [0,1,1], [1,0,1], [1,1,1]]
+        pp = list(itertools.product(range(4), repeat=3))
         for exp in pp:
             for p in pp:
                 l,m,n = exp