Skip to content

Commit

Permalink
Release v5.1.6
Browse files Browse the repository at this point in the history
Add 4-center 1-electron integral int4c1e
  • Loading branch information
sunqm committed Aug 27, 2022
1 parent 99b2514 commit 8f06588
Show file tree
Hide file tree
Showing 8 changed files with 56 additions and 63 deletions.
3 changes: 2 additions & 1 deletion .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand All @@ -34,3 +34,4 @@ jobs:
python test_3c2e.py --quick
python test_cart2sph.py
python test_c2s.py
python test_cint4c1e.py
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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}")
Expand Down
3 changes: 3 additions & 0 deletions ChangeLog
Original file line number Diff line number Diff line change
@@ -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

Expand Down
4 changes: 2 additions & 2 deletions README.rst
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
libcint
=======

version 5.1.5
2022-08-24
version 5.1.6
2022-08-26


What is libcint
Expand Down
33 changes: 7 additions & 26 deletions src/cint3c1e.c
Original file line number Diff line number Diff line change
Expand Up @@ -454,34 +454,15 @@ 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)
{
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,
Expand All @@ -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);
}

Expand All @@ -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);
}

Expand All @@ -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,
Expand All @@ -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);
}

Expand All @@ -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);
}

Expand Down
53 changes: 30 additions & 23 deletions src/cint4c1e.c
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down Expand Up @@ -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 *
Expand All @@ -156,23 +157,26 @@ 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 {
fac1l = common_fac;
*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 {
Expand All @@ -181,27 +185,30 @@ 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 {
fac1j = fac1k;
*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) {
Expand All @@ -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: ;
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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,
Expand All @@ -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);
}

Expand Down
18 changes: 10 additions & 8 deletions src/g4c1e.c
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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));
Expand Down
3 changes: 1 addition & 2 deletions testsuite/test_cint4c1e.py
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand Down Expand Up @@ -233,4 +232,4 @@ def by_pp(shls, shape):
print('cint4c1e_sph pass')

test_int2c1e_sph()
#test_int4c1e_sph()
test_int4c1e_sph()

0 comments on commit 8f06588

Please sign in to comment.