Fixing bug in FB and adding Fgamma-cached variant of repulsion integr…

…al evaluation

examples/foster_boys_quick.py

@@ -0,0 +1,27 @@
+from pyqint import Molecule, HF, COHP, FosterBoys
+import numpy as np
+
+d = 1.145414
+mol = Molecule()
+mol.add_atom('C', 0.0, 0.0, -d/2, unit='angstrom')
+mol.add_atom('O', 0.0, 0.0,  d/2, unit='angstrom')
+
+res = HF().rhf(mol, 'sto3g')
+cohp = COHP(res).run(res['orbc'], 0, 1)
+
+resfb = FosterBoys(res).run(nr_runners=10)
+cohp_fb = COHP(res).run(resfb['orbc'], 0, 1)
+
+print('COHP values of canonical Hartree-Fock orbitals')
+for i,(e,chi) in enumerate(zip(res['orbe'], cohp)):
+    print('%3i %12.4f %12.4f' % (i+1,e,chi))
+print()
+
+print('COHP values after Foster-Boys localization')
+for i,(e,chi) in enumerate(zip(resfb['orbe'], cohp_fb)):
+    print('%3i %12.4f %12.4f' % (i+1,e,chi))
+print()
+
+print('Sum of COHP coefficient canonical orbitals: ', np.sum(cohp[:7]))
+print('Sum of COHP coefficient Foster-Boys orbitals: ', np.sum(cohp_fb[:7]))
+print('Result FB: ', resfb['r2start'], resfb['r2final'])

meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: "pyqint"
-  version: "0.13.0.0"
+  version: "0.13.1.0"
 
 source:
   path: .

pyqint/_version.py

@@ -1,2 +1,2 @@
-__version__ = "0.14.0.0"
+__version__ = "0.13.1.0"
 

pyqint/foster_boys.py

@@ -41,7 +41,7 @@ def run(self, nr_runners=1):
             if res['r2final'] > bestr2:
                 bestres = res
 
