Add Gaussian type geminals

CMakeLists.txt

@@ -82,6 +82,12 @@ if(WITH_F12)
   message("Enabled WITH_F12")
 endif(WITH_F12)
 
+if(WITH_GTG)
+  set(cintSrc ${cintSrc} src/g2e_gtg.c src/cint2e_gtg.c)
+  add_definitions(-DWITH_GTG)
+  message("Enabled WITH_GTG")
+endif(WITH_GTG)
+
 option(BUILD_SHARED_LIBS "build shared libraries" 1)
 option(ENABLE_STATIC "Enforce static library build" 0)
 if(ENABLE_STATIC)

include/cint.h.in

@@ -30,6 +30,7 @@
 #define PTR_RANGE_OMEGA         8
 // Yukawa potential and slater-type geminal e^{-zeta r}
 #define PTR_F12_ZETA            9
+#define PTR_GTG_ZETA            10
 #define PTR_ENV_START           20
 
 // slots of atm

src/cint2e_f12.c

@@ -69,8 +69,8 @@ void c##NAME##_sph_optimizer(CINTOpt **opt, int *atm, int natm, \
         NAME##_optimizer(opt, atm, natm, bas, nbas, env); \
 }
 
-ALL_CINT(int2e_yp)
-ALL_CINT(int2e_stg)
+ALL_CINT(int2e_yp);
+ALL_CINT(int2e_stg);
 
 
 

src/cint2e_gtg.c

@@ -0,0 +1,56 @@
+/*
+ * Qcint is a general GTO integral library for computational chemistry
+ * Copyright (C) 2014- Qiming Sun <[email protected]>
+ *
+ * This file is part of Qcint.
+ *
+ * Qcint is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+#include <stdlib.h>
+#include "cint_bas.h"
+#include "g2e.h"
+#include "optimizer.h"
+#include "cint2e.h"
+
+int int2e_gtg_sph(double *out, int *dims, int *shls, int *atm, int natm,
+                  int *bas, int nbas, double *env, CINTOpt *opt, double *cache)
+{
+        int ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
+        CINTEnvVars envs;
+        CINTinit_int2e_gtg_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
+        envs.f_gout = &CINTgout2e;
+        envs.f_gout_simd1 = &CINTgout2e_simd1;
+        return CINT2e_spheric_drv(out, dims, &envs, opt, cache);
+}
+void int2e_gtg_optimizer(CINTOpt **opt, int *atm, int natm,
+                         int *bas, int nbas, double *env)
+{
+        int ng[] = {0, 0, 0, 0, 0, 1, 1, 1};
+        CINTall_2e_gtg_optimizer(opt, ng, atm, natm, bas, nbas, env);
+}
+
+#define ALL_CINT(NAME) \
+int c##NAME##_sph(double *out, int *shls, int *atm, int natm, \
+            int *bas, int nbas, double *env, CINTOpt *opt) { \
+        return NAME##_sph(out, NULL, shls, atm, natm, bas, nbas, env, opt, NULL); \
+} \
+void c##NAME##_sph_optimizer(CINTOpt **opt, int *atm, int natm, \
+                         int *bas, int nbas, double *env) { \
+        NAME##_optimizer(opt, atm, natm, bas, nbas, env); \
+}
+
+ALL_CINT(int2e_gtg);
+

src/cint_const.h

@@ -28,6 +28,7 @@
 #define PTR_RANGE_OMEGA         8
 // Yukawa potential and slater-type geminal e^{-zeta r}
 #define PTR_F12_ZETA            9
+#define PTR_GTG_ZETA            10
 #define PTR_ENV_START           20
 
 // slots of atm

src/g2e.c

@@ -439,7 +439,7 @@ MM_STORE(p0z+j*SIMDD, MM_LOAD(c00z+i*SIMDD) * MM_LOAD(gz+j*SIMDD) + r0 * MM_LOAD
  * g0[i,k,l,j] = < ik | lj > = ( i j | k l )
  */
 /* 2d is based on l,j */
-void CINTg0_lj2d_4d(double *g, CINTEnvVars *envs)
+void CINTg0_lj_4d(double *g, CINTEnvVars *envs)
 {
         int li = envs->li_ceil;
         int lk = envs->lk_ceil;
@@ -1817,7 +1817,7 @@ void CINTg0_2e_lj2d4d(double *g, Rys2eT *bc, CINTEnvVars *envs)
         default:
 _g0_4d_default:
                         CINTg0_2e_2d(g, bc, envs);
-                        CINTg0_lj2d_4d(g, envs);
+                        CINTg0_lj_4d(g, envs);
                         return;
         }
 error:
@@ -1842,15 +1842,6 @@ void CINTg0_2e_il2d4d(double *g, Rys2eT *bc, CINTEnvVars *envs)
         CINTg0_il_4d(g, envs);
 }
 
