Version 3.1.0

Add p orbital convention py, pz, px
sunqm · Aug 25, 2020 · 1e3d00a · 1e3d00a
1 parent 97dcffd
commit 1e3d00a
Show file tree

Hide file tree

Showing 7 changed files with 59 additions and 26 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required (VERSION 2.6)
 project (qcint C)
-set(qcint_VERSION "3.0.19")
+set(qcint_VERSION "3.1.0")
 set(qcint_SOVERSION "3")
 
 if ("${CMAKE_BUILD_TYPE}" STREQUAL "")
@@ -89,6 +89,11 @@ if(WITH_GTG)
   message("Enabled WITH_GTG")
 endif(WITH_GTG)
 
+if(PYPZPX)
+  add_definitions(-DPYPZPX)
+  message("P orbitals are sorted to (py, pz, px)")
+endif(PYPZPX)
+
 option(BUILD_SHARED_LIBS "build shared libraries" 1)
 option(ENABLE_STATIC "Enforce static library build" 0)
 if(ENABLE_STATIC)

diff --git a/ChangeLog b/ChangeLog
@@ -1,3 +1,5 @@
+Version 3.1.0  (2020-08-25):
+	* Add convention py, pz, px
 Version 3.0.19  (2019-11-30):
 	* support fractional nuclear charges
 Version 3.0.18  (2019-07-02):

diff --git a/README.md b/README.md
@@ -3,9 +3,9 @@ qcint (quick libcint)
 
 An optimized libcint branch for X86 platform
 
-version 3.0.19
+version 3.1.0
 
-2019-11-30
+2020-08-25
 
 
 What is qcint

diff --git a/src/cart2sph.c b/src/cart2sph.c
@@ -40,6 +40,22 @@
 static const double g_trans_cart2sph[] = {
         1, /* factors of s and p are moved to CINTcommon_fac_sp */
         // px
+#ifdef PYPZPX
+        // py
+        0,
+        1,
+        0,
+        // pz
+        0,
+        0,
+        1,
+        // px
+        1,
+        0,
+        0,
+#else
+// by default, p orbitals are ordered px, py, pz
+        // px
         1,
         0,
         0,
@@ -51,6 +67,7 @@ static const double g_trans_cart2sph[] = {
         0,
         0,
         1,
+#endif
         // dxy
         0,
         1.092548430592079070,
@@ -2913,40 +2930,50 @@ static double *s_ket_cart2spheric1(double *gsph, double *gcart,
 // transform p function from cartesian to spheric
 static double *p_bra_cart2spheric(double *gsph, int nket, double *gcart, int l)
 {
-        /*
-        double *pgcart = gcart;
+#ifdef PYPZPX
         int i;
         for (i = 0; i < nket; i++) {
-                gsph[0] = gcart[0];
-                gsph[1] = gcart[1];
-                gsph[2] = gcart[2];
-                gsph += 3;
-                gcart += 3;
+                gsph[i*nket+0] = gcart[i*nket+1];  // py
+                gsph[i*nket+1] = gcart[i*nket+2];  // pz
+                gsph[i*nket+2] = gcart[i*nket+0];  // px
         }
-        return pgcart;*/
+        return gsph;
+#else
         return gcart;
+#endif
 }
 static double *p_ket_cart2spheric(double *gsph, double *gcart,
                                   int lds, int nbra, int l)
 {
-        /*
+#ifdef PYPZPX
         int i;
         for (i = 0; i < nbra; i++) {
-                gsph[0*nbra+i] = gcart[0*nbra+i];
-                gsph[1*nbra+i] = gcart[1*nbra+i];
-                gsph[2*nbra+i] = gcart[2*nbra+i];
-        }*/
+                gsph[0*nbra+i] = gcart[1*nbra+i];  // py
+                gsph[1*nbra+i] = gcart[2*nbra+i];  // pz
+                gsph[2*nbra+i] = gcart[0*nbra+i];  // px
+        }
+        return gsph;
+#else
         return gcart;
+#endif
 }
 static double *p_ket_cart2spheric1(double *gsph, double *gcart,
                                    int lds, int nbra, int l)
 {
         int i;
+#ifdef PYPZPX
+        for (i = 0; i < nbra; i++) {
+                gsph[0*nbra+i] = gcart[1*nbra+i];  // py
+                gsph[1*nbra+i] = gcart[2*nbra+i];  // pz
+                gsph[2*nbra+i] = gcart[0*nbra+i];  // px
+        }
+#else
         for (i = 0; i < nbra; i++) {
                 gsph[0*lds+i] = gcart[0*nbra+i];
                 gsph[1*lds+i] = gcart[1*nbra+i];
                 gsph[2*lds+i] = gcart[2*nbra+i];
         }
+#endif
         return gsph;
 }
 

diff --git a/src/cint2e.c b/src/cint2e.c
@@ -160,13 +160,13 @@
 //    double d;
 //    unsigned short s[4];
 //} type_IEEE754;
-static double approx_log(double x)
+static inline double approx_log(double x)
 {
         //type_IEEE754 y;
         //y.d = x;
         //return ((double)(y.s[3] >> 4) - 1023) * 0.7;
         //return log(x);
-        return 2.5;
+        return 3;
 }
 
 int CINT2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache)
@@ -244,8 +244,7 @@ int CINT2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache)
         ALIGNMM Rys2eT bc;
         double rr_kl = SQUARE(envs->rkrl);
         double log_rr_kl = (envs->lk_ceil+envs->ll_ceil+1)*approx_log(rr_kl+1)/2;
