Change SR integral algorithm for nroots <= 3

Refactor g2e.c Remove GTG integrals and coulerf function
sunqm · Sep 19, 2023 · 5c9fd60 · 5c9fd60
1 parent a88d65e
commit 5c9fd60
Show file tree

Hide file tree

Showing 23 changed files with 3,635 additions and 1,951 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -1,7 +1,7 @@
 cmake_minimum_required (VERSION 3.5)
 project (cint C)
-set(cint_VERSION_MAJOR "5")
-set(cint_VERSION_MINOR "5")
+set(cint_VERSION_MAJOR "6")
+set(cint_VERSION_MINOR "0")
 set(cint_VERSION_PATCH "0")
 set(cint_VERSION_TWEAK "0")
 set(cint_VERSION "${cint_VERSION_MAJOR}.${cint_VERSION_MINOR}.${cint_VERSION_PATCH}")
@@ -89,23 +89,17 @@ set(cintSrc
 if(WITH_RANGE_COULOMB)
 # defined in config.h
 #  add_definitions(-DWITH_RANGE_COULOMB)
-  message("Enabled WITH_RANGE_COULOMB")
-  if(WITH_POLYNOMIAL_FIT)
-    set(cintSrc ${cintSrc} src/polyfits.c src/sr_rys_polyfits.c)
-    add_definitions(-DWITH_POLYNOMIAL_FIT)
-    message("Enabled WITH_POLYNOMIAL_FIT")
-  endif(WITH_POLYNOMIAL_FIT)
+#  message("Enabled WITH_RANGE_COULOMB")
 endif(WITH_RANGE_COULOMB)
 
-if(WITH_EIG_ROOTFINDER)
-  add_definitions(-DWITH_EIG_ROOTFINDER)
-  message("Enabled linear equation solver with eigenvalue algorithm")
-endif(WITH_EIG_ROOTFINDER)
+if(WITH_POLYNOMIAL_FIT)
+  set(cintSrc ${cintSrc} src/polyfits.c src/sr_rys_polyfits.c)
+  add_definitions(-DWITH_POLYNOMIAL_FIT)
+  message("Enabled WITH_POLYNOMIAL_FIT")
+endif(WITH_POLYNOMIAL_FIT)
 
 if(WITH_COULOMB_ERF)
-  set(cintSrc ${cintSrc} src/g2e_coulerf.c src/cint2e_coulerf.c)
-  add_definitions(-DWITH_COULOMB_ERF)
-  message("Enabled WITH_COULOMB_ERF")
+  message("WITH_COULOMB_ERF is deprecated since v6.0")
 endif(WITH_COULOMB_ERF)
 
 if(WITH_F12)
@@ -115,12 +109,13 @@ if(WITH_F12)
 endif(WITH_F12)
 
 # Gaussian type geminal code are incorrect
-#if(WITH_GTG)
+if(WITH_GTG)
 #  set(cintSrc ${cintSrc} src/g2e_gtg.c src/cint2e_gtg.c src/cint3c2e_gtg.c
 #    src/cint2c2e_gtg.c)
 #  add_definitions(-DWITH_GTG)
 #  message("Enabled WITH_GTG. Note there are bugs in gtg type integrals")
-#endif(WITH_GTG)
+  message("WITH_GTG is deprecated since v6.0")
+endif(WITH_GTG)
 
 if(WITH_4C1E)
   set(cintSrc ${cintSrc} src/cint4c1e.c src/g4c1e.c)

diff --git a/ChangeLog b/ChangeLog
@@ -1,7 +1,9 @@
+Version 6.0.0 (2023-09-14):
+	* Add polynomial root finder using the eigenvalue algorithm
+	* Fitting Rys roots of SR integrals with Chebyshev polynomials
+
 Version 5.5.0 (2023-09-14):
 	* Fix memory leak issue in int4c1e
-	* Add polynomial root finder with eigenvalue algorithm
-	* Fitting Rys roots of SR integrals with Chebyshev polynomials
 	* New integrals related to LRESC
 	* Supports Gaussian model for int1e_grids
 

diff --git a/include/cint.h.in b/include/cint.h.in
@@ -4,7 +4,8 @@
  * Parameters and function signature for libcint.
  */
 
-#define CINT_VERSION    @cint_VERSION@
+#define CINT_VERSION            @cint_VERSION@
+#define CINT_SOVERSION          @cint_SOVERSION
 
 #cmakedefine I8
 #ifdef I8
@@ -170,7 +171,7 @@ typedef struct {
         union {FINT nfk; FINT grids_offset;};
         union {FINT nfl; FINT ngrids;};
         FINT nf;  // = nfi*nfj*nfk*nfl;
-        FINT _padding;
+        FINT rys_order; // = nrys_roots for regular ERIs. can be nrys_roots/2 for SR ERIs
         FINT x_ctr[4];
 
         FINT gbits;

diff --git a/src/cint2c2e_gtg.c b/src/cint2c2e_gtg.c
diff --git a/src/cint2e.c b/src/cint2e.c
@@ -122,15 +122,14 @@ FINT CINT2e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache, FINT *empt
         double *rij;
         akl = ak[k_prim-1] + al[l_prim-1];
         log_rr_kl = 1.7 - 1.5 * approx_log(akl);
-#ifdef WITH_RANGE_COULOMB
         double omega = env[PTR_RANGE_OMEGA];
         if (omega < 0) {
                 // Normally the factor
                 //    (aj*d/aij+theta*R)^li * (ai*d/aij+theta*R)^lj * pi^1.5/aij^{(li+lj+3)/2}
                 // is a good approximation for polynomial parts in SR-ERIs.
                 //    <~ (aj*d/aij+theta*R)^li * (ai*d/aij+theta*R)^lj * (pi/aij)^1.5
                 //    <~ (d+theta*R)^li * (d+theta*R)^lj * (pi/aij)^1.5
-                if (envs->nrys_roots > 1) {
+                if (envs->rys_order > 1) {
                         double r_guess = 8.;
                         double omega2 = omega * omega;
                         int lij = envs->li_ceil + envs->lj_ceil;
@@ -152,11 +151,6 @@ FINT CINT2e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache, FINT *empt
                         log_rr_kl += lkl * approx_log(sqrt(rr_kl) + 1.);
                 }
         }
-#else
-        if (lkl > 0) {
-                log_rr_kl += lkl * approx_log(sqrt(rr_kl) + 1.);
-        }
-#endif
 
         FINT *idx;
         MALLOC_INSTACK(idx, nf * 3);
@@ -360,10 +354,9 @@ k_contracted: ;
                 CINTg2e_index_xyz(idx, envs); \
         }
 
-#ifdef WITH_RANGE_COULOMB
 #define ADJUST_CUTOFF      \
         double omega = env[PTR_RANGE_OMEGA]; \
-        if (omega < 0 && envs->nrys_roots > 1) { \
+        if (omega < 0 && envs->rys_order > 1) { \
                 double r_guess = 8.; \
                 double omega2 = omega * omega; \
                 int lij = envs->li_ceil + envs->lj_ceil; \
@@ -383,9 +376,6 @@ k_contracted: ;
                                 (dist_kl+theta*r_guess+1.)/(dist_kl+1.)); \
                 } \
         }
-#else
-#define ADJUST_CUTOFF
-#endif
 
 #define SET_RIJ(I,J)    \
         if (pdata_##I##J->cceij > e##I##J##cutoff) { \

diff --git a/src/cint2e_coulerf.c b/src/cint2e_coulerf.c
diff --git a/src/cint2e_gtg.c b/src/cint2e_gtg.c
diff --git a/src/cint3c2e.c b/src/cint3c2e.c
@@ -105,9 +105,8 @@ FINT CINT3c2e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache, FINT *em
         double expij, cutoff;
         double *rij;
         double *rkl = envs->rk;
-#ifdef WITH_RANGE_COULOMB
         double omega = env[PTR_RANGE_OMEGA];
-        if (omega < 0 && envs->nrys_roots > 1) {
+        if (omega < 0 && envs->rys_order > 1) {
                 double r_guess = 8.;
                 double omega2 = omega * omega;
                 int lij = envs->li_ceil + envs->lj_ceil;
@@ -123,7 +122,6 @@ FINT CINT3c2e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache, FINT *em
                         expcutoff += envs->lk_ceil * approx_log(theta*r_guess+1.);
                 }
         }
-#endif
 
         FINT *idx;
         MALLOC_INSTACK(idx, nf * 3);
@@ -279,10 +277,9 @@ i_contracted: ;
                 CINTg2e_index_xyz(idx, envs); \
         }
 
-#ifdef WITH_RANGE_COULOMB
 #define ADJUST_CUTOFF      \
         double omega = env[PTR_RANGE_OMEGA]; \
-        if (omega < 0 && envs->nrys_roots > 1) { \
+        if (omega < 0 && envs->rys_order > 1) { \
                 double r_guess = 8.; \
                 double omega2 = omega * omega; \
                 int lij = envs->li_ceil + envs->lj_ceil; \
@@ -298,9 +295,6 @@ i_contracted: ;
                         expcutoff += envs->lk_ceil * approx_log(theta*r_guess+1.); \
                 } \
         }
-#else
-#define ADJUST_CUTOFF
-#endif
 
 #define SET_RIJ    \
         if (pdata_ij->cceij > expcutoff) { \

diff --git a/src/cint3c2e_gtg.c b/src/cint3c2e_gtg.c
diff --git a/src/fmt.c b/src/fmt.c
@@ -334,15 +334,13 @@ void fmt_erfc_like(double *f, double t, double lower, int m)
 
         int i;
         double lower2 = lower * lower;
-#ifdef WITH_RANGE_COULOMB
         // F[m] < .5*sqrt(pi/t) * erfc(low*tt)
         if (t * lower2 > ERFC_bound) {
                 for (i = 0; i <= m; i++) {
                         f[i] = 0;
                 }
                 return;
         }
-#endif
 
         if (t < TURNOVER_POINT[m]) {
                 fmt1_erfc_like(f, t, lower, m);
@@ -372,15 +370,13 @@ void fmt_lerfc_like(long double *f, long double t, long double lower, int m)
 
         int i;
         long double lower2 = lower * lower;
-#ifdef WITH_RANGE_COULOMB
         // F[m] < .5*sqrt(pi/t) * erfc(low*tt)
         if (t * lower2 > ERFC_bound) {
                 for (i = 0; i <= m; i++) {
                         f[i] = 0;
                 }
                 return;
         }
-#endif
 
         if (t < TURNOVER_POINT[m]) {
                 fmt1_lerfc_like(f, t, lower, m);
@@ -492,15 +488,13 @@ void fmt_qerfc_like(__float128 *f, __float128 t, __float128 lower, int m)
 
         int i;
         __float128 lower2 = lower * lower;
-#ifdef WITH_RANGE_COULOMB
         // F[m] < .5*sqrt(pi/t) * erfc(low*tt)
         if (t * lower2 > ERFC_bound) {
                 for (i = 0; i <= m; i++) {
                         f[i] = 0;
                 }
                 return;
         }
-#endif
 
         if (t < TURNOVER_POINT[m]) {
                 fmt1_qerfc_like(f, t, lower, m);