diff --git a/include/cint_funcs.h b/include/cint_funcs.h index 9d65fb21..03ace0db 100644 --- a/include/cint_funcs.h +++ b/include/cint_funcs.h @@ -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; + /* */ extern CINTOptimizerFunction int1e_ipipipnuc_optimizer; extern CINTIntegralFunction int1e_ipipipnuc_cart; diff --git a/scripts/auto_intor.cl b/scripts/auto_intor.cl index ad80afae..a8a24b9b 100644 --- a/scripts/auto_intor.cl +++ b/scripts/auto_intor.cl @@ -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" diff --git a/src/autocode/int3c1e.c b/src/autocode/int3c1e.c index 9065c860..db1c6461 100644 --- a/src/autocode/int3c1e.c +++ b/src/autocode/int3c1e.c @@ -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) diff --git a/testsuite/test_int3c1e.py b/testsuite/test_int3c1e.py index cc3afeb7..38ebc6c6 100644 --- a/testsuite/test_int3c1e.py +++ b/testsuite/test_int3c1e.py @@ -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) diff --git a/testsuite/test_int3c1e_hess.py b/testsuite/test_int3c1e_hess.py new file mode 100644 index 00000000..d2e55e40 --- /dev/null +++ b/testsuite/test_int3c1e_hess.py @@ -0,0 +1,377 @@ +""" +Finite-difference validation of int3c1e second-derivative integrals. + +Sign convention: libcint's nabla = d/dr (electron coord derivative). +Nuclear coordinate derivative: d/dR = -d/dr. +Therefore: FD of ip1 w.r.t. nuclear coord = -analytical_second_derivative. + +Component ordering: for (nabla_A nabla_B), component index = alpha_A * 3 + beta_B +where alpha_A is the direction of the first nabla and beta_B is the second. +""" + +import os +import numpy as np +from pyscf import gto + + +def _make_mol(): + """Create a test molecule with 3 atoms and high angular momentum (up to g). + + Uses a custom basis with s, p, d, f, g shells to stress-test the + G1E_D_I/G1E_D_J recursion at high l on bra/ket centers. + """ + basis_str = gto.basis.parse(''' +H S + 5.0250000 1.0000000 +H S + 1.0130000 1.0000000 +H P + 4.1330000 0.0868660 0.0000000 + 0.3827000 0.5010080 1.0000000 +H D + 1.0970000 1.0000000 +H F + 0.7610000 1.0000000 +H G + 0.8827000 1.0000000 + ''') + return gto.M( + atom=''' + H 0.0 0.0 0.0 + H 1.5 0.3 0.2 + H 0.5 1.8 0.4 + ''', + basis={'H': basis_str}, + unit='Bohr', + spin=1, + ) + + +def _make_auxmol(): + """Aux mol with s, p, d, f shells on a single atom (separate from mol atoms).""" + basis_str = { + 'He': gto.basis.parse(''' +He S + 1.5 1.0 +He P + 0.8 1.0 +He D + 1.2 1.0 +He F + 0.9 1.0 + '''), + } + return gto.M( + atom='He 2.0 0.5 1.0', + basis=basis_str, + unit='Bohr', + ) + + +def _compute_3c1e(mol, auxmol, intor, comp=1): + """Compute a 3c1e integral using mol+auxmol as supermol.""" + supermol = mol + auxmol + nbas_mol = mol.nbas + nbas_aux = auxmol.nbas + shls_slice = (0, nbas_mol, 0, nbas_mol, nbas_mol, nbas_mol + nbas_aux) + return supermol.intor(intor, comp=comp, shls_slice=shls_slice) + + +def test_int3c1e_ipip1(): + """ + Verify int3c1e_ipip1 via FD of int3c1e_ip1 w.r.t. bra atom. + + For shell triple (i on atom A, j NOT on A, k NOT on A): + d/dR_A,β [ip1[α, i, j, k]] = -ipip1[α*3+β, i, j, k] + + Note: Co-centered pairs (i and j on the same atom) are excluded because + displacing atom A moves both bra and ket shells simultaneously, mixing + ipip1 and ipvip1 contributions in the FD. Co-centered elements are + validated indirectly via test_translational_invariance (ipip1 + ipvip1 + + ip1ip2 = 0 for ALL elements including co-centered ones). + """ + mol = _make_mol() + auxmol = _make_auxmol() + delta = 1e-5 + + analytic = _compute_3c1e(mol, auxmol, 'int3c1e_ipip1', comp=9) + nao = mol.nao + naux = auxmol.nao + assert analytic.shape == (9, nao, nao, naux) + + max_err = 0.0 + aoslice = mol.aoslice_by_atom() + + for iatm in range(mol.natm): + i0, i1 = int(aoslice[iatm, 2]), int(aoslice[iatm, 3]) + # j indices NOT on iatm + j_mask = np.ones(nao, dtype=bool) + j_mask[i0:i1] = False + + for beta in range(3): + coords_p = mol.atom_coords(unit='Bohr').copy() + coords_m = mol.atom_coords(unit='Bohr').copy() + coords_p[iatm, beta] += delta + coords_m[iatm, beta] -= delta + mol_p = mol.set_geom_(coords_p, unit='Bohr', inplace=False) + mol_m = mol.set_geom_(coords_m, unit='Bohr', inplace=False) + ip1_p = _compute_3c1e(mol_p, auxmol, 'int3c1e_ip1', comp=3) + ip1_m = _compute_3c1e(mol_m, auxmol, 'int3c1e_ip1', comp=3) + dip1 = (ip1_p - ip1_m) / (2 * delta) + + for alpha in range(3): + comp_idx = alpha * 3 + beta + # FD = -ipip1 + ana_block = analytic[comp_idx, i0:i1][:, j_mask, :] + fd_block = -dip1[alpha, i0:i1][:, j_mask, :] + err = np.max(np.abs(ana_block - fd_block)) + max_err = max(max_err, err) + + print(f"int3c1e_ipip1: max FD error = {max_err:.2e}") + assert max_err < 1e-7, f"int3c1e_ipip1 FD error too large: {max_err}" + + +def test_int3c1e_ip1ip2(): + """ + Verify int3c1e_ip1ip2 by displacing aux atom. + + Since aux atoms are distinct from mol atoms: + d/dR_C,β [ip1[α, i, j, k]] = -ip1ip2[α*3+β, i, j, k] for k on C + """ + mol = _make_mol() + auxmol = _make_auxmol() + delta = 1e-5 + + analytic = _compute_3c1e(mol, auxmol, 'int3c1e_ip1ip2', comp=9) + nao = mol.nao + naux = auxmol.nao + assert analytic.shape == (9, nao, nao, naux) + + max_err = 0.0 + aoslice_aux = auxmol.aoslice_by_atom() + + for katm in range(auxmol.natm): + k0, k1 = int(aoslice_aux[katm, 2]), int(aoslice_aux[katm, 3]) + for beta in range(3): + coords_p = auxmol.atom_coords(unit='Bohr').copy() + coords_m = auxmol.atom_coords(unit='Bohr').copy() + coords_p[katm, beta] += delta + coords_m[katm, beta] -= delta + auxmol_p = auxmol.set_geom_(coords_p, unit='Bohr', inplace=False) + auxmol_m = auxmol.set_geom_(coords_m, unit='Bohr', inplace=False) + ip1_p = _compute_3c1e(mol, auxmol_p, 'int3c1e_ip1', comp=3) + ip1_m = _compute_3c1e(mol, auxmol_m, 'int3c1e_ip1', comp=3) + dip1 = (ip1_p - ip1_m) / (2 * delta) + + for alpha in range(3): + comp_idx = alpha * 3 + beta + ana_block = analytic[comp_idx, :, :, k0:k1] + fd_block = -dip1[alpha, :, :, k0:k1] + err = np.max(np.abs(ana_block - fd_block)) + max_err = max(max_err, err) + + print(f"int3c1e_ip1ip2: max FD error = {max_err:.2e}") + assert max_err < 1e-7, f"int3c1e_ip1ip2 FD error too large: {max_err}" + + +def test_int3c1e_ipvip1(): + """ + Verify int3c1e_ipvip1 by displacing ket atom. + + For shell triple (i NOT on B, j on B, k NOT on B): + d/dR_B,β [ip1[α, i, j, k]] = -ipvip1[α*3+β, i, j, k] + + Note: Co-centered pairs (i and j on the same atom) are excluded because + displacing atom B moves both bra and ket shells simultaneously, mixing + ipvip1 and ipip1 contributions in the FD. See test_translational_invariance + for validation of co-centered elements. + """ + mol = _make_mol() + auxmol = _make_auxmol() + delta = 1e-5 + + analytic = _compute_3c1e(mol, auxmol, 'int3c1e_ipvip1', comp=9) + nao = mol.nao + naux = auxmol.nao + assert analytic.shape == (9, nao, nao, naux) + + max_err = 0.0 + aoslice = mol.aoslice_by_atom() + + for jatm in range(mol.natm): + j0, j1 = int(aoslice[jatm, 2]), int(aoslice[jatm, 3]) + # i indices NOT on jatm + i_mask = np.ones(nao, dtype=bool) + i_mask[j0:j1] = False + + for beta in range(3): + coords_p = mol.atom_coords(unit='Bohr').copy() + coords_m = mol.atom_coords(unit='Bohr').copy() + coords_p[jatm, beta] += delta + coords_m[jatm, beta] -= delta + mol_p = mol.set_geom_(coords_p, unit='Bohr', inplace=False) + mol_m = mol.set_geom_(coords_m, unit='Bohr', inplace=False) + ip1_p = _compute_3c1e(mol_p, auxmol, 'int3c1e_ip1', comp=3) + ip1_m = _compute_3c1e(mol_m, auxmol, 'int3c1e_ip1', comp=3) + dip1 = (ip1_p - ip1_m) / (2 * delta) + + for alpha in range(3): + comp_idx = alpha * 3 + beta + ana_block = analytic[comp_idx][i_mask][:, j0:j1, :] + fd_block = -dip1[alpha][i_mask][:, j0:j1, :] + err = np.max(np.abs(ana_block - fd_block)) + max_err = max(max_err, err) + + print(f"int3c1e_ipvip1: max FD error = {max_err:.2e}") + assert max_err < 1e-7, f"int3c1e_ipvip1 FD error too large: {max_err}" + + +def test_int3c1e_ipip2(): + """ + Verify int3c1e_ipip2 (d²/dR₃²) via double FD of int3c1e. + + Since aux atoms are separate from mol, no contamination. + Sign: ipip2 = (nabla_k)^2 S. FD gives d²S/dR_C^2. + d²S/dR_C,α dR_C,β = (-d/dr_k,α)(-d/dr_k,β) S = ipip2[α*3+β] + So FD should directly match ipip2. + """ + mol = _make_mol() + auxmol = _make_auxmol() + delta = 1e-4 + + analytic = _compute_3c1e(mol, auxmol, 'int3c1e_ipip2', comp=9) + nao = mol.nao + naux = auxmol.nao + assert analytic.shape == (9, nao, nao, naux) + + max_err = 0.0 + aoslice_aux = auxmol.aoslice_by_atom() + + for katm in range(auxmol.natm): + k0, k1 = int(aoslice_aux[katm, 2]), int(aoslice_aux[katm, 3]) + for alpha in range(3): + for beta in range(3): + coords_pp = auxmol.atom_coords(unit='Bohr').copy() + coords_pm = auxmol.atom_coords(unit='Bohr').copy() + coords_mp = auxmol.atom_coords(unit='Bohr').copy() + coords_mm = auxmol.atom_coords(unit='Bohr').copy() + coords_pp[katm, alpha] += delta + coords_pp[katm, beta] += delta + coords_pm[katm, alpha] += delta + coords_pm[katm, beta] -= delta + coords_mp[katm, alpha] -= delta + coords_mp[katm, beta] += delta + coords_mm[katm, alpha] -= delta + coords_mm[katm, beta] -= delta + + auxmol_pp = auxmol.set_geom_(coords_pp, unit='Bohr', inplace=False) + auxmol_pm = auxmol.set_geom_(coords_pm, unit='Bohr', inplace=False) + auxmol_mp = auxmol.set_geom_(coords_mp, unit='Bohr', inplace=False) + auxmol_mm = auxmol.set_geom_(coords_mm, unit='Bohr', inplace=False) + + s_pp = _compute_3c1e(mol, auxmol_pp, 'int3c1e', comp=1) + s_pm = _compute_3c1e(mol, auxmol_pm, 'int3c1e', comp=1) + s_mp = _compute_3c1e(mol, auxmol_mp, 'int3c1e', comp=1) + s_mm = _compute_3c1e(mol, auxmol_mm, 'int3c1e', comp=1) + + d2s = (s_pp - s_pm - s_mp + s_mm) / (4 * delta**2) + + comp_idx = alpha * 3 + beta + err = np.max(np.abs( + analytic[comp_idx, :, :, k0:k1] - d2s[:, :, k0:k1] + )) + max_err = max(max_err, err) + + print(f"int3c1e_ipip2: max FD error = {max_err:.2e}") + assert max_err < 1e-6, f"int3c1e_ipip2 FD error too large: {max_err}" + + +def test_translational_invariance(): + """ + Verify: ipip1 + ipvip1 + ip1ip2 = 0 + + From d/dR₁ (ip1) + d/dR₂ (ip1) + d/dR₃ (ip1) = 0 + where each d/dR introduces a sign, but all three nabla operators + are electron-coordinate derivatives, so the relation holds directly. + """ + mol = _make_mol() + auxmol = _make_auxmol() + + ipip1 = _compute_3c1e(mol, auxmol, 'int3c1e_ipip1', comp=9) + ip1ip2 = _compute_3c1e(mol, auxmol, 'int3c1e_ip1ip2', comp=9) + ipvip1 = _compute_3c1e(mol, auxmol, 'int3c1e_ipvip1', comp=9) + + total = ipip1 + ipvip1 + ip1ip2 + max_err = np.max(np.abs(total)) + print(f"Translational invariance (ipip1+ipvip1+ip1ip2=0): max error = {max_err:.2e}") + assert max_err < 1e-10, \ + f"Translational invariance violated: {max_err}" + + +def test_shapes(): + """Verify all new integrals return correct shapes.""" + mol = _make_mol() + auxmol = _make_auxmol() + nao = mol.nao + naux = auxmol.nao + + for name in ['int3c1e_ipip1', 'int3c1e_ip1ip2', + 'int3c1e_ipvip1', 'int3c1e_ipip2']: + result = _compute_3c1e(mol, auxmol, name, comp=9) + expected = (9, nao, nao, naux) + assert result.shape == expected, \ + f"{name}: shape {result.shape} != expected {expected}" + print(f"{name}: shape OK {result.shape}") + + +def test_ipip1_symmetry(): + """Verify d²/dR_A,α dR_A,β == d²/dR_A,β dR_A,α for ipip1.""" + mol = _make_mol() + auxmol = _make_auxmol() + + result = _compute_3c1e(mol, auxmol, 'int3c1e_ipip1', comp=9) + # comp = alpha*3+beta + for alpha in range(3): + for beta in range(alpha + 1, 3): + c_ab = result[alpha * 3 + beta] + c_ba = result[beta * 3 + alpha] + max_diff = np.max(np.abs(c_ab - c_ba)) + print(f"ipip1 symmetry ({alpha},{beta}): max diff = {max_diff:.2e}") + assert max_diff < 1e-12, \ + f"ipip1 not symmetric for ({alpha},{beta}): {max_diff}" + + +def test_ipip2_symmetry(): + """Verify d²/dR_C,α dR_C,β == d²/dR_C,β dR_C,α for ipip2.""" + mol = _make_mol() + auxmol = _make_auxmol() + + result = _compute_3c1e(mol, auxmol, 'int3c1e_ipip2', comp=9) + for alpha in range(3): + for beta in range(alpha + 1, 3): + c_ab = result[alpha * 3 + beta] + c_ba = result[beta * 3 + alpha] + max_diff = np.max(np.abs(c_ab - c_ba)) + print(f"ipip2 symmetry ({alpha},{beta}): max diff = {max_diff:.2e}") + assert max_diff < 1e-12, \ + f"ipip2 not symmetric for ({alpha},{beta}): {max_diff}" + + +if __name__ == '__main__': + print("\n=== Shape tests ===") + test_shapes() + + print("\n=== Symmetry tests ===") + test_ipip1_symmetry() + test_ipip2_symmetry() + + print("\n=== Translational invariance ===") + test_translational_invariance() + + print("\n=== Finite-difference tests ===") + test_int3c1e_ipip1() + test_int3c1e_ip1ip2() + test_int3c1e_ipvip1() + test_int3c1e_ipip2() + + print("\n=== All tests passed ===")