-#ifdef WITH_F12
-void CINTg0_2e_stg_lj2d4d(double *g, Rys2eT *bc, CINTEnvVars *envs)
-{
-        CINTg0_2e_2d(g, bc, envs);
-        CINTg0_lj2d_4d(g, envs);
-}
-#endif
-
-
 /*
  * g[i,k,l,j] = < ik | lj > = ( i j | k l )
  */
@@ -2039,13 +2030,12 @@ void CINTg0_2e(double *g, Rys2eT *bc, CINTEnvVars *envs, int count)
         r0 = MM_LOAD(a0);
         r1 = MM_LOAD(a1);
         r2 = MM_SET1(.5);
-        r3 = MM_SET1(2.);
+        r3 = MM_SET1(1.);
         for (i = 0; i < nroots; i++) {
                 r4 = MM_MUL(r0, MM_LOAD(u+i*SIMDD));
-                r5 = MM_DIV(r2, MM_FMA(r4, ra, r1));
-                MM_STORE(tmp4+i*SIMDD, r5);
-                r6 = MM_MUL(MM_MUL(r3, r4), r5);
-                MM_STORE(tmp1+i*SIMDD, r6);
+                r5 = MM_DIV(r3, MM_FMA(r4, ra, r1));
+                MM_STORE(tmp4+i*SIMDD, MM_MUL(r2, r5));
+                MM_STORE(tmp1+i*SIMDD, MM_MUL(r4, r5));
         }
         ra = MM_SET1(.5);
         r2 = MM_LOAD(akl);

src/g2e.h

@@ -109,6 +109,18 @@ void CINTx1k_2e_simd1(double *f, double *g, double *rk,
 void CINTx1l_2e_simd1(double *f, double *g, double *rl,
                       int li, int lj, int lk, int ll, CINTEnvVars *envs);
 
+#ifdef WITH_F12
+void CINTinit_int2e_yp_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
+                               int *atm, int natm, int *bas, int nbas, double *env);
+void CINTinit_int2e_stg_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
+                                int *atm, int natm, int *bas, int nbas, double *env);
+#endif
+
+#ifdef WITH_GTG
+void CINTinit_int2e_gtg_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
+                                int *atm, int natm, int *bas, int nbas, double *env);
+#endif
+
 
 #define G2E_D_I(f, g, li, lj, lk, ll)   CINTnabla1i_2e(f, g, li, lj, lk, ll, envs)
 #define G2E_D_J(f, g, li, lj, lk, ll)   CINTnabla1j_2e(f, g, li, lj, lk, ll, envs)

src/g2e_f12.c

@@ -33,8 +33,10 @@ void CINTg0_2e_yp(double *g, Rys2eT *bc, CINTEnvVars *envs, int count);
 void CINTg0_2e_yp_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs, int idsimd);
 void CINTg0_2e_stg(double *g, Rys2eT *bc, CINTEnvVars *envs, int count);
 void CINTg0_2e_stg_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs, int idsimd);
-void CINTg0_2e_stg_lj2d4d(double *g, Rys2eT *bc, const CINTEnvVars *envs);
-void CINTg0_2e_stg_lj2d4d_simd1(double *g, Rys2eT *bc, const CINTEnvVars *envs);
+void CINTg0_2e_stg_lj2d4d(double *g, Rys2eT *bc, CINTEnvVars *envs);
+void CINTg0_2e_stg_lj2d4d_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs);
+void CINTg0_lj_4d(double *g, CINTEnvVars *envs);
+void CINTg0_lj_4d_simd1(double *g, CINTEnvVars *envs);
 
 void CINTinit_int2e_yp_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
                                int *atm, int natm, int *bas, int nbas, double *env)
@@ -167,6 +169,18 @@ void CINTinit_int2e_stg_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
         envs->f_g0_2e_simd1 = &CINTg0_2e_stg_simd1;
 }
 
