diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2698818a..c3b93945 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -18,7 +18,7 @@ jobs: - name: Compile run: | env - cmake -DWITH_RANGE_COULOMB=1 -DWITH_COULOMB_ERF=1 -DWITH_F12=1 -Bbuild -DKEEP_GOING=1 . + cmake -DWITH_CINT2_INTERFACE=1 -DWITH_RANGE_COULOMB=1 -DWITH_COULOMB_ERF=1 -DWITH_F12=1 -DWITH_4C1E=1 -Bbuild -DKEEP_GOING=1 . cmake --build build - name: Test for rys-roots if: startsWith(matrix.os, 'ubuntu') @@ -34,3 +34,4 @@ jobs: python test_3c2e.py --quick python test_cart2sph.py python test_c2s.py + python test_cint4c1e.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 79eb793d..cd95018f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required (VERSION 3.5) project (cint C) set(cint_VERSION_MAJOR "5") set(cint_VERSION_MINOR "1") -set(cint_VERSION_PATCH "5") +set(cint_VERSION_PATCH "6") set(cint_VERSION_TWEAK "0") set(cint_VERSION "${cint_VERSION_MAJOR}.${cint_VERSION_MINOR}.${cint_VERSION_PATCH}") set(cint_SOVERSION "${cint_VERSION_MAJOR}") diff --git a/ChangeLog b/ChangeLog index 96a3a831..d6690e78 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,6 @@ +Version 5.1.6 (2022-08-26): + * Add 4-center 1-electron intergral int4c1e + Version 5.1.5 (2022-08-24): * New integral int3c1e_ip1 diff --git a/README.rst b/README.rst index 47248e7c..f84859a3 100644 --- a/README.rst +++ b/README.rst @@ -1,8 +1,8 @@ libcint ======= -version 5.1.5 -2022-08-24 +version 5.1.6 +2022-08-26 What is libcint diff --git a/src/cint3c1e.c b/src/cint3c1e.c index 4084e8e1..1a1f5a2d 100644 --- a/src/cint3c1e.c +++ b/src/cint3c1e.c @@ -454,26 +454,7 @@ CACHE_SIZE_T CINT3c1e_spinor_drv(double complex *out, FINT *dims, CINTEnvVars *e exit(1); } -void CINTgout3c1e(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) -{ - FINT ix, iy, iz, n; - - if (gout_empty) { - for (n = 0; n < envs->nf; n++, idx+=3) { - ix = idx[0]; - iy = idx[1]; - iz = idx[2]; - gout[n] = g[ix] * g[iy] * g[iz]; - } - } else { - for (n = 0; n < envs->nf; n++, idx+=3) { - ix = idx[0]; - iy = idx[1]; - iz = idx[2]; - gout[n] += g[ix] * g[iy] * g[iz]; - } - } -} +void CINTgout1e(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT empty); CACHE_SIZE_T int3c1e_sph(double *out, FINT *dims, FINT *shls, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) @@ -481,7 +462,7 @@ CACHE_SIZE_T int3c1e_sph(double *out, FINT *dims, FINT *shls, FINT *atm, FINT na FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, 0, 0); } void int3c1e_optimizer(CINTOpt **opt, FINT *atm, FINT natm, @@ -496,7 +477,7 @@ CACHE_SIZE_T int3c1e_cart(double *out, FINT *dims, FINT *shls, FINT *atm, FINT n FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, 0, 0); } @@ -506,7 +487,7 @@ CACHE_SIZE_T int3c1e_spinor(double complex *out, FINT *dims, FINT *shls, FINT *a FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0, 0); } @@ -516,7 +497,7 @@ CACHE_SIZE_T int3c1e_rinv_sph(double *out, FINT *dims, FINT *shls, FINT *atm, FI FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, 1, 0); } void int3c1e_rinv_optimizer(CINTOpt **opt, FINT *atm, FINT natm, @@ -531,7 +512,7 @@ CACHE_SIZE_T int3c1e_rinv_cart(double *out, FINT *dims, FINT *shls, FINT *atm, F FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, 1, 0); } @@ -541,7 +522,7 @@ CACHE_SIZE_T int3c1e_rinv_spinor(double complex *out, FINT *dims, FINT *shls, FI FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 1, 0); } diff --git a/src/cint4c1e.c b/src/cint4c1e.c index 39ca5edb..1fa81d3b 100644 --- a/src/cint4c1e.c +++ b/src/cint4c1e.c @@ -31,13 +31,12 @@ } \ *ctrsymb##empty = 0 -void CINTg4c1e_ovlp(double *g, double fac, CINTEnvVars *envs); +void CINTg4c1e_ovlp(double *g, CINTEnvVars *envs, double *cache); void CINTg4c1e_index_xyz(FINT *idx, CINTEnvVars *envs); CACHE_SIZE_T CINTinit_int4c1e_EnvVars(CINTEnvVars *envs, FINT *ng, FINT *shls, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env); -void CINTgout3c1e(double *g, double *gout, FINT *idx, - CINTEnvVars *envs, FINT gout_empty); +void CINTgout1e(double *gout, double *g, FINT *idx, CINTEnvVars *envs, FINT empty); CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) { @@ -146,8 +145,10 @@ CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) gout = g1; } - double eij, ekl, expijkl, eijkl, a0; + double eij, ekl, expijkl, eijkl, a0, aij, akl; double rijrkl[3]; + double rij[3]; + double rkl[3]; double dist_ij = SQUARE(envs->rirj); double dist_kl = SQUARE(envs->rkrl); double common_fac = envs->common_factor * SQRTPI * M_PI * @@ -156,7 +157,7 @@ CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) *lempty = 1; for (lp = 0; lp < l_prim; lp++) { - envs->al = al[lp]; + envs->al[0] = al[lp]; if (l_ctr == 1) { fac1l = common_fac * cl[lp]; } else { @@ -164,15 +165,18 @@ CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) *kempty = 1; } for (kp = 0; kp < k_prim; kp++) { - envs->ak = ak[kp]; - envs->akl = ak[kp] + al[lp]; - ekl = dist_kl * ak[kp] * al[lp] / envs->akl; + envs->ak[0] = ak[kp]; + akl = ak[kp] + al[lp]; + ekl = dist_kl * ak[kp] * al[lp] / akl; if (ekl > EXPCUTOFF) { goto k_contracted; } - envs->rkl[0] = (ak[kp]*rk[0] + al[lp]*rl[0]) / envs->akl; - envs->rkl[1] = (ak[kp]*rk[1] + al[lp]*rl[1]) / envs->akl; - envs->rkl[2] = (ak[kp]*rk[2] + al[lp]*rl[2]) / envs->akl; + 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; + envs->rkl[0] = rkl[0]; + envs->rkl[1] = rkl[1]; + envs->rkl[2] = rkl[2]; if (k_ctr == 1) { fac1k = fac1l * ck[kp]; } else { @@ -181,7 +185,7 @@ CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) } for (jp = 0; jp < j_prim; jp++) { - envs->aj = aj[jp]; + envs->aj[0] = aj[jp]; if (j_ctr == 1) { fac1j = fac1k * cj[jp]; } else { @@ -189,19 +193,22 @@ CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) *iempty = 1; } for (ip = 0; ip < i_prim; ip++) { - envs->ai = ai[ip]; - envs->aij = ai[ip] + aj[jp]; - eij = dist_ij * ai[ip] * aj[jp] / envs->aij; + envs->ai[0] = ai[ip]; + aij = ai[ip] + aj[jp]; + eij = dist_ij * ai[ip] * aj[jp] / aij; if (eij+ekl > EXPCUTOFF) { goto i_contracted; } - envs->rij[0] = (ai[ip]*ri[0] + aj[jp]*rj[0]) / envs->aij; - envs->rij[1] = (ai[ip]*ri[1] + aj[jp]*rj[1]) / envs->aij; - envs->rij[2] = (ai[ip]*ri[2] + aj[jp]*rj[2]) / envs->aij; + rij[0] = (ai[ip]*ri[0] + aj[jp]*rj[0]) / aij; + rij[1] = (ai[ip]*ri[1] + aj[jp]*rj[1]) / aij; + rij[2] = (ai[ip]*ri[2] + aj[jp]*rj[2]) / aij; + envs->rij[0] = rij[0]; + envs->rij[1] = rij[1]; + envs->rij[2] = rij[2]; rijrkl[0] = rij[0] - rkl[0]; rijrkl[1] = rij[1] - rkl[1]; rijrkl[2] = rij[2] - rkl[2]; - a0 = envs->aij * envs->akl / (envs->aij + envs->akl); + a0 = aij * akl / (aij + akl); eijkl = a0 * SQUARE(rijrkl); eijkl += eij + ekl; if (eijkl > EXPCUTOFF) { @@ -214,7 +221,7 @@ CACHE_SIZE_T CINT4c1e_loop_nopt(double *gctr, CINTEnvVars *envs, double *cache) fac1i = fac1j*expijkl; } envs->fac[0] = fac1i; - CINTg4c1e_ovlp(g, envs); + CINTg4c1e_ovlp(g, envs, cache); (*envs->f_gout)(gout, g, idx, envs, *gempty); PRIM2CTR0(i, gout, envs->nf*n_comp); i_contracted: ; @@ -284,7 +291,7 @@ CACHE_SIZE_T CINT4c1e_drv(double *out, FINT *dims, CINTEnvVars *envs, CINTOpt *o has_value = CINT4c1e_loop_nopt(gctr, envs, cache); FINT counts[4]; - if (f_c2s == &c2s_sph_1e) { + if (f_c2s == &c2s_sph_2e1) { counts[0] = (envs->i_l*2+1) * x_ctr[0]; counts[1] = (envs->j_l*2+1) * x_ctr[1]; counts[2] = (envs->k_l*2+1) * x_ctr[2]; @@ -320,7 +327,7 @@ CACHE_SIZE_T int4c1e_sph(double *out, FINT *dims, FINT *shls, FINT *atm, FINT na FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int4c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT4c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_2e1); } void int4c1e_optimizer(CINTOpt **opt, FINT *atm, FINT natm, @@ -335,7 +342,7 @@ CACHE_SIZE_T int4c1e_cart(double *out, FINT *dims, FINT *shls, FINT *atm, FINT n FINT ng[] = {0, 0, 0, 0, 0, 1, 1, 1}; CINTEnvVars envs; CINTinit_int4c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env); - envs.f_gout = &CINTgout3c1e; + envs.f_gout = &CINTgout1e; return CINT4c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_2e1); } diff --git a/src/g4c1e.c b/src/g4c1e.c index 6a0e29ec..054db8cc 100644 --- a/src/g4c1e.c +++ b/src/g4c1e.c @@ -141,7 +141,7 @@ void CINTg4c1e_index_xyz(FINT *idx, const CINTEnvVars *envs) } -void CINTg4c1e_ovlp(double *g, const CINTEnvVars *envs) +void CINTg4c1e_ovlp(double *g, CINTEnvVars *envs, double *cache) { const FINT li = envs->li_ceil; const FINT lj = envs->lj_ceil; @@ -150,19 +150,22 @@ void CINTg4c1e_ovlp(double *g, const CINTEnvVars *envs) const FINT nmax = li + lj; const FINT mmax = lk + ll; const FINT db = nmax + mmax + 1; - const FINT b_size = db*db; - double buf[3*b_size]; - double *bufx = buf; + const FINT b_size = db * (MAX(nmax, mmax) + 1); + double *bufx = cache; double *bufy = bufx + b_size; double *bufz = bufy + b_size; - double aij = envs->aij; - double akl = envs->akl; const double *rij = envs->rij; const double *rkl = envs->rkl; const double *ri = envs->rx_in_rijrx; const double *rk = envs->rx_in_rklrx; - double aijkl; + double ai = envs->ai[0]; + double aj = envs->aj[0]; + double ak = envs->ak[0]; + double al = envs->al[0]; + double aij = ai + aj; + double akl = ak + al; + double aijkl = aij + akl; double r1r12[3]; double r1r2[3]; FINT i, j, ptr; @@ -173,7 +176,6 @@ void CINTg4c1e_ovlp(double *g, const CINTEnvVars *envs) const FINT dn = envs->g2d_ijmax; const FINT dm = envs->g2d_klmax; - aijkl = aij + akl; bufx[0] = 1; bufy[0] = 1; bufz[0] = envs->fac[0] / (aijkl * sqrt(aijkl)); diff --git a/testsuite/test_cint4c1e.py b/testsuite/test_cint4c1e.py index bd6e839e..9a3e7fd8 100644 --- a/testsuite/test_cint4c1e.py +++ b/testsuite/test_cint4c1e.py @@ -171,7 +171,6 @@ def by_pp(shls, shape): ref+= (buf[:,:,0] + buf[:,:,4] + buf[:,:,8]).transpose(1,0) return ref * (-.25/numpy.pi) - #intor = _cint.cint4c1e_sph ao_loc = mol.ao_loc_nr() for nucid in range(mol.natm): mol.set_rinv_orig(coords[nucid]) @@ -233,4 +232,4 @@ def by_pp(shls, shape): print('cint4c1e_sph pass') test_int2c1e_sph() -#test_int4c1e_sph() +test_int4c1e_sph()