-        double akl, ekl, expijkl;
-        int ccekl;
+        double akl, ekl, expijkl, ccekl, eijcutoff;
         ALIGNMM double rkl[4];
         int non0ctri[i_prim];
         int non0ctrj[j_prim];
@@ -284,10 +283,11 @@ int CINT2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache)
                 for (kp = 0; kp < k_prim; kp++) {
                         akl = 1 / (ak[kp] + al[lp]);
                         ekl = rr_kl * ak[kp] * al[lp] * akl;
-                        ccekl = (int)(ekl - log_rr_kl);
+                        ccekl = ekl - log_rr_kl;
                         if (ccekl > CUTOFF15) {
                                 goto k_contracted;
                         }
+                        eijcutoff = CUTOFF15 - ccekl;
                         rkl[0] = (ak[kp]*rk[0] + al[lp]*rl[0]) * akl;
                         rkl[1] = (ak[kp]*rk[1] + al[lp]*rl[1]) * akl;
                         rkl[2] = (ak[kp]*rk[2] + al[lp]*rl[2]) * akl;
@@ -298,8 +298,7 @@ int CINT2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache)
                         for (jp = 0; jp < j_prim; jp++) {
                                 INIT_GCTR_ADDR(i, j, fac1k);
                                 for (ip = 0; ip < i_prim; ip++, pdata_ij++) {
-                                        if (pdata_ij->cceij > CUTOFF15 ||
-                                            pdata_ij->cceij+ccekl > CUTOFF15) {
+                                        if (pdata_ij->cceij > eijcutoff) {
                                                 goto i_contracted;
                                         }
                                         expijkl = pdata_ij->eij * ekl;

diff --git a/src/optimizer.c b/src/optimizer.c
@@ -459,13 +459,13 @@ void CINTOpt_set_non0coeff(CINTOpt *opt, int *atm, int natm,
 //    double d;
 //    unsigned short s[4];
 //} type_IEEE754;
-static double approx_log(double x)
+static inline double approx_log(double x)
 {
         //type_IEEE754 y;
         //y.d = x;
         //return ((double)(y.s[3] >> 4) - 1023) * 0.7;
         //return log(x);
-        return 2.5;
+        return 3;
 }
 
 int CINTset_pairdata(PairData *pdata, double *ai, double *aj, double *ri, double *rj,

diff --git a/src/rys_roots.c b/src/rys_roots.c
@@ -1739,7 +1739,7 @@ static void R_qnode(__float128 *a, __float128 *rt, int order)
                         xi = x0 + (x0 - x1) / (p1 - p0) * p0;
                 }
                 n = 0;
-                while (x1 > accrt+x0 || x0 > x1+accrt) {
+                while (fabs(x1 - x0) > accrt) {
                         n++;
                         if (n > 600) {
                                 fprintf(stderr, "libcint::rys_roots NO CONV. IN R_qnode\n");