Release v4.4.0

Fix int1e_grids code generator
Refactor cint1e_grids
Refactor cart2sph

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required (VERSION 3.5)
 project (cint C)
 set(cint_VERSION_MAJOR "4")
-set(cint_VERSION_MINOR "3")
+set(cint_VERSION_MINOR "4")
 set(cint_VERSION_PATCH "0")
 set(cint_VERSION_TWEAK "0")
 set(cint_VERSION "${cint_VERSION_MAJOR}.${cint_VERSION_MINOR}.${cint_VERSION_PATCH}")

ChangeLog

@@ -1,5 +1,10 @@
+Version 4.4.0  (2021-05-04):
+	* Fix bugs in c2s transfomration for int1e_grids
+	* Fix bugs int1e_grids code generator
+	* Refactor cart2sph
+	* Refactor int1e drivers
 Version 4.3.0  (2021-04-24):
-	* Add int1e_grids
+	* Add new integral type int1e_grids
 	* Fix cache size type
 Version 4.1.3  (2021-04-13):
 	* Fix memory address int32 overflow for heavily contracted basis

README

@@ -1,8 +1,8 @@
 libcint
 =======
 
-version 4.3.0
-2021-04-24
+version 4.4.0
+2021-05-04
 
 
 What is libcint

include/cint.h.in

@@ -22,21 +22,27 @@
 #define CACHE_SIZE_T FINT
 #endif
 
-#define PTR_EXPCUTOFF   0
-#define PTR_COMMON_ORIG 1
-#define PTR_RINV_ORIG   4
-#define PTR_RINV_ZETA   7
+// global parameters in env
+// Overall cutoff for integral prescreening, value needs to be ~ln(threshold)
+#define PTR_EXPCUTOFF           0
+// R_C of (r-R_C) in dipole, GIAO operators
+#define PTR_COMMON_ORIG         1
+// R_O in 1/|r-R_O|
+#define PTR_RINV_ORIG           4
+// ZETA parameter for Gaussian charge distribution (Gaussian nuclear model)
+#define PTR_RINV_ZETA           7
 // omega parameter in range-separated coulomb operator
 // LR interaction: erf(omega*r12)/r12 if omega > 0
 // SR interaction: erfc(omega*r12)/r12 if omega < 0
-#define PTR_RANGE_OMEGA 8
-// Yukawa potential and slater-type geminal e^{-zeta r}
-#define PTR_F12_ZETA    9
+#define PTR_RANGE_OMEGA         8
+// Yukawa potential and Slater-type geminal e^{-zeta r}
+#define PTR_F12_ZETA            9
 // Gaussian type geminal e^{-zeta r^2}
-#define PTR_GTG_ZETA    10
-#define NGRIDS          11
-#define PTR_GRIDS       12
-#define PTR_ENV_START   20
+#define PTR_GTG_ZETA            10
+#define NGRIDS                  11
+#define PTR_GRIDS               12
+#define PTR_ENV_START           20
+
 
 // slots of atm
 #define CHARGE_OF       0
@@ -113,8 +119,15 @@
 #define TSRZY       7
 #define TSRZZ       8
 
-// some other boundaries
-#define ANG_MAX         12 // l = 0..11
+// other boundaries
+#define MXRYSROOTS      32 // > ANG_MAX*2+1 for 4c2e
+#define ANG_MAX         15 // l = 0..15
+#define LMAX1           16 // > ANG_MAX
+#define CART_MAX        136 // > (ANG_MAX*(ANG_MAX+1)/2)
+#define SHLS_MAX        1048576
+#define NPRIM_MAX       64
+#define NCTR_MAX        64
+
 #define POINT_NUC       1
 #define GAUSSIAN_NUC    2
 #define FRAC_CHARGE_NUC 3
@@ -168,32 +181,29 @@ double *CINTc2s_ket_sph1(double *sph, double *cart, FINT lds, FINT ldc, FINT l);
 double CINTgto_norm(FINT n, double a);
 
 
