Skip to content

Second derivative integrals for int3c1e #129

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
maxscheurer wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Second derivative integrals for int3c1e #129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
04ef48d
chore: open PR for int3c1e second-derivative integrals
maxscheurer May 28, 2026
637c01e
Add int3c1e second-derivative integrals (ipip1, ip1ip2, ipvip1, ipip2)
maxscheurer May 28, 2026
1796016
Improve int3c1e hessian tests: high-l basis, ctypes tests, docs
maxscheurer May 28, 2026
56e784d
Merge pull request #1 from maxscheurer/feature/issue-2-int3c1e-second…
maxscheurer May 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 24 additions & 0 deletions include/cint_funcs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1066,6 +1066,30 @@ extern CINTIntegralFunction int3c1e_ip1_cart;
extern CINTIntegralFunction int3c1e_ip1_sph;
extern CINTIntegralFunction int3c1e_ip1_spinor;

/* 3-center 1-electron integral <(NABLA NABLA i) (j) (k)> */
extern CINTOptimizerFunction int3c1e_ipip1_optimizer;
extern CINTIntegralFunction int3c1e_ipip1_cart;
extern CINTIntegralFunction int3c1e_ipip1_sph;
extern CINTIntegralFunction int3c1e_ipip1_spinor;

/* 3-center 1-electron integral <(NABLA i) (j) (NABLA k)> */
extern CINTOptimizerFunction int3c1e_ip1ip2_optimizer;
extern CINTIntegralFunction int3c1e_ip1ip2_cart;
extern CINTIntegralFunction int3c1e_ip1ip2_sph;
extern CINTIntegralFunction int3c1e_ip1ip2_spinor;

/* 3-center 1-electron integral <(NABLA i) (NABLA j) (k)> */
extern CINTOptimizerFunction int3c1e_ipvip1_optimizer;
extern CINTIntegralFunction int3c1e_ipvip1_cart;
extern CINTIntegralFunction int3c1e_ipvip1_sph;
extern CINTIntegralFunction int3c1e_ipvip1_spinor;

/* 3-center 1-electron integral <(i) (j) (NABLA NABLA k)> */
extern CINTOptimizerFunction int3c1e_ipip2_optimizer;
extern CINTIntegralFunction int3c1e_ipip2_cart;
extern CINTIntegralFunction int3c1e_ipip2_sph;
extern CINTIntegralFunction int3c1e_ipip2_spinor;

/* <NABLA NABLA NABLA i|NUC |j> */
extern CINTOptimizerFunction int1e_ipipipnuc_optimizer;
extern CINTIntegralFunction int1e_ipipipnuc_cart;
Expand Down
4 changes: 4 additions & 0 deletions scripts/auto_intor.cl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,6 +216,10 @@
'("int3c1e_p2" ( \, \, p dot p))
'("int3c1e_iprinv" ( p \, \| rinv \| ))
'("int3c1e_ip1" ( nabla \, \,))
'("int3c1e_ipip1" ( nabla nabla \, \,))
'("int3c1e_ip1ip2" ( nabla \, \, nabla))
'("int3c1e_ipvip1" ( nabla \, nabla \,))
'("int3c1e_ipip2" ( \, \, nabla nabla))
)