-        return res
+        return bestres
 
     def __single_runner(self):
         """

pyqint/gamma.cpp

@@ -38,7 +38,7 @@
 
 #include "gamma.h"
 
-double GammaInc::Fgamma(const double m, double x) const {
+double GammaInc::Fgamma(double m, double x) const {
     double tiny = 0.00000001;
     x = std::max(std::abs(x),tiny);
     double val = gamm_inc(m+0.5,x);
@@ -58,7 +58,7 @@ double GammaInc::Fgamma(const double m, double x) const {
  *
  * returns double value of the incomplete Gamma Function
  */
-double GammaInc::gamm_inc(const double a, const double x) const {
+double GammaInc::gamm_inc(double a, double x) const {
     double gammap = gammp(a,x);
     double gln = gammln(a);
     return exp(gln) * gammap;
@@ -76,7 +76,7 @@ double GammaInc::gamm_inc(const double a, const double x) const {
  *
  * returns double value of the incomplete Gamma Function
  */
-double GammaInc::gammp(const double a, double x) const {
+double GammaInc::gammp(double a, double x) const {
     static const int ASWITCH = 100;
 
     if(x < 0.0 || a <= 0.0) {
@@ -104,7 +104,7 @@ double GammaInc::gammp(const double a, double x) const {
  *
  * returns double value of the incomplete Gamma Function
  */
-double GammaInc::gser(const double a, const double x) const {
+double GammaInc::gser(double a, double x) const {
     double EPS = std::numeric_limits<double>::epsilon();
 
     double sum, del, ap;
@@ -122,10 +122,10 @@ double GammaInc::gser(const double a, const double x) const {
     }
 }
 
-double GammaInc::gammln(const double xx) const {
+double GammaInc::gammln(double xx) const {
     int j;
     double x, tmp, y, ser;
-    static const double cof[14]={57.1562356658629235,-59.5979603554754912,
+    static double cof[14]={57.1562356658629235,-59.5979603554754912,
         14.1360979747417471,-0.491913816097620199,.339946499848118887e-4,
         .465236289270485756e-4,-.983744753048795646e-4,.158088703224912494e-3,
         -.210264441724104883e-3,.217439618115212643e-3,-.164318106536763890e-3,
@@ -152,7 +152,7 @@ double GammaInc::gammln(const double xx) const {
  *
  * returns double value of the incomplete Gamma Function
  */
-double GammaInc::gcf(const double a, const double x) const {
+double GammaInc::gcf(double a, double x) const {
     double EPS = std::numeric_limits<double>::epsilon();
     double FPMIN = std::numeric_limits<double>::min() / EPS;
 

pyqint/gamma.h

@@ -44,7 +44,7 @@
 
 class GammaInc {
 public:
-    double Fgamma(const double m, double x) const;
+    double Fgamma(double m, double x) const;
 
     /*
      * @fn gamm_inc
@@ -59,7 +59,7 @@ class GammaInc {
      *
      * returns double value of the incomplete Gamma Function
      */
-    double gamm_inc(const double a, const double x) const;
+    double gamm_inc(double a, double x) const;
 
     /*
      * @fn gamm_inc
@@ -73,7 +73,7 @@ class GammaInc {
      *
      * returns double value of the incomplete Gamma Function
      */
-    double gammp(const double m, const double x) const;
+    double gammp(double m, double x) const;
 
 private:
     /*
@@ -85,10 +85,10 @@ class GammaInc {
      *
      * returns double value of the incomplete Gamma Function
      */
-    double gser(const double a, const double x) const;
+    double gser(double a, double x) const;
 
 
-    double gammln(const double xx) const;
+    double gammln(double xx) const;
 
     /*
      * @fn gcf
@@ -99,7 +99,7 @@ class GammaInc {
      *
      * returns double value of the incomplete Gamma Function
      */
-    double gcf(const double a, const double x) const;
+    double gcf(double a, double x) const;
 
     /*
      * @fn gammpapprox

pyqint/integrals.cpp

@@ -849,10 +849,10 @@ double Integrator::repulsion_deriv(const CGF &cgf1, const CGF &cgf2, const CGF &
  */
 double Integrator::repulsion(const GTO &gto1, const GTO &gto2, const GTO &gto3, const GTO &gto4) const {
 
-    double rep = repulsion(gto1.get_position(), gto1.get_l(), gto1.get_m(), gto1.get_n(), gto1.get_alpha(),
-                           gto2.get_position(), gto2.get_l(), gto2.get_m(), gto2.get_n(), gto2.get_alpha(),
-                           gto3.get_position(), gto3.get_l(), gto3.get_m(), gto3.get_n(), gto3.get_alpha(),
-                           gto4.get_position(), gto4.get_l(), gto4.get_m(), gto4.get_n(), gto4.get_alpha());
+    double rep = this->repulsion(gto1.get_position(), gto1.get_l(), gto1.get_m(), gto1.get_n(), gto1.get_alpha(),
+                                 gto2.get_position(), gto2.get_l(), gto2.get_m(), gto2.get_n(), gto2.get_alpha(),
+                                 gto3.get_position(), gto3.get_l(), gto3.get_m(), gto3.get_n(), gto3.get_alpha(),
+                                 gto4.get_position(), gto4.get_l(), gto4.get_m(), gto4.get_n(), gto4.get_alpha());
 
     return rep;
 }
@@ -1096,6 +1096,23 @@ double Integrator::A_term(const int i, const int r, const int u, const int l1, c
             std::pow(0.25/gamma,r+u)/boost::math::factorial<double>(r)/boost::math::factorial<double>(u)/boost::math::factorial<double>(i-2*r-2*u);
 }
 
+/**
+ * @brief Performs nuclear integral evaluation
+ *
+ * @param vec3 a            Center of the Gaussian orbital of the first GTO
+ * @param unsigned int l1   Power of x component of the polynomial of the first GTO
+ * @param unsigned int m1   Power of y component of the polynomial of the first GTO
+ * @param unsigned int n1   Power of z component of the polynomial of the first GTO
+ * @param double alpha1     Gaussian exponent of the first GTO
+ * @param vec3 b            Center of the Gaussian orbital of the second GTO
+ * @param unsigned int l2   Power of x component of the polynomial of the second GTO
+ * @param unsigned int m2   Power of y component of the polynomial of the second GTO
+ * @param unsigned int n2   Power of z component of the polynomial of the second GTO
+ * @param double alpha2     Gaussian exponent of the second GTO
+ * @param vec3 c
+ *
+ * @return double value of the nuclear integral
+ */
 double Integrator::repulsion(const vec3 &a, const int la, const int ma, const int na, const double alphaa,
                              const vec3 &b, const int lb, const int mb, const int nb, const double alphab,
                              const vec3 &c, const int lc, const int mc, const int nc, const double alphac,
@@ -1132,6 +1149,70 @@ double Integrator::repulsion(const vec3 &a, const int la, const int ma, const in
                  std::exp(-alphac*alphad*rcd2/gamma2)*sum;
 }
 
+/**
+ * @brief Performs nuclear integral evaluation
+ *
+ * This function uses function-level caching of the Fgamma function
+ *
+ * @param vec3 a            Center of the Gaussian orbital of the first GTO
+ * @param unsigned int l1   Power of x component of the polynomial of the first GTO
+ * @param unsigned int m1   Power of y component of the polynomial of the first GTO
+ * @param unsigned int n1   Power of z component of the polynomial of the first GTO
+ * @param double alpha1     Gaussian exponent of the first GTO
+ * @param vec3 b            Center of the Gaussian orbital of the second GTO
+ * @param unsigned int l2   Power of x component of the polynomial of the second GTO
+ * @param unsigned int m2   Power of y component of the polynomial of the second GTO
+ * @param unsigned int n2   Power of z component of the polynomial of the second GTO
+ * @param double alpha2     Gaussian exponent of the second GTO
+ * @param vec3 c
+ *
+ * @return double value of the nuclear integral
+ */
+double Integrator::repulsion_fgamma_cached(const vec3 &a, const int la, const int ma, const int na, const double alphaa,
+                                           const vec3 &b, const int lb, const int mb, const int nb, const double alphab,
+                                           const vec3 &c, const int lc, const int mc, const int nc, const double alphac,
+                                           const vec3 &d, const int ld, const int md, const int nd, const double alphad) const {
+
+    static const double pi = boost::math::constants::pi<double>();
+    double rab2 = (a-b).squaredNorm();
+    double rcd2 = (c-d).squaredNorm();
+
+    vec3 p = gaussian_product_center(alphaa, a, alphab, b);
+    vec3 q = gaussian_product_center(alphac, c, alphad, d);
+
+    double rpq2 = (p-q).squaredNorm();
+
+    double gamma1 = alphaa + alphab;
+    double gamma2 = alphac + alphad;
+    double delta = 0.25 * (1.0 / gamma1 + 1.0 / gamma2);
+
+    std::vector<double> bx = B_array(la, lb, lc, ld, p[0], a[0], b[0], q[0], c[0], d[0], gamma1, gamma2, delta);
+    std::vector<double> by = B_array(ma, mb, mc, md, p[1], a[1], b[1], q[1], c[1], d[1], gamma1, gamma2, delta);
+    std::vector<double> bz = B_array(na, nb, nc, nd, p[2], a[2], b[2], q[2], c[2], d[2], gamma1, gamma2, delta);
+
+    // pre-calculate Fgamma values and cache them inside a vector
+    const unsigned int imax = la+lb+lc+ld;
+    const unsigned int jmax = ma+mb+mc+md;
+    const unsigned int kmax = na+nb+nc+nd;
+    std::vector<double> Fgamma_lookup(imax + jmax + kmax + 1, 0.0);
+    for(unsigned int i=0; i<Fgamma_lookup.size(); i++) {
+        Fgamma_lookup[i] = this->gamma_inc.Fgamma(i,0.25*rpq2/delta);
+    }
+
+    double sum = 0.0;
+    for(int i=0; i<=(la+lb+lc+ld); i++) {
+        for(int j=0; j<=(ma+mb+mc+md); j++) {
+            for(int k=0; k<=(na+nb+nc+nd); k++) {
+                sum += bx[i]*by[j]*bz[k]*Fgamma_lookup[i+j+k];
+            }
+        }
+    }
+
+    return 2.0 * std::pow(pi,2.5)/(gamma1*gamma2*std::sqrt(gamma1+gamma2))*
+                 std::exp(-alphaa*alphab*rab2/gamma1)*
+                 std::exp(-alphac*alphad*rcd2/gamma2)*sum;
+}
+
 std::vector<double> Integrator::B_array(const int l1,const int l2,const int l3,const int l4,
         const double p, const double a, const double b, const double q, const double c, const double d,
         const double g1, const double g2, const double delta) const {

pyqint/integrals.h

@@ -593,6 +593,32 @@ class Integrator {
                      const vec3 &c, const int lc, const int mc, const int nc, const double alphac,
                      const vec3 &d, const int ld, const int md, const int nd, const double alphad) const;
 
+    /**
+     * @brief Performs nuclear integral evaluation including caching of Fgamma
+     *
+     * This function uses function-level caching of the Fgamma function; this implementation
+     * was suggested in https://github.com/ifilot/hfcxx/issues/8, but explicit unit testing
+     * actually shows not appreciable difference in speed.
+     *
+     * @param vec3 a            Center of the Gaussian orbital of the first GTO
+     * @param unsigned int l1   Power of x component of the polynomial of the first GTO
+     * @param unsigned int m1   Power of y component of the polynomial of the first GTO
+     * @param unsigned int n1   Power of z component of the polynomial of the first GTO
+     * @param double alpha1     Gaussian exponent of the first GTO
+     * @param vec3 b            Center of the Gaussian orbital of the second GTO
+     * @param unsigned int l2   Power of x component of the polynomial of the second GTO
+     * @param unsigned int m2   Power of y component of the polynomial of the second GTO
+     * @param unsigned int n2   Power of z component of the polynomial of the second GTO
+     * @param double alpha2     Gaussian exponent of the second GTO
+     * @param vec3 c
+     *
+     * @return double value of the nuclear integral
+     */
+    double repulsion_fgamma_cached(const vec3 &a, const int la, const int ma, const int na, const double alphaa,
+                                   const vec3 &b, const int lb, const int mb, const int nb, const double alphab,
+                                   const vec3 &c, const int lc, const int mc, const int nc, const double alphac,
+                                   const vec3 &d, const int ld, const int md, const int nd, const double alphad) const;
+
     /**
      * @brief Calculates one dimensional overlap integral
      *