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);
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