-void CINTinit_2e_optimizer(CINTOpt **opt, const FINT *atm, const FINT natm,
-                           const FINT *bas, const FINT nbas, const double *env);
-void CINTinit_optimizer(CINTOpt **opt, const FINT *atm, const FINT natm,
-                        const FINT *bas, const FINT nbas, const double *env);
+void CINTinit_2e_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
+                           FINT *bas, FINT nbas, double *env);
+void CINTinit_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
+                        FINT *bas, FINT nbas, double *env);
 void CINTdel_2e_optimizer(CINTOpt **opt);
 void CINTdel_optimizer(CINTOpt **opt);
 
 
-CACHE_SIZE_T cint2e_cart(double *opijkl, const FINT *shls,
-                const FINT *atm, const FINT natm,
-                const FINT *bas, const FINT nbas, const double *env,
-                const CINTOpt *opt);
-void cint2e_cart_optimizer(CINTOpt **opt, const FINT *atm, const FINT natm,
-                           const FINT *bas, const FINT nbas, const double *env);
-CACHE_SIZE_T cint2e_sph(double *opijkl, const FINT *shls,
-               const FINT *atm, const FINT natm,
-               const FINT *bas, const FINT nbas, const double *env,
-               const CINTOpt *opt);
-void cint2e_sph_optimizer(CINTOpt **opt, const FINT *atm, const FINT natm,
-                          const FINT *bas, const FINT nbas, const double *env);
-CACHE_SIZE_T cint2e(double *opijkl, const FINT *shls,
-           const FINT *atm, const FINT natm,
-           const FINT *bas, const FINT nbas, const double *env,
-           const CINTOpt *opt);
-void cint2e_optimizer(CINTOpt **opt, const FINT *atm, const FINT natm,
-                      const FINT *bas, const FINT nbas, const double *env);
+FINT cint2e_cart(double *opijkl, FINT *shls,
+                FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env,
+                CINTOpt *opt);
+void cint2e_cart_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
+                           FINT *bas, FINT nbas, double *env);
+FINT cint2e_sph(double *opijkl, FINT *shls,
+               FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env,
+               CINTOpt *opt);
+void cint2e_sph_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
+                          FINT *bas, FINT nbas, double *env);
+FINT cint2e(double *opijkl, FINT *shls,
+           FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env,
+           CINTOpt *opt);
+void cint2e_optimizer(CINTOpt **opt, FINT *atm, FINT natm,
+                      FINT *bas, FINT nbas, double *env);
 
 #ifndef __cplusplus
 #include <complex.h>

include/cint_funcs.h

@@ -1072,3 +1072,15 @@ extern CINTIntegralFunction int1e_grids_ip_cart;
 extern CINTIntegralFunction int1e_grids_ip_sph;
 extern CINTIntegralFunction int1e_grids_ip_spinor;
 
+/* <NABLA i| 1/r_{grids} |NABLA j> */
+extern CINTOptimizerFunction int1e_grids_ipvip_optimizer;
+extern CINTIntegralFunction int1e_grids_ipvip_cart;
+extern CINTIntegralFunction int1e_grids_ipvip_sph;
+extern CINTIntegralFunction int1e_grids_ipvip_spinor;
+
+/* <SIGMA DOT P i| 1/r_{grids} |SIGMA DOT P j> */
+extern CINTOptimizerFunction int1e_grids_spvsp_optimizer;
+extern CINTIntegralFunction int1e_grids_spvsp_cart;
+extern CINTIntegralFunction int1e_grids_spvsp_sph;
+extern CINTIntegralFunction int1e_grids_spvsp_spinor;
+

scripts/auto_intor.cl

