Skip to content

Commit

Permalink
Change SR integral algorithm for nroots <= 3
Browse files Browse the repository at this point in the history 
Refactor g2e.c

Remove GTG integrals and coulerf function

Add root finder using eigenvalue algorithm
  • Loading branch information
@sunqm
sunqm committed Sep 19, 2023
1 parent 839f0ec commit c5eb4a7
Show file tree
Hide file tree
Showing 21 changed files with 6,617 additions and 2,415 deletions.
32 changes: 12 additions & 20 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
cmake_minimum_required (VERSION 3.5)
project (qcint C)
set(qcint_VERSION_MAJOR "5")
set(qcint_VERSION_MINOR "5")
set(qcint_VERSION_MAJOR "6")
set(qcint_VERSION_MINOR "0")
set(qcint_VERSION_PATCH "0")
set(qcint_VERSION_TWEAK "0")
set(qcint_VERSION "${qcint_VERSION_MAJOR}.${qcint_VERSION_MINOR}.${qcint_VERSION_PATCH}")
set(qcint_SOVERSION "${qcint_VERSION_MAJOR}")
set(cint_SOVERSION "${qcint_VERSION_MAJOR}")

#set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -g -O2 -DNDEBUG")
if ("${CMAKE_BUILD_TYPE}" STREQUAL "")
Expand Down Expand Up @@ -71,7 +71,7 @@ set(cintSrc
src/g1e.c src/g2e.c src/g2e_simd1.c src/g3c1e.c src/g3c2e.c
src/gout2e.c src/misc.c src/optimizer.c
src/fmt.c src/rys_wheeler.c src/eigh.c src/rys_roots.c src/find_roots.c
src/polyfits.c src/sr_rys_polyfits.c
src/polyfits.c
src/cint1e_a.c src/cint3c1e_a.c
src/cint1e_grids.c src/g1e_grids.c
src/autocode/breit1.c src/autocode/dkb.c src/autocode/gaunt1.c
Expand All @@ -84,17 +84,17 @@ set(cintSrc
if(WITH_RANGE_COULOMB)
# defined in config.h
# add_definitions(-DWITH_RANGE_COULOMB)
message("Enabled WITH_RANGE_COULOMB")
if(WITH_POLYNOMIAL_FIT)
add_definitions(-DWITH_POLYNOMIAL_FIT)
message("Enabled WITH_POLYNOMIAL_FIT")
endif(WITH_POLYNOMIAL_FIT)
# message("Enabled WITH_RANGE_COULOMB")
endif(WITH_RANGE_COULOMB)

if(WITH_POLYNOMIAL_FIT)
set(cintSrc ${cintSrc} src/sr_rys_polyfits.c)
add_definitions(-DWITH_POLYNOMIAL_FIT)
message("Enabled WITH_POLYNOMIAL_FIT")
endif(WITH_POLYNOMIAL_FIT)

if(WITH_COULOMB_ERF)
set(cintSrc ${cintSrc} src/g2e_coulerf.c src/cint2e_coulerf.c)
add_definitions(-DWITH_COULOMB_ERF)
message("Enabled WITH_COULOMB_ERF")
message("WITH_COULOMB_ERF is deprecated since v6.0")
endif(WITH_COULOMB_ERF)

if(WITH_F12)
Expand All @@ -103,14 +103,6 @@ if(WITH_F12)
message("Enabled WITH_F12")
endif(WITH_F12)

if(WITH_GTG)
set(cintSrc ${cintSrc} src/g2e_gtg.c src/cint2e_gtg.c src/cint3c2e_gtg.c
src/cint2c2e_gtg.c)
add_definitions(-DWITH_GTG)
message("Enabled WITH_GTG")
message("Enabled WITH_GTG. Note there are bugs in gtg type integrals")
endif(WITH_GTG)

if(PYPZPX)
add_definitions(-DPYPZPX)
message("P orbitals convention (py, pz, px)")
Expand Down
5 changes: 3 additions & 2 deletions include/cint.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@
* Parameters and function prototypes for libcint.
*/

#define QCINT_VERSION @qcint_VERSION@
#define QCINT_VERSION @qcint_VERSION@
#define CINT_SOVERSION @cint_SOVERSION

#cmakedefine CACHE_SIZE_I8
#ifdef CACHE_SIZE_I8
Expand Down Expand Up @@ -203,7 +204,7 @@ typedef struct {
union {int nfk; int grids_offset;};
union {int nfl; int ngrids;};
int nf; // = nfi*nfj*nfk*nfl;
int _padding;
int rys_order; // = nrys_roots for regular ERIs. can be nrys_roots/2 for SR ERIs
int x_ctr[4];

int gbits;
Expand Down
115 changes: 77 additions & 38 deletions src/cint2c2e.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,8 @@
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <assert.h>
#include "cint_bas.h"
#include "misc.h"
#include "g2e.h"
Expand Down Expand Up @@ -65,15 +65,38 @@
(*(fp2c[i]))(gctr[it], gprim[i], coeff[it]+im, \
ngp[it], x_prim[it], x_ctr[it], \
non0ctr[it][im], non0idx[it]+im*x_ctr[it]); \
empty_overall = 0; \
} \
cum = 0; \
np2c = 0;

#define POP_PRIM2CTR_AND_SET0 \
for (i = 0; i < np2c; i++) { \
it = shltyp[i]; \
if (it != SHLTYPi) { \
im = iprim[i]; \
(*(fp2c[i]))(gctr[i], gprim[i], coeff[it]+im, \
ngp[it], x_prim[it], x_ctr[it], \
non0ctr[it][im], non0idx[it]+im*x_ctr[it]); \
empty_overall = 0; \
} else if (fp2c[i] == CINTiprim_to_ctr_0) { \
double *pout = gctr[it]; \
for (int k = 0; k < nf; k++) { \
pout[k] = 0.; \
} \
} \
} \
cum = 0; \
np2c = 0;

#define PUSH \
if (cum == SIMDD) { \
(*envs->f_g0_2e)(g, cutoff, &bc, envs, cum); \
(*envs->f_gout)(gout, g, idx, envs); \
POP_PRIM2CTR; \
if ((*envs->f_g0_2e)(g, cutoff, &bc, envs, cum)) { \
(*envs->f_gout)(gout, g, idx, envs); \
POP_PRIM2CTR; \
} else { \
POP_PRIM2CTR_AND_SET0; \
} \
} \
envs->ai[cum] = ai[ip]; \
envs->ak[cum] = ak[kp]; \
Expand Down Expand Up @@ -111,26 +134,30 @@

#define RUN_REST \
if (cum == 1) { \
(*envs->f_g0_2e_simd1)(g, cutoff, &bc, envs, 0); \
(*envs->f_gout_simd1)(gout, g, idx, envs); \
if ((*envs->f_g0_2e_simd1)(g, cutoff, &bc, envs, 0)) { \
(*envs->f_gout_simd1)(gout, g, idx, envs); \
POP_PRIM2CTR; \
} else { \
POP_PRIM2CTR_AND_SET0; \
} \
} else if (cum > 1) { \
r1 = MM_SET1(1.); \
for (i = 0; i < envs->nrys_roots; i++) { \
MM_STORE(bc.u+i*SIMDD, r1); \
MM_STORE(bc.w+i*SIMDD, r1); \
if ((*envs->f_g0_2e)(g, cutoff, &bc, envs, cum)) { \
(*envs->f_gout)(gout, g, idx, envs); \
POP_PRIM2CTR; \
} else { \
POP_PRIM2CTR_AND_SET0; \
} \
(*envs->f_g0_2e)(g, cutoff, &bc, envs, cum); \
(*envs->f_gout)(gout, g, idx, envs); \
} \
POP_PRIM2CTR;
} else { \
assert(np2c == 0); \
}

#define TRANSPOSE(a) \
if (*empty) { \
CINTdmat_transpose(out, a, nf*nc, n_comp); \
*empty = 0; \
} else { \
CINTdplus_transpose(out, a, nf*nc, n_comp); \
} \
*empty = 0;
}

int CINT2c2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache, int *empty)
{
Expand Down Expand Up @@ -158,6 +185,7 @@ int CINT2c2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache, int *empty
int _empty[2] = {1, 1};
int *iempty = _empty + 0;
int *kempty = _empty + 1;
int empty_overall = 1;
int ngp[2];
ngp[0] = nf * n_comp;
ngp[1] = ngp[0] * i_ctr;
Expand All @@ -171,7 +199,7 @@ int CINT2c2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache, int *empty
if (n_comp == 1) {
MALLOC_INSTACK(g1, lenk);
gctr[SHLTYPk] = out;
kempty = empty;
*kempty = *empty;
} else {
MALLOC_INSTACK(gctr[SHLTYPk], lenk);
g1 = out;
Expand All @@ -195,8 +223,7 @@ int CINT2c2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache, int *empty
CINTOpt_non0coeff_byshell(non0idxk, non0ctrk, coeff[1], k_prim, k_ctr);
int *non0ctr[2] = {non0ctri, non0ctrk};
int *non0idx[2] = {non0idxi, non0idxk};
double common_factor = envs->common_factor * (M_PI*M_PI*M_PI)*2/SQRTPI
* CINTcommon_fac_sp(envs->i_l) * CINTcommon_fac_sp(envs->k_l);
double common_factor = envs->common_factor;
INITSIMD;

for (kp = 0; kp < k_prim; kp++) {
Expand All @@ -215,11 +242,12 @@ int CINT2c2e_loop_nopt(double *out, CINTEnvVars *envs, double *cache, int *empty
} // end loop k_prim
RUN_REST;

if (n_comp > 1 && !*kempty) {
if (n_comp > 1 && !empty_overall) {
int nc = i_ctr * k_ctr;
TRANSPOSE(gctr[SHLTYPk]);
}
return !*empty;
*empty &= empty_overall;
return !empty_overall;
}


Expand Down Expand Up @@ -249,6 +277,7 @@ int CINT2c2e_loop(double *out, CINTEnvVars *envs, double *cache, int *empty)
int _empty[2] = {1, 1};
int *iempty = _empty + 0;
int *kempty = _empty + 1;
int empty_overall = 1;
int ngp[2];
ngp[0] = nf * n_comp;
ngp[1] = ngp[0] * i_ctr;
Expand All @@ -262,7 +291,7 @@ int CINT2c2e_loop(double *out, CINTEnvVars *envs, double *cache, int *empty)
if (n_comp == 1) {
MALLOC_INSTACK(g1, lenk);
gctr[SHLTYPk] = out;
kempty = empty;
*kempty = *empty;
} else {
MALLOC_INSTACK(gctr[SHLTYPk], lenk);
g1 = out;
Expand All @@ -273,8 +302,7 @@ int CINT2c2e_loop(double *out, CINTEnvVars *envs, double *cache, int *empty)

ALIGNMM Rys2eT bc;
ALIGNMM double cutoff[SIMDD];
double common_factor = envs->common_factor * (M_PI*M_PI*M_PI)*2/SQRTPI
* CINTcommon_fac_sp(envs->i_l) * CINTcommon_fac_sp(envs->k_l);
double common_factor = envs->common_factor;

CINTOpt *opt = envs->opt;
int *idx = opt->index_xyz_array[envs->i_l*LMAX1+envs->k_l];
Expand All @@ -299,11 +327,12 @@ int CINT2c2e_loop(double *out, CINTEnvVars *envs, double *cache, int *empty)
} // end loop k_prim
RUN_REST;

if (n_comp > 1 && !*kempty) {
if (n_comp > 1 && !empty_overall) {
int nc = i_ctr * k_ctr;
TRANSPOSE(gctr[SHLTYPk]);
}
return !*empty;
*empty &= empty_overall;
return !empty_overall;
}

void CINTinit_int2c2e_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
Expand Down Expand Up @@ -331,7 +360,12 @@ void CINTinit_int2c2e_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
envs->nfk = (envs->k_l+1)*(envs->k_l+2)/2;
envs->nfl = 1;
envs->nf = envs->nfi * envs->nfk;
envs->common_factor = 1;

envs->ri = env + atm(PTR_COORD, bas(ATOM_OF, i_sh));
envs->rk = env + atm(PTR_COORD, bas(ATOM_OF, k_sh));

envs->common_factor = (M_PI*M_PI*M_PI)*2/SQRTPI
* CINTcommon_fac_sp(envs->i_l) * CINTcommon_fac_sp(envs->k_l);
if (env[PTR_EXPCUTOFF] == 0) {
envs->expcutoff = EXPCUTOFF;
} else {
Expand All @@ -347,20 +381,21 @@ void CINTinit_int2c2e_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
envs->lj_ceil = 0;
envs->lk_ceil = envs->k_l + ng[KINC];
envs->ll_ceil = 0;

envs->ri = env + atm(PTR_COORD, bas(ATOM_OF, i_sh));
envs->rk = env + atm(PTR_COORD, bas(ATOM_OF, k_sh));

int nroots = (envs->li_ceil + envs->lk_ceil)/2 + 1;
envs->nrys_roots = nroots;
assert(nroots < MXRYSROOTS);
int rys_order =(envs->li_ceil + envs->lk_ceil)/2 + 1;
int nrys_roots = rys_order;
double omega = env[PTR_RANGE_OMEGA];
if (omega < 0 && rys_order <= 3) {
nrys_roots *= 2;
}
envs->rys_order = rys_order;
envs->nrys_roots = nrys_roots;

int dli = envs->li_ceil + 1;
int dlk = envs->lk_ceil + 1;
envs->g_stride_i = nroots;
envs->g_stride_k = nroots * dli;
envs->g_stride_i = nrys_roots;
envs->g_stride_k = nrys_roots * dli;
envs->g_stride_l = envs->g_stride_k;
envs->g_size = nroots * dli * dlk;
envs->g_size = nrys_roots * dli * dlk;

MM_STORE(envs->aj, MM_SET1(0.));
MM_STORE(envs->al, MM_SET1(0.));
Expand All @@ -381,9 +416,13 @@ void CINTinit_int2c2e_EnvVars(CINTEnvVars *envs, int *ng, int *shls,
envs->rx_in_rklrx = envs->rk;
envs->rx_in_rijrx = envs->ri;

if (nroots <= 2) {
if (rys_order <= 2) {
envs->f_g0_2d4d = &CINTg0_2e_2d4d_unrolled;
envs->f_g0_2d4d_simd1 = &CINTg0_2e_2d4d_unrolled_simd1;
if (rys_order != nrys_roots) {
envs->f_g0_2d4d = &CINTsrg0_2e_2d4d_unrolled;
envs->f_g0_2d4d_simd1 = &CINTsrg0_2e_2d4d_unrolled_simd1;
}
} else {
envs->f_g0_2d4d = &CINTg0_2e_2d;
envs->f_g0_2d4d_simd1 = &CINTg0_2e_2d_simd1;
Expand Down
Loading

0 comments on commit c5eb4a7

Please sign in to comment.