+void CINTg0_2e_stg_lj2d4d(double *g, Rys2eT *bc, CINTEnvVars *envs)
+{
+        CINTg0_2e_2d(g, bc, envs);
+        CINTg0_lj_4d(g, envs);
+}
+
+void CINTg0_2e_stg_lj2d4d_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs)
+{
+        CINTg0_2e_2d_simd1(g, bc, envs);
+        CINTg0_lj_4d_simd1(g, envs);
+}
+
 void CINTg0_2e_yp(double *g, Rys2eT *bc, CINTEnvVars *envs, int count)
 {
         ALIGNMM double aij[SIMDD];
@@ -273,13 +287,12 @@ void CINTg0_2e_yp(double *g, Rys2eT *bc, CINTEnvVars *envs, int count)
         r0 = MM_LOAD(a0);
         r1 = MM_LOAD(a1);
         r2 = MM_SET1(.5);
-        r3 = MM_SET1(2.);
+        r3 = MM_SET1(1.);
         for (i = 0; i < nroots; i++) {
                 r4 = MM_MUL(r0, MM_LOAD(u+i*SIMDD));
-                r5 = MM_DIV(r2, MM_FMA(r4, ra, r1));
-                MM_STORE(tmp4+i*SIMDD, r5);
-                r6 = MM_MUL(MM_MUL(r3, r4), r5);
-                MM_STORE(tmp1+i*SIMDD, r6);
+                r5 = MM_DIV(r3, MM_FMA(r4, ra, r1));
+                MM_STORE(tmp4+i*SIMDD, MM_MUL(r2, r5));
+                MM_STORE(tmp1+i*SIMDD, MM_MUL(r4, r5));
         }
         ra = MM_SET1(.5);
         r2 = MM_LOAD(akl);
@@ -320,7 +333,11 @@ void CINTg0_2e_yp(double *g, Rys2eT *bc, CINTEnvVars *envs, int count)
 
 void CINTg0_2e_yp_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs, int idsimd)
 {
-        double aij, akl, a0, a1, fac1;
+        const double aij = envs->ai[idsimd] + envs->aj[idsimd];
+        const double akl = envs->ak[idsimd] + envs->al[idsimd];
+        const double zeta = envs->env[PTR_F12_ZETA];
+        const int nroots = envs->nrys_roots;
+        double a0, a1, fac1;
         ALIGNMM double x[SIMDD];
         ALIGNMM double ua[SIMDD];
         double *rij = envs->rij;
@@ -330,12 +347,8 @@ void CINTg0_2e_yp_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs, int idsimd)
         double rklrx[3];
         double *u = bc->u;
         double *w = bc->w;
-        double zeta = envs->env[PTR_F12_ZETA];
-        int nroots = envs->nrys_roots;
         int i;
 
-        aij = envs->ai[idsimd] + envs->aj[idsimd];
-        akl = envs->ak[idsimd] + envs->al[idsimd];
         a1 = aij * akl;
         a0 = a1 / (aij + akl);
         //fac1 = sqrt(a0 / (a1 * a1 * a1)) * envs->fac[idsimd];
@@ -527,13 +540,12 @@ void CINTg0_2e_stg(double *g, Rys2eT *bc, CINTEnvVars *envs, int count)
         r0 = MM_LOAD(a0);
         r1 = MM_LOAD(a1);
         r2 = MM_SET1(.5);
-        r3 = MM_SET1(2.);
+        r3 = MM_SET1(1.);
         for (i = 0; i < nroots; i++) {
                 r4 = MM_MUL(r0, MM_LOAD(u+i*SIMDD));
-                r5 = MM_DIV(r2, MM_FMA(r4, ra, r1));
-                MM_STORE(tmp4+i*SIMDD, r5);
-                r6 = MM_MUL(MM_MUL(r3, r4), r5);
-                MM_STORE(tmp1+i*SIMDD, r6);
+                r5 = MM_DIV(r3, MM_FMA(r4, ra, r1));
+                MM_STORE(tmp4+i*SIMDD, MM_MUL(r2, r5));
+                MM_STORE(tmp1+i*SIMDD, MM_MUL(r4, r5));
         }
         ra = MM_SET1(.5);
         r2 = MM_LOAD(akl);
@@ -574,7 +586,11 @@ void CINTg0_2e_stg(double *g, Rys2eT *bc, CINTEnvVars *envs, int count)
 
 void CINTg0_2e_stg_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs, int idsimd)
 {
-        double aij, akl, a0, a1, fac1;
+        const double aij = envs->ai[idsimd] + envs->aj[idsimd];
+        const double akl = envs->ak[idsimd] + envs->al[idsimd];
+        const double zeta = envs->env[PTR_F12_ZETA];
+        const int nroots = envs->nrys_roots;
+        double a0, a1, fac1;
         ALIGNMM double x[SIMDD];
         ALIGNMM double ua[SIMDD];
         double *rij = envs->rij;
@@ -584,12 +600,8 @@ void CINTg0_2e_stg_simd1(double *g, Rys2eT *bc, CINTEnvVars *envs, int idsimd)
         double rklrx[3];
         double *u = bc->u;
         double *w = bc->w;
-        double zeta = envs->env[PTR_F12_ZETA];
-        int nroots = envs->nrys_roots;
         int i;
 
-        aij = envs->ai[idsimd] + envs->aj[idsimd];
-        akl = envs->ak[idsimd] + envs->al[idsimd];
         a1 = aij * akl;
         a0 = a1 / (aij + akl);
         //fac1 = sqrt(a0 / (a1 * a1 * a1)) * envs->fac[idsimd];