@@ -4,8 +4,8 @@
 (load "gen-code.cl")
 
 (gen-cint "intor1.c"
-  '("int1e_ovlp"                ( \| ))
-  '("int1e_nuc"                 ( \| nuc \| ))
+  ;'("int1e_ovlp"                ( \| ))
+  ;'("int1e_nuc"                 ( \| nuc \| ))
   '("int1e_kin"                 (.5 \| p dot p))
   '("int1e_ia01p"               (#C(0 1) \| nabla-rinv \| cross p))
   '("int1e_giao_irjxp"          (#C(0 1) \| r cross p))
@@ -222,5 +222,7 @@
 )
 
 (gen-cint "int1e_grids1.c"
-  '("int1e_grids_ip"           ( nabla \| grids \| ))
+  '("int1e_grids_ip"            ( nabla \| grids \| ))
+  '("int1e_grids_ipvip"         ( nabla \| grids \| nabla ))
+  '("int1e_grids_spvsp"         ( sigma dot p \| grids \| sigma dot p ))
 )

scripts/cart2sph.py

@@ -265,13 +265,29 @@ def sph2spinor(l):
     ua, ub = sph2spinor(l)
     for k in range(l * 2):
         print(f'// j = {l*2-1}/2, mj = {k*2+1-l*2}/2')