(gen-cint "deriv3.c"
Expand Down
308 changes: 308 additions & 0 deletions src/autocode/int3c1e.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -185,3 +185,311 @@ return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0, 0);
} // int3c1e_ip1_spinor
ALL_CINT(int3c1e_ip1)
ALL_CINT_FORTRAN_(int3c1e_ip1)
static void CINTgout1e_int3c1e_ipip1(double *gout,
double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
FINT nf = envs->nf;
FINT ix, iy, iz, n;
double *g0 = g;
double *g1 = g0 + envs->g_size * 3;
double *g2 = g1 + envs->g_size * 3;
double *g3 = g2 + envs->g_size * 3;
double s[9];
G1E_D_I(g1, g0, envs->i_l+1, envs->j_l, envs->k_l);
G1E_D_I(g2, g0, envs->i_l+0, envs->j_l, envs->k_l);
G1E_D_I(g3, g1, envs->i_l+0, envs->j_l, envs->k_l);
for (n = 0; n < nf; n++) {
ix = idx[0+n*3];
iy = idx[1+n*3];
iz = idx[2+n*3];
s[0] = + g3[ix+0]*g0[iy+0]*g0[iz+0];
s[1] = + g2[ix+0]*g1[iy+0]*g0[iz+0];
s[2] = + g2[ix+0]*g0[iy+0]*g1[iz+0];
s[3] = + g1[ix+0]*g2[iy+0]*g0[iz+0];
s[4] = + g0[ix+0]*g3[iy+0]*g0[iz+0];
s[5] = + g0[ix+0]*g2[iy+0]*g1[iz+0];
s[6] = + g1[ix+0]*g0[iy+0]*g2[iz+0];
s[7] = + g0[ix+0]*g1[iy+0]*g2[iz+0];
s[8] = + g0[ix+0]*g0[iy+0]*g3[iz+0];
if (gout_empty) {
gout[n*9+0] = + s[0];
gout[n*9+1] = + s[3];
gout[n*9+2] = + s[6];
gout[n*9+3] = + s[1];
gout[n*9+4] = + s[4];
gout[n*9+5] = + s[7];
gout[n*9+6] = + s[2];
gout[n*9+7] = + s[5];
gout[n*9+8] = + s[8];
} else {
gout[n*9+0] += + s[0];
gout[n*9+1] += + s[3];
gout[n*9+2] += + s[6];
gout[n*9+3] += + s[1];
gout[n*9+4] += + s[4];
gout[n*9+5] += + s[7];
gout[n*9+6] += + s[2];
gout[n*9+7] += + s[5];
gout[n*9+8] += + s[8];
}}}
void int3c1e_ipip1_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {
FINT ng[] = {2, 0, 0, 0, 2, 1, 1, 9};
CINTall_3c1e_optimizer(opt, ng, atm, natm, bas, nbas, env);
}
CACHE_SIZE_T int3c1e_ipip1_cart(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {2, 0, 0, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipip1;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, 0, 0);
}
// int3c1e_ipip1_cart
CACHE_SIZE_T int3c1e_ipip1_sph(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {2, 0, 0, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipip1;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, 0, 0);
} // int3c1e_ipip1_sph
CACHE_SIZE_T int3c1e_ipip1_spinor(double complex *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {2, 0, 0, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipip1;
return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0, 0);
} // int3c1e_ipip1_spinor
ALL_CINT(int3c1e_ipip1)
ALL_CINT_FORTRAN_(int3c1e_ipip1)
static void CINTgout1e_int3c1e_ip1ip2(double *gout,
double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
FINT nf = envs->nf;
FINT ix, iy, iz, n;
double *g0 = g;
double *g1 = g0 + envs->g_size * 3;
double *g2 = g1 + envs->g_size * 3;
double *g3 = g2 + envs->g_size * 3;
double s[9];
G1E_D_K(g1, g0, envs->i_l+1, envs->j_l+0, envs->k_l+0);
G1E_D_I(g2, g0, envs->i_l+0, envs->j_l, envs->k_l);
G1E_D_I(g3, g1, envs->i_l+0, envs->j_l, envs->k_l);
for (n = 0; n < nf; n++) {
ix = idx[0+n*3];
iy = idx[1+n*3];
iz = idx[2+n*3];
s[0] = + g3[ix+0]*g0[iy+0]*g0[iz+0];
s[1] = + g2[ix+0]*g1[iy+0]*g0[iz+0];
s[2] = + g2[ix+0]*g0[iy+0]*g1[iz+0];
s[3] = + g1[ix+0]*g2[iy+0]*g0[iz+0];
s[4] = + g0[ix+0]*g3[iy+0]*g0[iz+0];
s[5] = + g0[ix+0]*g2[iy+0]*g1[iz+0];
s[6] = + g1[ix+0]*g0[iy+0]*g2[iz+0];
s[7] = + g0[ix+0]*g1[iy+0]*g2[iz+0];
s[8] = + g0[ix+0]*g0[iy+0]*g3[iz+0];
if (gout_empty) {
gout[n*9+0] = + s[0];
gout[n*9+1] = + s[1];
gout[n*9+2] = + s[2];
gout[n*9+3] = + s[3];
gout[n*9+4] = + s[4];
gout[n*9+5] = + s[5];
gout[n*9+6] = + s[6];
gout[n*9+7] = + s[7];
gout[n*9+8] = + s[8];
} else {
gout[n*9+0] += + s[0];
gout[n*9+1] += + s[1];
gout[n*9+2] += + s[2];
gout[n*9+3] += + s[3];
gout[n*9+4] += + s[4];
gout[n*9+5] += + s[5];
gout[n*9+6] += + s[6];
gout[n*9+7] += + s[7];
gout[n*9+8] += + s[8];
}}}
void int3c1e_ip1ip2_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {
FINT ng[] = {1, 0, 1, 0, 2, 1, 1, 9};
CINTall_3c1e_optimizer(opt, ng, atm, natm, bas, nbas, env);
}
CACHE_SIZE_T int3c1e_ip1ip2_cart(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {1, 0, 1, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ip1ip2;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, 0, 0);
}
// int3c1e_ip1ip2_cart
CACHE_SIZE_T int3c1e_ip1ip2_sph(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {1, 0, 1, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ip1ip2;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, 0, 0);
} // int3c1e_ip1ip2_sph
CACHE_SIZE_T int3c1e_ip1ip2_spinor(double complex *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {1, 0, 1, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ip1ip2;
return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0, 0);
} // int3c1e_ip1ip2_spinor
ALL_CINT(int3c1e_ip1ip2)
ALL_CINT_FORTRAN_(int3c1e_ip1ip2)
static void CINTgout1e_int3c1e_ipvip1(double *gout,
double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
FINT nf = envs->nf;
FINT ix, iy, iz, n;
double *g0 = g;
double *g1 = g0 + envs->g_size * 3;
double *g2 = g1 + envs->g_size * 3;
double *g3 = g2 + envs->g_size * 3;
double s[9];
G1E_D_J(g1, g0, envs->i_l+1, envs->j_l+0, envs->k_l);
G1E_D_I(g2, g0, envs->i_l+0, envs->j_l, envs->k_l);
G1E_D_I(g3, g1, envs->i_l+0, envs->j_l, envs->k_l);
for (n = 0; n < nf; n++) {
ix = idx[0+n*3];
iy = idx[1+n*3];
iz = idx[2+n*3];
s[0] = + g3[ix+0]*g0[iy+0]*g0[iz+0];
s[1] = + g2[ix+0]*g1[iy+0]*g0[iz+0];
s[2] = + g2[ix+0]*g0[iy+0]*g1[iz+0];
s[3] = + g1[ix+0]*g2[iy+0]*g0[iz+0];
s[4] = + g0[ix+0]*g3[iy+0]*g0[iz+0];
s[5] = + g0[ix+0]*g2[iy+0]*g1[iz+0];
s[6] = + g1[ix+0]*g0[iy+0]*g2[iz+0];
s[7] = + g0[ix+0]*g1[iy+0]*g2[iz+0];
s[8] = + g0[ix+0]*g0[iy+0]*g3[iz+0];
if (gout_empty) {
gout[n*9+0] = + s[0];
gout[n*9+1] = + s[1];
gout[n*9+2] = + s[2];
gout[n*9+3] = + s[3];
gout[n*9+4] = + s[4];
gout[n*9+5] = + s[5];
gout[n*9+6] = + s[6];
gout[n*9+7] = + s[7];
gout[n*9+8] = + s[8];
} else {
gout[n*9+0] += + s[0];
gout[n*9+1] += + s[1];
gout[n*9+2] += + s[2];
gout[n*9+3] += + s[3];
gout[n*9+4] += + s[4];
gout[n*9+5] += + s[5];
gout[n*9+6] += + s[6];
gout[n*9+7] += + s[7];
gout[n*9+8] += + s[8];
}}}
void int3c1e_ipvip1_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {
FINT ng[] = {1, 1, 0, 0, 2, 1, 1, 9};
CINTall_3c1e_optimizer(opt, ng, atm, natm, bas, nbas, env);
}
CACHE_SIZE_T int3c1e_ipvip1_cart(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {1, 1, 0, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipvip1;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, 0, 0);
}
// int3c1e_ipvip1_cart
CACHE_SIZE_T int3c1e_ipvip1_sph(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {1, 1, 0, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipvip1;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, 0, 0);
} // int3c1e_ipvip1_sph
CACHE_SIZE_T int3c1e_ipvip1_spinor(double complex *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {1, 1, 0, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipvip1;
return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0, 0);
} // int3c1e_ipvip1_spinor
ALL_CINT(int3c1e_ipvip1)
ALL_CINT_FORTRAN_(int3c1e_ipvip1)
static void CINTgout1e_int3c1e_ipip2(double *gout,
double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
FINT nf = envs->nf;
FINT ix, iy, iz, n;
double *g0 = g;
double *g1 = g0 + envs->g_size * 3;
double *g2 = g1 + envs->g_size * 3;
double *g3 = g2 + envs->g_size * 3;
double s[9];
G1E_D_K(g1, g0, envs->i_l+0, envs->j_l+0, envs->k_l+1);
G1E_D_K(g2, g0, envs->i_l+0, envs->j_l+0, envs->k_l+0);
G1E_D_K(g3, g1, envs->i_l+0, envs->j_l+0, envs->k_l+0);
for (n = 0; n < nf; n++) {
ix = idx[0+n*3];
iy = idx[1+n*3];
iz = idx[2+n*3];
s[0] = + g3[ix+0]*g0[iy+0]*g0[iz+0];
s[1] = + g2[ix+0]*g1[iy+0]*g0[iz+0];
s[2] = + g2[ix+0]*g0[iy+0]*g1[iz+0];
s[3] = + g1[ix+0]*g2[iy+0]*g0[iz+0];
s[4] = + g0[ix+0]*g3[iy+0]*g0[iz+0];
s[5] = + g0[ix+0]*g2[iy+0]*g1[iz+0];
s[6] = + g1[ix+0]*g0[iy+0]*g2[iz+0];
s[7] = + g0[ix+0]*g1[iy+0]*g2[iz+0];
s[8] = + g0[ix+0]*g0[iy+0]*g3[iz+0];
if (gout_empty) {
gout[n*9+0] = + s[0];
gout[n*9+1] = + s[3];
gout[n*9+2] = + s[6];
gout[n*9+3] = + s[1];
gout[n*9+4] = + s[4];
gout[n*9+5] = + s[7];
gout[n*9+6] = + s[2];
gout[n*9+7] = + s[5];
gout[n*9+8] = + s[8];
} else {
gout[n*9+0] += + s[0];
gout[n*9+1] += + s[3];
gout[n*9+2] += + s[6];
gout[n*9+3] += + s[1];
gout[n*9+4] += + s[4];
gout[n*9+5] += + s[7];
gout[n*9+6] += + s[2];
gout[n*9+7] += + s[5];
gout[n*9+8] += + s[8];
}}}
void int3c1e_ipip2_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {
FINT ng[] = {0, 0, 2, 0, 2, 1, 1, 9};
CINTall_3c1e_optimizer(opt, ng, atm, natm, bas, nbas, env);
}
CACHE_SIZE_T int3c1e_ipip2_cart(double *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {0, 0, 2, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipip2;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_cart_3c1e, 0, 0);
}
// int3c1e_ipip2_cart
CACHE_SIZE_T int3c1e_ipip2_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, 2, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipip2;
return CINT3c1e_drv(out, dims, &envs, opt, cache, &c2s_sph_3c1e, 0, 0);
} // int3c1e_ipip2_sph
CACHE_SIZE_T int3c1e_ipip2_spinor(double complex *out, FINT *dims, FINT *shls,
FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
FINT ng[] = {0, 0, 2, 0, 2, 1, 1, 9};
CINTEnvVars envs;
CINTinit_int3c1e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
envs.f_gout = &CINTgout1e_int3c1e_ipip2;
return CINT3c1e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_3c2e1, 0, 0);
} // int3c1e_ipip2_spinor
ALL_CINT(int3c1e_ipip2)
ALL_CINT_FORTRAN_(int3c1e_ipip2)
4 changes: 4 additions & 0 deletions testsuite/test_int3c1e.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,3 +103,7 @@ def run(intor, comp=1, suffix='_sph', thr=1e-7):
run('int3c1e_r2_origk')
run('int3c1e_r4_origk')
run('int3c1e_r6_origk', thr=1e-5)
run('int3c1e_ipip1', comp=9)
run('int3c1e_ip1ip2', comp=9)
run('int3c1e_ipvip1', comp=9)
run('int3c1e_ipip2', comp=9)
Loading