-        ca = [f'{mpmath.nstr(c.real, 18)} + {mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ua[:,k]]
-        cb = [f'{mpmath.nstr(c.real, 18)} + {mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ub[:,k]]
+        ca = [f'{mpmath.nstr(c.real, 18)}' for c in ua[:,k]]
+        cb = [f'{mpmath.nstr(c.real, 18)}' for c in ub[:,k]]
         print(f'{", ".join(ca)},')
         print(f'{", ".join(cb)},')
     for k in range(l * 2, l * 4 + 2):
         print(f'// j = {l*2+1}/2, mj = {k*2-1-l*6}/2')
-        ca = [f'{mpmath.nstr(c.real, 18)} + {mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ua[:,k]]
-        cb = [f'{mpmath.nstr(c.real, 18)} + {mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ub[:,k]]
+        ca = [f'{mpmath.nstr(c.real, 18)}' for c in ua[:,k]]
+        cb = [f'{mpmath.nstr(c.real, 18)}' for c in ub[:,k]]
+        print(f'{", ".join(ca)},')
+        print(f'{", ".join(cb)},')
+
+    l = 4
+    ncart = (l + 1) * (l + 2) // 2
+    ua, ub = sph2spinor(l)
+    for k in range(l * 2):
+        print(f'// j = {l*2-1}/2, mj = {k*2+1-l*2}/2')
+        ca = [f'{mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ua[:,k]]
+        cb = [f'{mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ub[:,k]]
+        print(f'{", ".join(ca)},')
+        print(f'{", ".join(cb)},')
+    for k in range(l * 2, l * 4 + 2):
+        print(f'// j = {l*2+1}/2, mj = {k*2-1-l*6}/2')
+        ca = [f'{mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ua[:,k]]
+        cb = [f'{mpmath.nstr(c.imag, 18)}*_Complex_I' for c in ub[:,k]]
         print(f'{", ".join(ca)},')
         print(f'{", ".join(cb)},')

scripts/gen-code.cl

@@ -386,7 +386,44 @@ ix = idx[0+n*3];
 iy = idx[1+n*3];
 iz = idx[2+n*3];~%")
       (dump-s-for-nroots fout tot-bits 1)
-      (gen-c-block+ fout flat-script)
+      (gen-c-block-with-empty fout flat-script)
+      (format fout "}}~%")
+      goutinc)))
+
+(defun gen-code-gout1e-rinv (fout intname raw-infix flat-script)
+  (destructuring-bind (op bra-i ket-j bra-k ket-l)
+    (split-int-expression raw-infix)
+    (let* ((i-rev (effect-keys bra-i)) ;<i| already in reverse order
+           (j-rev (reverse (effect-keys ket-j)))
+           (op-rev (reverse (effect-keys op)))
+           (i-len (length i-rev))
+           (j-len (length j-rev))
+           (op-len (length op-rev))
+           (tot-bits (+ i-len j-len op-len))
+           (goutinc (length flat-script)))
+      (format fout "void CINTgout1e_~a" intname)
+      (format fout "(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
+FINT nf = envs->nf;
+FINT nrys_roots = envs->nrys_roots;
+FINT ix, iy, iz, n, i;
+double *g0 = g;~%")
+      (loop
+        for i in (range (num-g-intermediates tot-bits op i-len j-len)) do
+        (format fout "double *g~a = g~a + envs->g_size * 3;~%" (1+ i) i))
+      (dump-declare-dri-for-rc fout bra-i "i")
+      (dump-declare-dri-for-rc fout (append op ket-j) "j")
+      (dump-declare-giao-ij fout bra-i (append op ket-j))
+;;; generate g_(bin)
+;;; for the operators act on the |ket>, the reversed scan order and r_combinator
+;;; is required; for the operators acto on the <bra|, the normal scan order and
+      (let ((fmt-i "G2E_~aI(g~a, g~a, envs->i_l+~a, envs->j_l, 0, 0);~%")
+            (fmt-op (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0, 0);
+G2E_~aI(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0, 0);
+for (ix = 0; ix < envs->g_size * 3; ix++) {g~a[ix] += g~a[ix];}~%"))
+            (fmt-j (mkstr "G2E_~aJ(g~a, g~a, envs->i_l+~d, envs->j_l+~a, 0, 0);~%")))
+        (dump-combo-braket fout fmt-i fmt-op fmt-j i-rev op-rev j-rev 0))
+      (dump-s-2e fout tot-bits 0)
+      (gen-c-block-with-empty fout flat-script)
       (format fout "}}~%")
       goutinc)))
 
@@ -415,8 +452,11 @@ iz = idx[2+n*3];~%")
                     (if (or (member 'nuc raw-infix)
                             (member 'rinv raw-infix)
                             (member 'nabla-rinv raw-infix))
-                      (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%"
-                              i-len (+ op-len j-len) tot-bits e1comps tensors)
+                      (if (intersection *act-left-right* op)
+                          (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%"
+                                  (1+ i-len) (+ op-len j-len) tot-bits e1comps tensors)
+                          (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%"
+                                  i-len (+ op-len j-len) tot-bits e1comps tensors))
                       (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 1, ~d};~%"
                               i-len (+ op-len j-len) tot-bits e1comps tensors))))
            (envs-common (with-output-to-string (tmpout)
@@ -429,7 +469,11 @@ iz = idx[2+n*3];~%")
       (write-functype-header
         t intname (format nil "/* <~{~a ~}i|~{~a ~}|~{~a ~}j> */" bra-i op ket-j))
       (format fout "/* <~{~a ~}i|~{~a ~}|~{~a ~}j> */~%" bra-i op ket-j)
-      (gen-code-gout1e fout intname raw-infix flat-script)
+      (cond ((or (member 'nuc raw-infix)
+                 (member 'rinv raw-infix)
+                 (member 'nabla-rinv raw-infix))
+             (gen-code-gout1e-rinv fout intname raw-infix flat-script))
+            (t (gen-code-gout1e fout intname raw-infix flat-script)))
       (format fout "void ~a_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {~%" intname)
       (format fout ngdef)
       (format fout "CINTall_1e_optimizer(opt, ng, atm, natm, bas, nbas, env);~%}~%")
@@ -1290,8 +1334,11 @@ for (ig = 0; ig < bgrids; ig++) {~%")
            (e1comps (if (eql sf 'sf) 1 4))
            (tensors (if (eql sf 'sf) goutinc (/ goutinc 4)))
            (ngdef (with-output-to-string (tmpout)
-                    (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%"
-                            i-len (+ op-len j-len) tot-bits e1comps tensors)))
+                    (if (intersection *act-left-right* op)
+                        (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%"
+                                (1+ i-len) (+ op-len j-len) tot-bits e1comps tensors)
+                        (format tmpout "FINT ng[] = {~d, ~d, 0, 0, ~d, ~d, 0, ~d};~%"
+                                i-len (+ op-len j-len) tot-bits e1comps tensors))))
            (envs-common (with-output-to-string (tmpout)
                           (format tmpout ngdef)
                           (format tmpout "CINTEnvVars envs;~%")