From 75963cef2789c4c59e6047db4db40e2331e4c88b Mon Sep 17 00:00:00 2001 From: "henry.wang1" Date: Thu, 24 Sep 2026 11:53:58 +0800 Subject: [PATCH 1/6] Support stress for Becke weight, support spherical basis for PBCeval_gto_strain_tensor kernel --- gpu4pyscf/lib/pbc/nr_eval_gto.cu | 461 +++++++++++++++--- gpu4pyscf/pbc/dft/gen_grid.py | 5 +- gpu4pyscf/pbc/dft/tests/test_pbc_grids.py | 23 +- gpu4pyscf/pbc/grad/krks.py | 1 + gpu4pyscf/pbc/grad/krks_stress.py | 4 +- gpu4pyscf/pbc/grad/kuks.py | 1 + .../pbc/grad/tests/test_pbc_krks_stress.py | 58 ++- .../pbc/grad/tests/test_pbc_rks_stress.py | 42 ++ 8 files changed, 505 insertions(+), 90 deletions(-) diff --git a/gpu4pyscf/lib/pbc/nr_eval_gto.cu b/gpu4pyscf/lib/pbc/nr_eval_gto.cu index 5eecda901..9c9c90cac 100644 --- a/gpu4pyscf/lib/pbc/nr_eval_gto.cu +++ b/gpu4pyscf/lib/pbc/nr_eval_gto.cu @@ -151,8 +151,8 @@ void _cart_deriv1_strain_tensor( } } -template __device__ -void _eval_cart_deriv1_strain_tensor( +template __device__ +void _eval_deriv1_strain_tensor( double *out, double *img_coords, double *env, double xi, double yi, double zi, double rrcutoff, int *bas, int nimgs, int nao, int ngrids) @@ -198,12 +198,203 @@ void _eval_cart_deriv1_strain_tensor( ao, gx, gy, gz, a2, rx, ry, rz, Rx, Ry, Rz, n0); } } - for (int n = 0; n < min(n_cart-n0, 2); ++n) { - for (int x = 0; x < 9; ++x) { - for (int s = 0; s < 4; ++s) { - out[(x*4+s)*naog+(n0+n)*ngrids] = ao[n*36+s*9+x]; + if constexpr (cart) { + for (int n = 0; n < min(n_cart-n0, 2); ++n) { + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+(n0+n)*ngrids] = ao[n*36+s*9+x]; + } } } + } else { + switch (ANG) { + case 0: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog] = ao[s*9+x]; + } + } + break; + case 1: + switch (n0) { + case 0: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+0*ngrids] = ao[0*36+s*9+x]; + out[(x*4+s)*naog+1*ngrids] = ao[1*36+s*9+x]; + } + } + break; + case 2: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+2*ngrids] = ao[0*36+s*9+x]; + } + } + break; + } + break; + case 2: + switch (n0) { + case 0: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+2*ngrids] = -0.315391565252520002 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+4*ngrids] = 0.546274215296039535 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+0*ngrids] = 1.092548430592079070 * ao[1*36+s*9+x]; + } + } + break; + case 2: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+3*ngrids] = 1.092548430592079070 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+2*ngrids]+= -0.315391565252520002 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+4*ngrids]+= -0.546274215296039535 * ao[1*36+s*9+x]; + } + } + break; + case 4: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+1*ngrids] = 1.092548430592079070 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+2*ngrids]+= 0.630783130505040012 * ao[1*36+s*9+x]; + } + } + break; + } + break; + case 3: + switch (n0) { + case 0: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+4*ngrids] = -0.457045799464465739 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+6*ngrids] = 0.590043589926643510 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+0*ngrids] = 1.770130769779930531 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+2*ngrids] = -0.457045799464465739 * ao[1*36+s*9+x]; + } + } + break; + case 2: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+3*ngrids] = -1.119528997770346170 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+5*ngrids] = 1.445305721320277020 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+4*ngrids]+= -0.457045799464465739 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+6*ngrids]+= -1.770130769779930530 * ao[1*36+s*9+x]; + } + } + break; + case 4: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+1*ngrids] = 2.890611442640554055 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+4*ngrids]+= 1.828183197857862944 * ao[1*36+s*9+x]; + } + } + break; + case 6: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+0*ngrids]+= -0.590043589926643510 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+2*ngrids]+= -0.457045799464465739 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+3*ngrids]+= -1.119528997770346170 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+5*ngrids]+= -1.445305721320277020 * ao[1*36+s*9+x]; + } + } + break; + case 8: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+2*ngrids]+= 1.828183197857862944 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+3*ngrids]+= 0.746352665180230782 * ao[1*36+s*9+x]; + } + } + break; + } + break; + case 4: + switch (n0) { + case 0: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+4*ngrids] = 0.317356640745612911 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+6*ngrids] = -0.473087347878780009 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+8*ngrids] = 0.625835735449176134 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+0*ngrids] = 2.503342941796704538 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+2*ngrids] = -0.946174695757560014 * ao[1*36+s*9+x]; + } + } + break; + case 2: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+5*ngrids] = -2.007139630671867500 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+7*ngrids] = 1.770130769779930531 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+4*ngrids]+= 0.634713281491225822 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+8*ngrids]+= -3.755014412695056800 * ao[1*36+s*9+x]; + } + } + break; + case 4: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+1*ngrids] = 5.310392309339791593 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+3*ngrids] = -2.007139630671867500 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+4*ngrids]+= -2.538853125964903290 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+6*ngrids]+= 2.838524087272680054 * ao[1*36+s*9+x]; + } + } + break; + case 6: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+0*ngrids]+= -2.503342941796704538 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+2*ngrids]+= -0.946174695757560014 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+5*ngrids]+= -2.007139630671867500 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+7*ngrids]+= -5.310392309339791590 * ao[1*36+s*9+x]; + } + } + break; + case 8: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+2*ngrids]+= 5.677048174545360108 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+5*ngrids]+= 2.676186174229156671 * ao[1*36+s*9+x]; + } + } + break; + case 10: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+4*ngrids]+= 0.317356640745612911 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+6*ngrids]+= 0.473087347878780009 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+8*ngrids]+= 0.625835735449176134 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+1*ngrids]+= -1.770130769779930530 * ao[1*36+s*9+x]; + out[(x*4+s)*naog+3*ngrids]+= -2.007139630671867500 * ao[1*36+s*9+x]; + } + } + break; + case 12: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+4*ngrids]+= -2.538853125964903290 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+6*ngrids]+= -2.838524087272680054 * ao[0*36+s*9+x]; + out[(x*4+s)*naog+3*ngrids]+= 2.676186174229156671 * ao[1*36+s*9+x]; + } + } + break; + case 14: + for (int x = 0; x < 9; ++x) { + for (int s = 0; s < 4; ++s) { + out[(x*4+s)*naog+4*ngrids]+= 0.846284375321634430 * ao[0*36+s*9+x]; + } + } + break; + } + break; + } } } } @@ -1158,8 +1349,9 @@ static void _sph_ip2_kernel(double *out, PBCIntEnvVars envs, double *grids, } } +template __global__ -static void _cart_deriv0_strain_tensor_kernel( +static void _deriv0_strain_tensor_kernel( double *out, PBCIntEnvVars envs, double *grids, size_t ngrids, int nao, double *rcut) { @@ -1375,83 +1567,186 @@ static void _cart_deriv0_strain_tensor_kernel( } } } } - if (li < 4) { - int *ao_loc = envs.ao_loc; + const int *ao_loc = envs.ao_loc; + out += (size_t)ao_loc[bas_id] * ngrids + grid_id; + const size_t naog = nao * ngrids; + if constexpr (cart) { int nf = (li + 1) * (li + 2) / 2; - out += (size_t)ao_loc[bas_id] * ngrids + grid_id; - size_t naog = nao * ngrids; - for (int n = 0; n < 10; ++n) { - if (n >= nf) break; - for (int x = 0; x < 9; ++x) { - out[x*naog+n*ngrids] = ao[n*9+x]; + if (li < 4) { + for (int n = 0; n < 10; ++n) { + if (n >= nf) break; + for (int x = 0; x < 9; ++x) { + out[x*naog+n*ngrids] = ao[n*9+x]; + } + } + } else { + for (int n = 0; n < 15; ++n) { + if (n >= nf) break; + for (int x = 0; x < 6; ++x) { + out[x*naog+n*ngrids] = ao[n*6+x]; + } + } + out += 6 * naog; // To process zx, zy, zz + + for (int n = 0; n < 45; ++n) { + ao[n] = 0; + } + for (int img = 0; img < nimgs; ++img) { + double ce = 0; + double ce_2a = 0; + double Rx = img_coords[img*3+0] + cell0_Rx; + double Ry = img_coords[img*3+1] + cell0_Ry; + double Rz = img_coords[img*3+2] + cell0_Rz; + double rx = xi - Rx; + double ry = yi - Ry; + double rz = zi - Rz; + double rr = rx * rx + ry * ry + rz * rz; + if (rr > rrcutoff) continue; + for (int ip = 0; ip < nprim; ++ip) { + double ai = expi[ip]; + double c_exp = ci[ip] * exp(-ai * rr); + ce += c_exp; + ce_2a -= c_exp * ai * 2; + } + if (fabs(ce) < 1e-18) continue; + double az = ce_2a * rz; + gto[0 ] = az * rx * rx * rx * rx; + gto[1 ] = az * rx * rx * rx * ry; + gto[2 ] = (az * rz + ce) * rx * rx * rx; + gto[3 ] = az * rx * rx * ry * ry; + gto[4 ] = (az * rz + ce) * rx * rx * ry; + gto[5 ] = (az * rz + 2 * ce) * rx * rx * rz; + gto[6 ] = az * rx * ry * ry * ry; + gto[7 ] = (az * rz + ce) * rx * ry * ry; + gto[8 ] = (az * rz + 2 * ce) * rx * ry * rz; + gto[9 ] = (az * rz + 3 * ce) * rx * rz * rz; + gto[10] = az * ry * ry * ry * ry; + gto[11] = (az * rz + ce) * ry * ry * ry; + gto[12] = (az * rz + 2 * ce) * ry * ry * rz; + gto[13] = (az * rz + 3 * ce) * ry * rz * rz; + gto[14] = (az * rz + 4 * ce) * rz * rz * rz; + for (int n = 0; n < 15; n++) { + ao[0+3*n] -= gto[n] * Rx; + ao[1+3*n] -= gto[n] * Ry; + ao[2+3*n] -= gto[n] * Rz; + } + } + for (int n = 0; n < 15; ++n) { + for (int x = 0; x < 3; ++x) { + out[x*naog+n*ngrids] = ao[n*3+x]; + } } } } else { - int *ao_loc = envs.ao_loc; - int nf = (li + 1) * (li + 2) / 2; - out += (size_t)ao_loc[bas_id] * ngrids + grid_id; - size_t naog = nao * ngrids; - for (int n = 0; n < 15; ++n) { - if (n >= nf) break; - for (int x = 0; x < 6; ++x) { - out[x*naog+n*ngrids] = ao[n*6+x]; + switch (li) { + case 0: + for (int n = 0; n < 9; ++n) { + out[n * naog] = ao[n]; } - } - out += 6 * naog; // To process zx, zy, zz + break; + case 1: + for (int n = 0; n < 9; ++n) { + out[n*naog+0*ngrids] = ao[0*9+n]; + out[n*naog+1*ngrids] = ao[1*9+n]; + out[n*naog+2*ngrids] = ao[2*9+n]; + } + break; + case 2: + for (int n = 0; n < 9; ++n) { + out[n*naog+0*ngrids] = 1.092548430592079070 * ao[1*9+n]; + out[n*naog+1*ngrids] = 1.092548430592079070 * ao[4*9+n]; + out[n*naog+2*ngrids] = 0.630783130505040012 * ao[5*9+n] - 0.315391565252520002 * (ao[0*9+n] + ao[3*9+n]); + out[n*naog+3*ngrids] = 1.092548430592079070 * ao[2*9+n]; + out[n*naog+4*ngrids] = 0.546274215296039535 * (ao[0*9+n] - ao[3*9+n]); + } + break; + case 3: + for (int n = 0; n < 9; ++n) { + out[n*naog+0*ngrids] = 1.770130769779930531 * ao[1*9+n] - 0.590043589926643510 * ao[6*9+n]; + out[n*naog+1*ngrids] = 2.890611442640554055 * ao[4*9+n]; + out[n*naog+2*ngrids] = 1.828183197857862944 * ao[8*9+n] - 0.457045799464465739 * (ao[1*9+n] + ao[6*9+n]); + out[n*naog+3*ngrids] = 0.746352665180230782 * ao[9*9+n] - 1.119528997770346170 * (ao[2*9+n] + ao[7*9+n]); + out[n*naog+4*ngrids] = 1.828183197857862944 * ao[5*9+n] - 0.457045799464465739 * (ao[0*9+n] + ao[3*9+n]); + out[n*naog+5*ngrids] = 1.445305721320277020 * (ao[2*9+n] - ao[7*9+n]); + out[n*naog+6*ngrids] = 0.590043589926643510 * ao[0*9+n] - 1.770130769779930530 * ao[3*9+n]; + } + break; + case 4: + for (int n = 0; n < 6; ++n) { + out[n*naog+0*ngrids] = 2.503342941796704538 * (ao[1*6+n] - ao[6*6+n]) ; + out[n*naog+1*ngrids] = 5.310392309339791593 * ao[4*6+n] - 1.770130769779930530 * ao[11*6+n]; + out[n*naog+2*ngrids] = 5.677048174545360108 * ao[8*6+n] - 0.946174695757560014 * (ao[1*6+n] + ao[6*6+n]); + out[n*naog+3*ngrids] = 2.676186174229156671 * ao[13*6+n]- 2.007139630671867500 * (ao[4*6+n] + ao[11*6+n]); + out[n*naog+4*ngrids] = 0.317356640745612911 * (ao[0*6+n] + ao[10*6+n]) + 0.634713281491225822 * ao[3*6+n] - 2.538853125964903290 * (ao[5*6+n] + ao[12*6+n]) + 0.846284375321634430 * ao[14*6+n]; + out[n*naog+5*ngrids] = 2.676186174229156671 * ao[9*6+n] - 2.007139630671867500 * (ao[2*6+n] + ao[7*6+n]); + out[n*naog+6*ngrids] = 2.838524087272680054 * (ao[5*6+n] - ao[12*6+n]) + 0.473087347878780009 * (ao[10*6+n]- ao[0*6+n]); + out[n*naog+7*ngrids] = 1.770130769779930531 * ao[2*6+n] - 5.310392309339791590 * ao[7*6+n]; + out[n*naog+8*ngrids] = 0.625835735449176134 * (ao[0*6+n] + ao[10*6+n]) - 3.755014412695056800 * ao[3*6+n]; + } + out += 6 * naog; // To process zx, zy, zz - for (int n = 0; n < 45; ++n) { - ao[n] = 0; - } - for (int img = 0; img < nimgs; ++img) { - double ce = 0; - double ce_2a = 0; - double Rx = img_coords[img*3+0] + cell0_Rx; - double Ry = img_coords[img*3+1] + cell0_Ry; - double Rz = img_coords[img*3+2] + cell0_Rz; - double rx = xi - Rx; - double ry = yi - Ry; - double rz = zi - Rz; - double rr = rx * rx + ry * ry + rz * rz; - if (rr > rrcutoff) continue; - for (int ip = 0; ip < nprim; ++ip) { - double ai = expi[ip]; - double c_exp = ci[ip] * exp(-ai * rr); - ce += c_exp; - ce_2a -= c_exp * ai * 2; + for (int n = 0; n < 45; ++n) { + ao[n] = 0; } - if (fabs(ce) < 1e-18) continue; - double az = ce_2a * rz; - gto[0 ] = az * rx * rx * rx * rx; - gto[1 ] = az * rx * rx * rx * ry; - gto[2 ] = (az * rz + ce) * rx * rx * rx; - gto[3 ] = az * rx * rx * ry * ry; - gto[4 ] = (az * rz + ce) * rx * rx * ry; - gto[5 ] = (az * rz + 2 * ce) * rx * rx * rz; - gto[6 ] = az * rx * ry * ry * ry; - gto[7 ] = (az * rz + ce) * rx * ry * ry; - gto[8 ] = (az * rz + 2 * ce) * rx * ry * rz; - gto[9 ] = (az * rz + 3 * ce) * rx * rz * rz; - gto[10] = az * ry * ry * ry * ry; - gto[11] = (az * rz + ce) * ry * ry * ry; - gto[12] = (az * rz + 2 * ce) * ry * ry * rz; - gto[13] = (az * rz + 3 * ce) * ry * rz * rz; - gto[14] = (az * rz + 4 * ce) * rz * rz * rz; - for (int n = 0; n < 15; n++) { - ao[0+3*n] -= gto[n] * Rx; - ao[1+3*n] -= gto[n] * Ry; - ao[2+3*n] -= gto[n] * Rz; + for (int img = 0; img < nimgs; ++img) { + double ce = 0; + double ce_2a = 0; + double Rx = img_coords[img*3+0] + cell0_Rx; + double Ry = img_coords[img*3+1] + cell0_Ry; + double Rz = img_coords[img*3+2] + cell0_Rz; + double rx = xi - Rx; + double ry = yi - Ry; + double rz = zi - Rz; + double rr = rx * rx + ry * ry + rz * rz; + if (rr > rrcutoff) continue; + for (int ip = 0; ip < nprim; ++ip) { + double ai = expi[ip]; + double c_exp = ci[ip] * exp(-ai * rr); + ce += c_exp; + ce_2a -= c_exp * ai * 2; + } + if (fabs(ce) < 1e-18) continue; + double az = ce_2a * rz; + gto[0 ] = az * rx * rx * rx * rx; + gto[1 ] = az * rx * rx * rx * ry; + gto[2 ] = (az * rz + ce) * rx * rx * rx; + gto[3 ] = az * rx * rx * ry * ry; + gto[4 ] = (az * rz + ce) * rx * rx * ry; + gto[5 ] = (az * rz + 2 * ce) * rx * rx * rz; + gto[6 ] = az * rx * ry * ry * ry; + gto[7 ] = (az * rz + ce) * rx * ry * ry; + gto[8 ] = (az * rz + 2 * ce) * rx * ry * rz; + gto[9 ] = (az * rz + 3 * ce) * rx * rz * rz; + gto[10] = az * ry * ry * ry * ry; + gto[11] = (az * rz + ce) * ry * ry * ry; + gto[12] = (az * rz + 2 * ce) * ry * ry * rz; + gto[13] = (az * rz + 3 * ce) * ry * rz * rz; + gto[14] = (az * rz + 4 * ce) * rz * rz * rz; + for (int n = 0; n < 15; n++) { + ao[0+3*n] -= gto[n] * Rx; + ao[1+3*n] -= gto[n] * Ry; + ao[2+3*n] -= gto[n] * Rz; + } } - } - for (int n = 0; n < 15; ++n) { - for (int x = 0; x < 3; ++x) { - out[x*naog+n*ngrids] = ao[n*3+x]; + for (int n = 0; n < 3; ++n) { + out[n*naog+0*ngrids] = 2.503342941796704538 * (ao[1*3+n] - ao[6*3+n]) ; + out[n*naog+1*ngrids] = 5.310392309339791593 * ao[4*3+n] - 1.770130769779930530 * ao[11*3+n]; + out[n*naog+2*ngrids] = 5.677048174545360108 * ao[8*3+n] - 0.946174695757560014 * (ao[1*3+n] + ao[6*3+n]); + out[n*naog+3*ngrids] = 2.676186174229156671 * ao[13*3+n]- 2.007139630671867500 * (ao[4*3+n] + ao[11*3+n]); + out[n*naog+4*ngrids] = 0.317356640745612911 * (ao[0*3+n] + ao[10*3+n]) + 0.634713281491225822 * ao[3*3+n] - 2.538853125964903290 * (ao[5*3+n] + ao[12*3+n]) + 0.846284375321634430 * ao[14*3+n]; + out[n*naog+5*ngrids] = 2.676186174229156671 * ao[9*3+n] - 2.007139630671867500 * (ao[2*3+n] + ao[7*3+n]); + out[n*naog+6*ngrids] = 2.838524087272680054 * (ao[5*3+n] - ao[12*3+n]) + 0.473087347878780009 * (ao[10*3+n]- ao[0*3+n]); + out[n*naog+7*ngrids] = 1.770130769779930531 * ao[2*3+n] - 5.310392309339791590 * ao[7*3+n]; + out[n*naog+8*ngrids] = 0.625835735449176134 * (ao[0*3+n] + ao[10*3+n]) - 3.755014412695056800 * ao[3*3+n]; } + break; } } } +template __global__ -static void _cart_deriv1_strain_tensor_kernel( +static void _deriv1_strain_tensor_kernel( double *out, PBCIntEnvVars envs, double *grids, size_t ngrids, int nao, double *rcut) { @@ -1477,19 +1772,19 @@ static void _cart_deriv1_strain_tensor_kernel( int nimgs = envs.nimgs; switch (li) { - case 0: _eval_cart_deriv1_strain_tensor<0>(out, img_coords, env, + case 0: _eval_deriv1_strain_tensor<0, cart>(out, img_coords, env, xi, yi, zi, rrcutoff, bas, nimgs, nao, ngrids); break; - case 1: _eval_cart_deriv1_strain_tensor<1>(out, img_coords, env, + case 1: _eval_deriv1_strain_tensor<1, cart>(out, img_coords, env, xi, yi, zi, rrcutoff, bas, nimgs, nao, ngrids); break; - case 2: _eval_cart_deriv1_strain_tensor<2>(out, img_coords, env, + case 2: _eval_deriv1_strain_tensor<2, cart>(out, img_coords, env, xi, yi, zi, rrcutoff, bas, nimgs, nao, ngrids); break; - case 3: _eval_cart_deriv1_strain_tensor<3>(out, img_coords, env, + case 3: _eval_deriv1_strain_tensor<3, cart>(out, img_coords, env, xi, yi, zi, rrcutoff, bas, nimgs, nao, ngrids); break; - case 4: _eval_cart_deriv1_strain_tensor<4>(out, img_coords, env, + case 4: _eval_deriv1_strain_tensor<4, cart>(out, img_coords, env, xi, yi, zi, rrcutoff, bas, nimgs, nao, ngrids); break; } @@ -1543,19 +1838,23 @@ int PBCeval_gto_strain_tensor(double *out, PBCIntEnvVars *envs, double *grids, int ngrids, int nao, int nbas, int deriv, int cart, double *rcut) { - if (!cart) { - fprintf(stderr, "PBCeval_gto_strain_tensor does not support spherical GTOs\n"); - return 1; - } constexpr int ngrids_per_block = THREADS; int threads = ngrids_per_block; dim3 blocks((ngrids+ngrids_per_block-1)/ngrids_per_block, nbas); switch (deriv) { case 0: - _cart_deriv0_strain_tensor_kernel<<>>(out, *envs, grids, ngrids, nao, rcut); + if (cart) { + _deriv0_strain_tensor_kernel <<>>(out, *envs, grids, ngrids, nao, rcut); + } else { + _deriv0_strain_tensor_kernel <<>>(out, *envs, grids, ngrids, nao, rcut); + } break; case 1: - _cart_deriv1_strain_tensor_kernel<<>>(out, *envs, grids, ngrids, nao, rcut); + if (cart) { + _deriv1_strain_tensor_kernel <<>>(out, *envs, grids, ngrids, nao, rcut); + } else { + _deriv1_strain_tensor_kernel <<>>(out, *envs, grids, ngrids, nao, rcut); + } break; default: fprintf(stderr, "PBCeval_gto_strain_tensor deriv = %d not supported\n", deriv); diff --git a/gpu4pyscf/pbc/dft/gen_grid.py b/gpu4pyscf/pbc/dft/gen_grid.py index 31e990669..c42a64b90 100644 --- a/gpu4pyscf/pbc/dft/gen_grid.py +++ b/gpu4pyscf/pbc/dft/gen_grid.py @@ -246,7 +246,10 @@ def get_becke_weight_derivative(grids, natm, grid_range = None): dweight_dA_unitcell = cp.zeros([natm, 3, ngrids]) cp.add.at(dweight_dA_unitcell, grids_supatm_to_atm_idx, dweight_dA_supercell) - return dweight_dA_unitcell + weight_stress = cp.einsum("Axg,Ay->xyg", dweight_dA_supercell, grids_supatm_coords) + + dweight_gradient_stress = cp.vstack((dweight_dA_unitcell, weight_stress)) + return dweight_gradient_stress class UniformGrids(lib.StreamObject): '''Uniform Grid class.''' diff --git a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py index 3824de845..6d4ec87f5 100644 --- a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py +++ b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py @@ -27,6 +27,7 @@ from gpu4pyscf.pbc.grad.kuks import get_vxc_full_response as unrestricted_get_vxc_full_response from gpu4pyscf.pbc.grad.kuks import get_vxc as unrestricted_get_vxc from gpu4pyscf.dft.tests.test_grids import find_matching_index_between_two_grids +from gpu4pyscf.pbc.grad.rhf import _finite_diff_cells def numerical_gradient_exc_becke(cell, xc, kpts, auxbasis, atom_grid, dm, unrestricted=False): def get_energy(cell): @@ -191,7 +192,7 @@ def test_becke_weight_derivative(self): test_dw_truncated = get_becke_weight_derivative(grids, cell.natm, truncation_range) dx = 1e-5 - reference_dw = cp.empty([cell.natm, 3, grids.coords.shape[0]]) + reference_dw = cp.empty([cell.natm + 3, 3, grids.coords.shape[0]]) cell_copy = cell.copy() for i_atom in range(cell.natm): for i_xyz in range(3): @@ -213,10 +214,26 @@ def test_becke_weight_derivative(self): reference_dw[i_atom, i_xyz, :] = (w_p - w_m) / (2 * dx) + for i_xyz in range(3): + for j_xyz in range(3): + cell_p, cell_m = _finite_diff_cells(cell, i_xyz, j_xyz, disp = dx) + grids.reset(cell_p) + grids.build() + w_p = grids.weights.copy() + + grids.reset(cell_m) + grids.build() + w_m = grids.weights.copy() + + reference_dw[cell.natm + i_xyz, j_xyz] = (w_p - w_m) / (2 * dx) + reference_dw_truncated = reference_dw[:, :, truncation_range[0] : truncation_range[1]] - assert cp.max(cp.abs(test_dw - reference_dw)) < 2e-9 - assert cp.max(cp.abs(test_dw_truncated - reference_dw_truncated)) < 2e-9 + assert cp.max(cp.abs(test_dw[:-3] - reference_dw[:-3])) < 2e-9 + assert cp.max(cp.abs(test_dw_truncated[:-3] - reference_dw_truncated[:-3])) < 2e-9 + + assert cp.max(cp.abs(test_dw[-3:] - reference_dw[-3:])) < 5e-9 + assert cp.max(cp.abs(test_dw_truncated[-3:] - reference_dw_truncated[-3:])) < 5e-9 def test_xc_gradient_lda_with_response(self): cell = pyscf.M( diff --git a/gpu4pyscf/pbc/grad/krks.py b/gpu4pyscf/pbc/grad/krks.py index 835ff77ef..547b508fe 100644 --- a/gpu4pyscf/pbc/grad/krks.py +++ b/gpu4pyscf/pbc/grad/krks.py @@ -115,6 +115,7 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): else: assert rho.ndim == 1 dweight_dA = get_becke_weight_derivative(grids, natm, (g0,g1)) + dweight_dA = dweight_dA[:-3] de_grid_response_weight += cp.einsum("Adg->Ad", dweight_dA * (rho * exc)) del dweight_dA, rho, exc assert g1 == ngrids diff --git a/gpu4pyscf/pbc/grad/krks_stress.py b/gpu4pyscf/pbc/grad/krks_stress.py index 25cfc4124..5a0e0af34 100644 --- a/gpu4pyscf/pbc/grad/krks_stress.py +++ b/gpu4pyscf/pbc/grad/krks_stress.py @@ -303,8 +303,8 @@ def _eval_ao_strain_derivatives(cell, coords, kpts=None, deriv=0, out=None, coords = cp.asarray(coords.T, order='C') bvk_ncells = opt.bvk_ncells comp = (deriv+1)*(deriv+2)*(deriv+3)//6 - nao = cell.nao_nr(cart=True) - cart = 1 + nao = cell.nao_nr() + cart = cell.cart out = cp.empty((3, 3, comp, bvk_ncells, nao, ngrids)) drv = libpbc.PBCeval_gto_strain_tensor diff --git a/gpu4pyscf/pbc/grad/kuks.py b/gpu4pyscf/pbc/grad/kuks.py index 53be3f39e..a4d1dbb28 100644 --- a/gpu4pyscf/pbc/grad/kuks.py +++ b/gpu4pyscf/pbc/grad/kuks.py @@ -131,6 +131,7 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): else: assert rho.ndim == 1 dweight_dA = get_becke_weight_derivative(grids, natm, (g0,g1)) + dweight_dA = dweight_dA[:-3] de_grid_response_weight += cp.einsum("Adg->Ad", dweight_dA * (rho * exc)) del dweight_dA, rho, exc assert g1 == ngrids diff --git a/gpu4pyscf/pbc/grad/tests/test_pbc_krks_stress.py b/gpu4pyscf/pbc/grad/tests/test_pbc_krks_stress.py index 34c5c65fe..e4de8db50 100644 --- a/gpu4pyscf/pbc/grad/tests/test_pbc_krks_stress.py +++ b/gpu4pyscf/pbc/grad/tests/test_pbc_krks_stress.py @@ -66,14 +66,40 @@ def _check_vs_finite_diff(dat, mf_scanner): assert abs(dat[i,j] - (e1-e2)/2/disp/vol) < 1e-6 class KnownValues(unittest.TestCase): - def test_eval_ao_kpts(self): + def test_eval_ao_cart_kpts(self): a = np.eye(3) * 5 np.random.seed(5) a += np.random.rand(3, 3) - .5 cell = gto.M(atom='He 1 1 1; He 2 1.5 2.4', basis=[[0, [.5, 1]], [1, [1.5, 1], [.5, 1]], - [2, [.8, 1]]], a=a, unit='Bohr', cart=True) + [2, [.8, 1]], + [4, [.8, 1]]], a=a, unit='Bohr', cart=True) + kmesh = [3, 1, 1] + kpts = cell.make_kpts(kmesh) + coords = np.random.rand(10, 3) + ao_value = _eval_ao_strain_derivatives(cell, coords, kpts) + ao_value = ao_value.get().transpose(0,1,2,3,5,4) + for (i, j) in [(0, 0), (0, 1), (0, 2), (2, 0), (2, 2)]: + cell1, cell2 = _finite_diff_cells(cell, i, j, disp=1e-5) + cell1.precision = 1e-10 + cell2.precision = 1e-10 + ao1 = dft.numint.eval_ao_kpts(cell1, coords, cell1.make_kpts(kmesh)) + ao2 = dft.numint.eval_ao_kpts(cell2, coords, cell2.make_kpts(kmesh)) + assert abs(ao_value[0][i,j,0] - (ao1[0] - ao2[0]) / 2e-5).max() < 1e-9 + assert abs(ao_value[1][i,j,0] - (ao1[1] - ao2[1]) / 2e-5).max() < 1e-9 + assert abs(ao_value[2][i,j,0] - (ao1[2] - ao2[2]) / 2e-5).max() < 1e-9 + + def test_eval_ao_sph_kpts(self): + a = np.eye(3) * 5 + np.random.seed(5) + a += np.random.rand(3, 3) - .5 + cell = gto.M(atom='He 1 1 1; He 2 1.5 2.4', + basis=[[0, [.5, 1]], + [1, [1.5, 1], [.5, 1]], + [2, [.8, 1]], + [4, [.8, 1]]], a=a, unit='Bohr') + assert not cell.cart kmesh = [3, 1, 1] kpts = cell.make_kpts(kmesh) coords = np.random.rand(10, 3) @@ -96,7 +122,8 @@ def test_eval_ao_deriv1_cart_kpts(self): cell = gto.M(atom='He 1 1 1; He 2 1.5 2.4', basis=[[0, [.5, 1]], [1, [1.5, 1], [.5, 1]], - [2, [.8, 1]]], a=a, unit='Bohr', cart=True) + [2, [.8, 1]], + [4, [.8, 1]]], a=a, unit='Bohr', cart=True) kmesh = [3, 1, 1] kpts = cell.make_kpts(kmesh) coords = np.random.rand(10, 3) @@ -112,6 +139,31 @@ def test_eval_ao_deriv1_cart_kpts(self): assert abs(ao_value[1][i,j] - (ao1[1] - ao2[1]) / 2e-5).max() < 1e-9 assert abs(ao_value[2][i,j] - (ao1[2] - ao2[2]) / 2e-5).max() < 1e-9 + def test_eval_ao_deriv1_sph_kpts(self): + a = np.eye(3) * 5 + np.random.seed(5) + a += np.random.rand(3, 3) - .5 + cell = gto.M(atom='He 1 1 1; He 2 1.5 2.4', + basis=[[0, [.5, 1]], + [1, [1.5, 1], [.5, 1]], + [2, [.8, 1]], + [4, [.8, 1]]], a=a, unit='Bohr') + assert not cell.cart + kmesh = [3, 1, 1] + kpts = cell.make_kpts(kmesh) + coords = np.random.rand(10, 3) + ao_value = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=1) + ao_value = ao_value.get().transpose(0,1,2,3,5,4) + for (i, j) in [(0, 0), (0, 1), (0, 2), (2, 0), (2, 2)]: + cell1, cell2 = _finite_diff_cells(cell, i, j, disp=1e-5) + cell1.precision = 1e-10 + cell2.precision = 1e-10 + ao1 = dft.numint.eval_ao_kpts(cell1, coords, cell1.make_kpts(kmesh), deriv=1) + ao2 = dft.numint.eval_ao_kpts(cell2, coords, cell2.make_kpts(kmesh), deriv=1) + assert abs(ao_value[0][i,j] - (ao1[0] - ao2[0]) / 2e-5).max() < 2e-9 + assert abs(ao_value[1][i,j] - (ao1[1] - ao2[1]) / 2e-5).max() < 2e-9 + assert abs(ao_value[2][i,j] - (ao1[2] - ao2[2]) / 2e-5).max() < 2e-9 + def test_get_vxc_lda(self): a = np.eye(3) * 5 np.random.seed(5) diff --git a/gpu4pyscf/pbc/grad/tests/test_pbc_rks_stress.py b/gpu4pyscf/pbc/grad/tests/test_pbc_rks_stress.py index fb72d43f2..7173d5225 100644 --- a/gpu4pyscf/pbc/grad/tests/test_pbc_rks_stress.py +++ b/gpu4pyscf/pbc/grad/tests/test_pbc_rks_stress.py @@ -94,6 +94,27 @@ def test_eval_ao_cart(self): ao2 = ni.eval_ao(cell2, coords) assert abs(ao_value[i,j,0] - (ao1 - ao2) / 2e-5).max() < 5e-9 + def test_eval_ao_sph(self): + a = np.eye(3) * 5 + np.random.seed(5) + a += np.random.rand(3, 3) - .5 + cell = gto.M(atom='He 1 1 1; He 2 1.5 2.4', + basis=[[0, [.5, 1]], + [1, [1.5, 1], [.5, 1]], + [2, [.8, 1]], + [3, [.7, 1]], + [4, [.6, 1]]], a=a, unit='Bohr') + assert not cell.cart + coords = np.random.rand(10, 3) + ao_value = _eval_ao_strain_derivatives(cell, coords) + ao_value = ao_value.get().transpose(0,1,2,3,5,4)[0] + ni = NumInt() + for (i, j) in [(0, 0), (0, 1), (0, 2), (2, 0), (2, 2)]: + cell1, cell2 = _finite_diff_cells(cell, i, j, disp=1e-5) + ao1 = ni.eval_ao(cell1, coords) + ao2 = ni.eval_ao(cell2, coords) + assert abs(ao_value[i,j,0] - (ao1 - ao2) / 2e-5).max() < 5e-9 + def test_eval_ao_deriv1_cart(self): a = np.eye(3) * 5 np.random.seed(5) @@ -114,6 +135,27 @@ def test_eval_ao_deriv1_cart(self): ao2 = ni.eval_ao(cell2, coords, deriv=1) assert abs(ao_value[i,j] - (ao1 - ao2) / 2e-4).max() < 1e-7 + def test_eval_ao_deriv1_sph(self): + a = np.eye(3) * 5 + np.random.seed(5) + a += np.random.rand(3, 3) - .5 + cell = gto.M(atom='He 1 1 1; He 2 1.5 2.4', + basis=[[0, [.5, 1]], + [1, [1.5, 1], [.5, 1]], + [2, [.8, 1]], + [3, [.7, 1]], + [4, [.6, 1]]], a=a, unit='Bohr') + assert not cell.cart + coords = np.random.rand(10, 3) + ao_value = _eval_ao_strain_derivatives(cell, coords, deriv=1) + ao_value = ao_value.get().transpose(0,1,2,3,5,4)[0] + ni = NumInt() + for (i, j) in [(0, 0), (0, 1), (0, 2), (2, 0), (2, 2)]: + cell1, cell2 = _finite_diff_cells(cell, i, j, disp=1e-4) + ao1 = ni.eval_ao(cell1, coords, deriv=1) + ao2 = ni.eval_ao(cell2, coords, deriv=1) + assert abs(ao_value[i,j] - (ao1 - ao2) / 2e-4).max() < 1e-7 + def test_get_vxc_lda(self): a = np.eye(3) * 5 np.random.seed(5) From c7e6d0617fb30ba3e95d41ef557aedeb87712251 Mon Sep 17 00:00:00 2001 From: "henry.wang1" Date: Thu, 24 Sep 2026 16:12:10 +0800 Subject: [PATCH 2/6] KRKS LDA and GGA stress supported (with or without grid response) --- gpu4pyscf/pbc/grad/krks.py | 158 ++++++++++++++++++++++++++---- gpu4pyscf/pbc/grad/krks_stress.py | 7 -- 2 files changed, 137 insertions(+), 28 deletions(-) diff --git a/gpu4pyscf/pbc/grad/krks.py b/gpu4pyscf/pbc/grad/krks.py index 547b508fe..96f51acf0 100644 --- a/gpu4pyscf/pbc/grad/krks.py +++ b/gpu4pyscf/pbc/grad/krks.py @@ -27,44 +27,115 @@ from gpu4pyscf.lib.cupy_helper import contract from gpu4pyscf.pbc.dft import multigrid, multigrid_v3, BeckeGrids from gpu4pyscf.pbc.dft.gen_grid import get_becke_weight_derivative +from gpu4pyscf.pbc.dft.numint import _GTOvalOpt +from gpu4pyscf.pbc.grad.krks_stress import _eval_ao_strain_derivatives __all__ = ['Gradients'] +XX, XY, XZ = 4, 5, 6 +YX, YY, YZ = 5, 7, 8 +ZX, ZY, ZZ = 6, 8, 9 + +def get_d2mu_dr2(ao_ks): + assert ao_ks.ndim == 4 + nkpts = ao_ks.shape[0] + ngrids = ao_ks.shape[2] + nao = ao_ks.shape[3] + + d2mu_dr2 = cp.empty([nkpts, 3, 3, ngrids, nao], dtype = ao_ks.dtype) + d2mu_dr2[:,0,0,:,:] = ao_ks[:,XX,:,:] + d2mu_dr2[:,0,1,:,:] = ao_ks[:,XY,:,:] + d2mu_dr2[:,1,0,:,:] = ao_ks[:,XY,:,:] + d2mu_dr2[:,0,2,:,:] = ao_ks[:,XZ,:,:] + d2mu_dr2[:,2,0,:,:] = ao_ks[:,XZ,:,:] + d2mu_dr2[:,1,1,:,:] = ao_ks[:,YY,:,:] + d2mu_dr2[:,1,2,:,:] = ao_ks[:,YZ,:,:] + d2mu_dr2[:,2,1,:,:] = ao_ks[:,YZ,:,:] + d2mu_dr2[:,2,2,:,:] = ao_ks[:,ZZ,:,:] + return d2mu_dr2 + def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): '''derivatives of the Exc per cell''' assert dm_kpts.ndim == 3 xctype = ni._xc_type(xc_code) nao = cell.nao nkpts = len(kpts) - vmat = cp.zeros((nkpts,3,nao,nao), dtype=dm_kpts.dtype) + if xctype == 'LDA': + ao_deriv = 0 + elif xctype == 'GGA': ao_deriv = 1 - for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts, - sort_grids=True): + elif xctype == 'MGGA': + ao_deriv = 1 + else: + raise NotImplementedError(f"Unrecognized xctype = {xctype}") + eval_gto_opt = _GTOvalOpt(cell, kpts, deriv=ao_deriv) + + vmat = cp.zeros((nkpts,3,nao,nao), dtype=dm_kpts.dtype) + de_stress_rho = cp.zeros((3,3)) + exc_sum = 0 + + if xctype == 'LDA': + for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho = ni.eval_rho(cell, ao_ks[:,0], dm_kpts, xctype=xctype, hermi=hermi) - vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=0)[1] + exc, vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=0)[:2] wv = weight * vxc[0] aow = cp.einsum('kpi,p->kpi', ao_ks[:,0], wv) for kn in range(nkpts): vmat[kn] += _d1_dot_(ao_ks[kn,1:4], aow[kn]) + del aow, vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain = ao_ks_strain[:,:,:,0] + ao_ks_strain += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) + dm_nu = contract('kpq,kgq->kpg', dm_kpts, ao_ks[:,0].conj()) + drho_stress = contract('kxypg,kpg->xyg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('xyg,g->xy', drho_stress, wv).real + + del ao_ks_strain, dm_nu, drho_stress, wv + + exc_sum += cp.sum(weight * (rho * exc)) + + del rho, exc + elif xctype == 'GGA': - ao_deriv = 2 - for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts, - sort_grids=True): + for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts, xctype=xctype, hermi=hermi) - vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=0)[1] + exc, vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=0)[:2] wv = weight * vxc wv[0] *= .5 for kn in range(nkpts): vmat[kn] += _gga_grad_sum_(ao_ks[kn], wv) + del vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, coords.T) + del d2ao + + wv[0] *= 2 + dm_nu = contract('kpq,kgq->kpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,kpg->xydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('xydg,dg->xy', dmu_stress, wv).real + del dmu_stress, dm_nu + dm_dmu = contract('kpq,kdgp->kdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,kdqg->xydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('xydg,dg->xy', dnu_stress, wv[1:4]).real + del dnu_stress, dm_dmu + + del ao_ks_strain, wv + + exc_sum += cp.sum(weight * (rho[0] * exc)) + + del rho, exc + elif xctype == 'MGGA': - ao_deriv = 2 - for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts, - sort_grids=True): + for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts, xctype=xctype, hermi=hermi) - vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=0)[1] + exc, vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=0)[:2] wv = weight * vxc wv[0] *= .5 wv[4] *= .5 # for the factor 1/2 in tau @@ -72,15 +143,23 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): vmat[kn] += _gga_grad_sum_(ao_ks[kn], wv[:4]) vmat[kn] += _tau_grad_dot_(ao_ks[kn], wv[4]) + de_stress_rho *= np.nan + + exc_sum += cp.sum(weight * (rho[0] * exc)) + elif xctype == 'HF': pass elif xctype == 'NLC': raise NotImplementedError("NLC") else: - raise NotImplementedError(xc_code) + raise NotImplementedError(f"Unrecognized xctype = {xctype}") + + de_stress_weight = exc_sum * cp.eye(3) - exc = krhf_grad.contract_h1e_dm(cell, vmat, dm_kpts, hermi=1) - exc *= -1.0 / nkpts + exc = np.zeros((cell.natm + 3, 3)) + exc[:-3] = -krhf_grad.contract_h1e_dm(cell, vmat, dm_kpts, hermi=1) + exc[:-3] *= 1.0 / nkpts + exc[-3:] = (de_stress_rho / nkpts + de_stress_weight).get() return exc def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): @@ -103,8 +182,9 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): ao_deriv = 1 else: raise NotImplementedError(f"Unrecognized xctype = {xctype}") + eval_gto_opt = _GTOvalOpt(cell, kpts, deriv=ao_deriv) - de_grid_response_weight = cp.zeros((natm, 3), dtype=cp.float64) + de_grid_response_weight = cp.zeros((natm + 3, 3), dtype=cp.float64) g1 = 0 for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts): g0, g1 = g1, g1 + weight.size @@ -115,13 +195,13 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): else: assert rho.ndim == 1 dweight_dA = get_becke_weight_derivative(grids, natm, (g0,g1)) - dweight_dA = dweight_dA[:-3] de_grid_response_weight += cp.einsum("Adg->Ad", dweight_dA * (rho * exc)) del dweight_dA, rho, exc assert g1 == ngrids dvmat_orbital_response = cp.zeros((nkpts,3,nao,nao), dtype=dm_kpts.dtype) de_grid_response_rho = cp.zeros((natm, 3), dtype=dm_kpts.dtype) + de_stress_rho = cp.zeros((3,3), dtype=cp.float64) g1 = 0 for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts): @@ -141,7 +221,20 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): dvmat_orbital_response[kn] += vtmp_k de_grid_response_rho[i_atom] += cp.einsum('xij,ji->x', vtmp_k, dm_kpts[kn]) * 2 del vtmp_k - del wv, aow, rho, vxc + + del aow, rho, vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain = ao_ks_strain[:,:,:,0] + associated_supatm_coords = grids.supatm_coords[grids.supatm_idx[g0:g1]] + ao_ks_strain += contract('kxgp,yg->kxypg', ao_ks[:,1:4], associated_supatm_coords.T) + del associated_supatm_coords + + dm_nu = contract('kpq,kgq->kpg', dm_kpts, ao_ks[:,0].conj()) + drho_stress = contract('kxypg,kpg->xyg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('xyg,g->xy', drho_stress, wv).real + + del ao_ks_strain, dm_nu, drho_stress, wv elif xctype == 'GGA': rho = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts, xctype=xctype, hermi=hermi) @@ -154,7 +247,26 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): dvmat_orbital_response[kn] += vtmp_k de_grid_response_rho[i_atom] += cp.einsum('xij,ji->x', vtmp_k, dm_kpts[kn]) * 2 del vtmp_k - del wv, rho, vxc + del rho, vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + associated_supatm_coords = grids.supatm_coords[grids.supatm_idx[g0:g1]] + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], associated_supatm_coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, associated_supatm_coords.T) + del d2ao, associated_supatm_coords + + wv[0] *= 2 + dm_nu = contract('kpq,kgq->kpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,kpg->xydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('xydg,dg->xy', dmu_stress, wv).real + del dmu_stress, dm_nu + dm_dmu = contract('kpq,kdgp->kdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,kdqg->xydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('xydg,dg->xy', dnu_stress, wv[1:4]).real + del dnu_stress, dm_dmu + + del ao_ks_strain, wv elif xctype == 'MGGA': rho = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts, xctype=xctype, hermi=hermi) @@ -170,12 +282,16 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): del vtmp_k del wv, rho, vxc + de_stress_rho *= np.nan + else: raise NotImplementedError(f"Unrecognized xctype = {xctype}") assert g1 == ngrids - exc = de_grid_response_rho.get().real - exc -= krhf_grad.contract_h1e_dm(cell, dvmat_orbital_response, dm_kpts, hermi=1) + exc = np.zeros((cell.natm + 3, 3), dtype=np.float64) + exc[:-3] = de_grid_response_rho.get().real + exc[:-3] -= krhf_grad.contract_h1e_dm(cell, dvmat_orbital_response, dm_kpts, hermi=1) + exc[-3:] = de_stress_rho.get() exc *= 1.0 / nkpts exc += de_grid_response_weight.get() return exc diff --git a/gpu4pyscf/pbc/grad/krks_stress.py b/gpu4pyscf/pbc/grad/krks_stress.py index 5a0e0af34..c94b43c70 100644 --- a/gpu4pyscf/pbc/grad/krks_stress.py +++ b/gpu4pyscf/pbc/grad/krks_stress.py @@ -100,13 +100,6 @@ def get_vxc(ks_grad, cell, dm_kpts, kpts, with_j=False, with_nuc=False): assert kpts.ndim == 2 assert dm_kpts.ndim == 3 - if not cell.cart: - c2s = asarray(cell.cart2sph_coeff()) - dm_kpts = sandwich_dot(dm_kpts, c2s.T) - # Ensure all AOs are evaluated in the Cartesian GTOs as ao_ks strain - # derivatives currently supports Cartesian format only - cell = cell.copy() - cell.cart = True nkpts, nao = dm_kpts.shape[:2] assert nkpts == len(kpts) From 37c1a551c7a7cc981e3cecacc85fb04e0c8d9a9e Mon Sep 17 00:00:00 2001 From: "henry.wang1" Date: Thu, 24 Sep 2026 18:02:35 +0800 Subject: [PATCH 3/6] Stress + tests with RKS/KRKS working --- .../df/tests/test_pbc_hcore_derivatives.py | 6 +- gpu4pyscf/pbc/dft/tests/test_pbc_grids.py | 116 +++++++++++------- gpu4pyscf/pbc/grad/krks.py | 74 ++++++++--- gpu4pyscf/pbc/grad/rks.py | 7 +- 4 files changed, 133 insertions(+), 70 deletions(-) diff --git a/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py b/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py index 8c3c2230f..d1a4c2546 100644 --- a/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py +++ b/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py @@ -221,7 +221,7 @@ def get_energy(cell): assert np.abs(test_energy - ref_energy) < 1e-9 assert np.max(np.abs(test_derivatives[:-3, :] - ref_derivatives[:-3, :])) < 3e-8 - # assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 3e-7 # TODO: Support Becke grid stress tensor + assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 1e-7 dm = mf.make_rdm1() kmesh = np.array([1,1,1]) @@ -336,7 +336,7 @@ def get_energy(cell): assert np.abs(test_energy - ref_energy) < 1e-9 assert np.max(np.abs(test_derivatives[:-3, :] - ref_derivatives[:-3, :])) < 5e-7 - # assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 5e-6 # TODO: Support Becke grid stress tensor + assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 1e-6 dm = mf.make_rdm1() kpts = cell.make_kpts(kmesh) @@ -458,7 +458,7 @@ def get_energy(cell): assert np.abs(test_energy - ref_energy) < 1e-9 assert np.max(np.abs(test_derivatives[:-3, :] - ref_derivatives[:-3, :])) < 1e-6 - # assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 1e-6 # TODO: Support Becke grid stress tensor + assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 1e-6 if __name__ == '__main__': print("Full Tests for PBC GDF Hcore gradient and stress tensor") diff --git a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py index 6d4ec87f5..a376ab350 100644 --- a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py +++ b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py @@ -29,8 +29,10 @@ from gpu4pyscf.dft.tests.test_grids import find_matching_index_between_two_grids from gpu4pyscf.pbc.grad.rhf import _finite_diff_cells -def numerical_gradient_exc_becke(cell, xc, kpts, auxbasis, atom_grid, dm, unrestricted=False): +def numerical_gradient_exc_becke(cell, xc, kmesh, auxbasis, atom_grid, dm, unrestricted=False): + assert np.array(kmesh).shape == (3,) def get_energy(cell): + kpts = cell.make_kpts(kmesh) if unrestricted: mf = kuks.KUKS(cell, xc=xc, kpts=kpts) else: @@ -47,7 +49,7 @@ def get_energy(cell): n, exc, vxc = mf._numint.nr_rks(cell, mf.grids, mf.xc, dm, 0, hermi=1, kpts=kpts, kpts_band=None) return exc - numerical_gradient = np.zeros((cell.natm, 3)) + numerical_gradient = np.zeros((cell.natm + 3, 3)) dx = 1e-4 cell_copy = cell.copy() for i_atom in range(cell.natm): @@ -69,6 +71,14 @@ def get_energy(cell): translation_invariance = np.sum(numerical_gradient, axis=0) assert np.max(np.abs(translation_invariance)) < 1e-8, "Bad numerical gradient" + for i_xyz in range(3): + for j_xyz in range(3): + cell_p, cell_m = _finite_diff_cells(cell, i_xyz, j_xyz, disp = dx) + E_p = get_energy(cell_p) + E_m = get_energy(cell_m) + + numerical_gradient[cell.natm + i_xyz, j_xyz] = (E_p - E_m) / (2 * dx) + # np.set_printoptions(precision=16, suppress=True, linewidth=np.inf) # print(repr(numerical_gradient)) return numerical_gradient @@ -265,16 +275,19 @@ def test_xc_gradient_lda_with_response(self): dm = dm[None,:,:] test_gradient = get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) - # ref_gradient = numerical_gradient_exc_becke(cell, "LDA0", kpts, 'def2-universal-jkfit', (50,194), dm) + # ref_gradient = numerical_gradient_exc_becke(cell, "LDA0", [1,1,1], 'def2-universal-jkfit', (50,194), dm) ref_gradient = np.array([ - [ 0.0000000094857455, 0.0000000095390362, 0.0000000095212727], - [-0.0002930449305438, -0.0002930448950167, -0.0002930448950167], - [-0.0000000263611355, -0.0000000263611355, 0.0000000985167503], - [ 0.0002930621789687, 0.0002930622144959, -0.0002931344589285], - [-0.0000000261479727, 0.000000098570041 , -0.0000000263966626], - [ 0.0002930622322594, -0.0002931345122192, 0.0002930622144959], - [ 0.0000000986766224, -0.0000000262012634, -0.0000000262012634], - [-0.0002931345122192, 0.0002930622144959, 0.0002930622144959], + [ 0.0000000095035091, 0.0000000094857455, 0.0000000095390362], + [-0.0002934405785027, -0.0002934406140298, -0.0002934405607391], + [-0.0000000263788991, -0.0000000263611355, 0.0000000984989867], + [ 0.0002934578446911, 0.0002934579157454, -0.0002935301246509], + [-0.0000000261657362, 0.0000000985878046, -0.0000000263788991], + [ 0.0002934578979819, -0.0002935301424145, 0.0002934578802183], + [ 0.0000000986233317, -0.0000000261834998, -0.0000000262190269], + [-0.0002935301779416, 0.0002934578802183, 0.0002934578979819], + [ 2.1212388574554097, 0.0000000038546943, 0.0000000030730973], + [ 0.0000000035527137, 2.1212388573843555, 0.000000003126388 ], + [ 0.0000000037658765, 0.0000000031974423, 2.121238857029084 ], ]) # It can match down to 1e-9, if the finite difference is computed using the same dm from SCF. @@ -292,7 +305,8 @@ def test_xc_gradient_gga_with_response(self): verbose = 0, ) - kpts = cell.make_kpts((1,2,3)) + kmesh = (1,2,3) + kpts = cell.make_kpts(kmesh) mf = krks.KRKS(cell, xc="HSE06", kpts=kpts).density_fit(auxbasis='def2-universal-jkfit') mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (50,194) @@ -306,13 +320,13 @@ def test_xc_gradient_gga_with_response(self): test_gradient = get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) # dm is not very stable, and numerical gradient is super fast - ref_gradient = numerical_gradient_exc_becke(cell, "HSE06", kpts, 'def2-universal-jkfit', (50,194), dm) + ref_gradient = numerical_gradient_exc_becke(cell, "HSE06", kmesh, 'def2-universal-jkfit', (50,194), dm) assert np.max(np.abs(test_gradient - ref_gradient)) < 1e-7 def test_xc_gradient_gga_without_response(self): cell = pyscf.M( - a = np.eye(3) * 3.5668, + a = np.eye(3) * 3.6668, atom = ''' C 0. 0. 0. C 0.8917 0.9017 0.8917 @@ -340,19 +354,22 @@ def test_xc_gradient_gga_without_response(self): dm = dm[None,:,:] test_gradient = get_vxc(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) - # ref_gradient = numerical_gradient_exc_becke(cell, "PBE", kpts, 'def2-universal-jkfit', (99,590), dm) + # ref_gradient = numerical_gradient_exc_becke(cell, "PBE", [1,1,1], 'def2-universal-jkfit', (99,590), dm) ref_gradient = np.array([ - [ 0.000429164934701 , -0.0022409300370896, 0.0004291650412824], - [-0.0000602672400873, -0.0009193305317012, -0.0000602672400873], - [ 0.0004409704956743, -0.0019386623151263, -0.0004303180389797], - [ 0.0000603361982598, 0.0037129129282221, -0.0000604878991339], - [-0.0004055164026795, -0.0022406373290096, -0.0004055164382066], - [ 0.0000604298833196, 0.0024061113990115, 0.0000604298477924], - [-0.0004303180745069, -0.0019386622795992, 0.0004409704601471], - [-0.0000604880057153, 0.0037129129637492, 0.0000603362693141], + [ 0.0057493621596905, 0.0040870682127547, 0.0057493626570704], + [-0.0347991650784252, -0.0320027493927455, -0.0347991649363166], + [ 0.0246147672555708, 0.0208271039525698, 0.0145351780389547], + [-0.0077134854237215, -0.0061064405798561, -0.014724848114156 ], + [ 0.0216872688696412, 0.0122791700718494, 0.0216872685854241], + [-0.0093490771035931, -0.0138048160280846, -0.0093490768904303], + [ 0.0145351779679004, 0.0208271048052211, 0.0246147670779351], + [-0.014724848647063 , -0.0061064410417089, -0.0077134864895356], + [ 2.679998687966645 , 0.204392930385211 , 0.179311419188366 ], + [ 0.204415101734412 , 2.6794902670701504, 0.204415101734412 ], + [ 0.1793114174120092, 0.2043929328365834, 2.6799986901693273], ]) - assert np.max(np.abs(test_gradient - ref_gradient)) < 3e-4 + assert np.max(np.abs(test_gradient - ref_gradient)) < 6e-4 def test_xc_gradient_mgga_with_response(self): cell = pyscf.M( @@ -364,7 +381,8 @@ def test_xc_gradient_mgga_with_response(self): verbose = 0, ) - kpts = cell.make_kpts((1,1,3)) + kmesh = (1,1,3) + kpts = cell.make_kpts(kmesh) mf = krks.KRKS(cell, xc="r2scan", kpts=kpts).density_fit(auxbasis='def2-universal-jkfit') mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (50,194) @@ -378,9 +396,10 @@ def test_xc_gradient_mgga_with_response(self): test_gradient = get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) # dm is not very stable, and numerical gradient is super fast - ref_gradient = numerical_gradient_exc_becke(cell, "r2scan", kpts, 'def2-universal-jkfit', (50,194), dm) + ref_gradient = numerical_gradient_exc_becke(cell, "r2scan", kmesh, 'def2-universal-jkfit', (50,194), dm) - assert np.max(np.abs(test_gradient - ref_gradient)) < 1e-9 + assert np.max(np.abs(test_gradient[:-3] - ref_gradient[:-3])) < 1e-9 + assert np.max(np.abs(test_gradient[-3:] - ref_gradient[-3:])) < 3e-8 def test_xc_gradient_mgga_without_response(self): cell = pyscf.M( @@ -393,7 +412,8 @@ def test_xc_gradient_mgga_without_response(self): verbose = 0, ) - kpts = cell.make_kpts((1,2,3)) + kmesh = (1,2,3) + kpts = cell.make_kpts(kmesh) mf = krks.KRKS(cell, xc="r2scan0", kpts=kpts).density_fit(auxbasis='def2-universal-jkfit') mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (120,590) @@ -407,9 +427,10 @@ def test_xc_gradient_mgga_without_response(self): test_gradient = get_vxc(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) # dm is not very stable, and numerical gradient is super fast - ref_gradient = numerical_gradient_exc_becke(cell, "r2scan0", kpts, 'def2-universal-jkfit', (120,590), dm) + ref_gradient = numerical_gradient_exc_becke(cell, "r2scan0", kmesh, 'def2-universal-jkfit', (120,590), dm) - assert np.max(np.abs(test_gradient - ref_gradient)) < 1e-4 + assert np.max(np.abs(test_gradient[:-3] - ref_gradient[:-3])) < 1e-4 + assert np.max(np.abs(test_gradient[-3:] - ref_gradient[-3:])) < 2e-3 def test_xc_gradient_unrestricted_no_k_without_response(self): cell = pyscf.M( @@ -434,10 +455,13 @@ def test_xc_gradient_unrestricted_no_k_without_response(self): dm = dm[:,None,:,:] test_gradient = unrestricted_get_vxc(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) - # ref_gradient = numerical_gradient_exc_becke(cell, "HSE06", kpts, 'def2-universal-jkfit', (99,590), dm, unrestricted=True) + # ref_gradient = numerical_gradient_exc_becke(cell, "HSE06", (1,1,1), 'def2-universal-jkfit', (99,590), dm, unrestricted=True) ref_gradient = np.array([ - [ 0.0000210273753964, -0.0175452356021566, 0.0000210258033206], - [-0.0000210273665147, 0.0175452356021566, -0.0000210257944389], + [ 0.0000210218686902, -0.0175452375472673, 0.0000210222950159], + [-0.0000210218686902, 0.0175452375383855, -0.0000210222772523], + [-0.5849730212226234, -0.0000565129454344, -0.0047818453818849], + [-0.0000557453816441, -0.5846170845913434, -0.0000557508350596], + [-0.0047818449822046, -0.0000565183100321, -0.5849730417928356], ]) assert np.max(np.abs(test_gradient - ref_gradient)) < 2e-4 @@ -473,16 +497,19 @@ def test_xc_gradient_unrestricted_no_k_with_response(self): dm = dm[:,None,:,:] test_gradient = unrestricted_get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) - # ref_gradient = numerical_gradient_exc_becke(cell, "r2scan", kpts, 'def2-universal-jkfit', (50,194), dm, unrestricted=True) + # ref_gradient = numerical_gradient_exc_becke(cell, "r2scan", (1,1,1), 'def2-universal-jkfit', (50,194), dm, unrestricted=True) ref_gradient = np.array([ - [ 0.0000003321254383, 0.0000003318678665, 0.0000003322675468], - [-0.0003303815532263, -0.0003303816153988, -0.0003303810025557], - [-0.0000001970335006, -0.0000001969180374, 0.0000003920330727], - [ 0.0003302472517674, 0.0003302459639087, -0.0003304415585603], - [-0.000000197513117 , 0.0000003920863634, -0.0000001963496032], - [ 0.0003302466300426, -0.0003304415407968, 0.0003302464079979], - [ 0.0000003911893032, -0.0000001970335006, -0.0000001977085162], - [-0.0003304409634808, 0.0003302473672306, 0.0003302460438448], + [ 0.000000331477068 , 0.0000003322497832, 0.0000003327649267], + [-0.0003304407325544, -0.0003304415763239, -0.0003304411144711], + [-0.0000001969979735, -0.0000001975841712, 0.0000003919886638], + [ 0.0003303068574212, 0.0003303051521186, -0.0003305007467702], + [-0.0000001974953534, 0.0000003921751812, -0.0000001968025742], + [ 0.0003303066264948, -0.0003305007645338, 0.000330305631735 ], + [ 0.0000003910827218, -0.0000001970779095, -0.0000001967759289], + [-0.0003305006668342, 0.0003303075324368, 0.0003303051876458], + [10.485175710082117 , -0.0000000322586402, -0.0000000312638804], + [-0.0000000314681614, 10.48517571377694 , -0.0000000326494387], + [-0.000000033377745 , -0.0000000317967874, 10.485175704362248 ], ]) # It can match down to 1e-9, if the finite difference is computed using the same dm from SCF. @@ -500,7 +527,8 @@ def test_xc_gradient_unrestricted_k_with_response(self): verbose = 0, ) - kpts = cell.make_kpts((1,3,1)) + kmesh = (1,3,1) + kpts = cell.make_kpts(kmesh) mf = kuks.KUKS(cell, xc="lda", kpts=kpts) mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (40,194) @@ -518,7 +546,7 @@ def test_xc_gradient_unrestricted_k_with_response(self): test_gradient = unrestricted_get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) # dm is not very stable, and numerical gradient is super fast - ref_gradient = numerical_gradient_exc_becke(cell, "lda", kpts, None, (40,194), dm, unrestricted=True) + ref_gradient = numerical_gradient_exc_becke(cell, "lda", kmesh, None, (40,194), dm, unrestricted=True) assert np.max(np.abs(test_gradient - ref_gradient)) < 1e-9 diff --git a/gpu4pyscf/pbc/grad/krks.py b/gpu4pyscf/pbc/grad/krks.py index 96f51acf0..c512d1209 100644 --- a/gpu4pyscf/pbc/grad/krks.py +++ b/gpu4pyscf/pbc/grad/krks.py @@ -86,6 +86,10 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): del aow, vxc + exc_sum += cp.sum(weight * (rho * exc)) + + del rho, exc + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) ao_ks_strain = ao_ks_strain[:,:,:,0] ao_ks_strain += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) @@ -95,10 +99,6 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): del ao_ks_strain, dm_nu, drho_stress, wv - exc_sum += cp.sum(weight * (rho * exc)) - - del rho, exc - elif xctype == 'GGA': for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts, xctype=xctype, hermi=hermi) @@ -110,6 +110,10 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): del vxc + exc_sum += cp.sum(weight * (rho[0] * exc)) + + del rho, exc + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) d2ao = get_d2mu_dr2(ao_ks) @@ -128,10 +132,6 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): del ao_ks_strain, wv - exc_sum += cp.sum(weight * (rho[0] * exc)) - - del rho, exc - elif xctype == 'MGGA': for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts, xctype=xctype, hermi=hermi) @@ -143,10 +143,34 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): vmat[kn] += _gga_grad_sum_(ao_ks[kn], wv[:4]) vmat[kn] += _tau_grad_dot_(ao_ks[kn], wv[4]) - de_stress_rho *= np.nan + del vxc exc_sum += cp.sum(weight * (rho[0] * exc)) + del rho, exc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, coords.T) + del d2ao + + wv[0] *= 2 + dm_nu = contract('kpq,kgq->kpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,kpg->xydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('xydg,dg->xy', dmu_stress, wv[0:4]).real + del dmu_stress, dm_nu + dm_dmu = contract('kpq,kdgp->kdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,kdqg->xydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('xydg,dg->xy', dnu_stress, wv[1:4]).real + del dnu_stress, dm_dmu + dm_dnu = contract('kpq,kdgq->kdpg', dm_kpts, ao_ks[:,1:4].conj()) + tau_stress = contract('kxydpg,kdpg->xyg', ao_ks_strain[:,:,:,1:4], dm_dnu) + de_stress_rho += 2 * contract('xyg,g->xy', tau_stress, wv[4]).real + del tau_stress, dm_dnu + + del ao_ks_strain, wv + elif xctype == 'HF': pass elif xctype == 'NLC': @@ -280,9 +304,30 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): dvmat_orbital_response[kn] += vtmp_k de_grid_response_rho[i_atom] += cp.einsum('xij,ji->x', vtmp_k, dm_kpts[kn]) * 2 del vtmp_k - del wv, rho, vxc + del rho, vxc - de_stress_rho *= np.nan + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + associated_supatm_coords = grids.supatm_coords[grids.supatm_idx[g0:g1]] + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], associated_supatm_coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, associated_supatm_coords.T) + del d2ao, associated_supatm_coords + + wv[0] *= 2 + dm_nu = contract('kpq,kgq->kpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,kpg->xydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('xydg,dg->xy', dmu_stress, wv[0:4]).real + del dmu_stress, dm_nu + dm_dmu = contract('kpq,kdgp->kdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,kdqg->xydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('xydg,dg->xy', dnu_stress, wv[1:4]).real + del dnu_stress, dm_dmu + dm_dnu = contract('kpq,kdgq->kdpg', dm_kpts, ao_ks[:,1:4].conj()) + tau_stress = contract('kxydpg,kdpg->xyg', ao_ks_strain[:,:,:,1:4], dm_dnu) + de_stress_rho += 2 * contract('xyg,g->xy', tau_stress, wv[4]).real + del tau_stress, dm_dnu + + del ao_ks_strain, wv else: raise NotImplementedError(f"Unrecognized xctype = {xctype}") @@ -350,12 +395,7 @@ def energy_ee(self, dm, kpts): fn = get_vxc_full_response else: fn = get_vxc - de[:-3] = fn(ni, cell, grids, xc, dm, kpts) - if isinstance(grids, BeckeGrids): - de[-3:] = np.nan - else: - de[-3:] = multigrid_v3.MultiGridNumInt(cell).energy_strain_gradient( - xc, dm, kpts=kpts, spin=0, with_j=False, with_nuc=False) + de = fn(ni, cell, grids, xc, dm, kpts) t0 = log.timer_debug1('vxc', *t0) if j_factor != 0 or k_sr != 0 or k_lr != 0: diff --git a/gpu4pyscf/pbc/grad/rks.py b/gpu4pyscf/pbc/grad/rks.py index 70521667e..45a7b49a9 100644 --- a/gpu4pyscf/pbc/grad/rks.py +++ b/gpu4pyscf/pbc/grad/rks.py @@ -63,12 +63,7 @@ def energy_ee(self, dm): else: fn = get_vxc cell = self.cell - de[:-3] = fn(ni, cell, grids, xc, dm[None], np.zeros((1, 3))) - if isinstance(grids, BeckeGrids): - de[-3:] = np.nan - else: - de[-3:] = multigrid_v3.MultiGridNumInt(cell).energy_strain_gradient( - xc, dm, spin=0, with_j=False, with_nuc=False) + de = fn(ni, cell, grids, xc, dm[None], np.zeros((1, 3))) if j_factor != 0 or k_sr != 0 or k_lr != 0: de += rhf._get_ejk_derivatives(mf, dm, None, j_factor, omega, k_lr, k_sr) From 4dabfe2c06d3edd4380adfb85bb927475ab35bc7 Mon Sep 17 00:00:00 2001 From: "henry.wang1" Date: Thu, 24 Sep 2026 20:46:23 +0800 Subject: [PATCH 4/6] Support UKS and KUKS --- .../df/tests/test_pbc_hcore_derivatives.py | 2 +- gpu4pyscf/pbc/dft/tests/test_pbc_grids.py | 33 ++- gpu4pyscf/pbc/grad/kuks.py | 193 +++++++++++++++--- gpu4pyscf/pbc/grad/uks.py | 8 +- 4 files changed, 182 insertions(+), 54 deletions(-) diff --git a/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py b/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py index d1a4c2546..b4a05f18f 100644 --- a/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py +++ b/gpu4pyscf/pbc/df/tests/test_pbc_hcore_derivatives.py @@ -458,7 +458,7 @@ def get_energy(cell): assert np.abs(test_energy - ref_energy) < 1e-9 assert np.max(np.abs(test_derivatives[:-3, :] - ref_derivatives[:-3, :])) < 1e-6 - assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 1e-6 + assert np.max(np.abs(test_derivatives[-3:, :] - ref_derivatives[-3:, :])) < 2e-6 if __name__ == '__main__': print("Full Tests for PBC GDF Hcore gradient and stress tensor") diff --git a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py index a376ab350..7c2fe6650 100644 --- a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py +++ b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py @@ -19,7 +19,7 @@ from pyscf.pbc.dft import gen_grid as gen_grid_cpu from gpu4pyscf.pbc.dft import gen_grid from pyscf.pbc.dft import rks as rks_cpu -from gpu4pyscf.pbc.dft import rks +from gpu4pyscf.pbc.dft import rks, uks from pyscf.pbc.dft import krks as krks_cpu from gpu4pyscf.pbc.dft import krks, kuks from gpu4pyscf.pbc.dft.gen_grid import get_becke_weight_derivative @@ -29,7 +29,7 @@ from gpu4pyscf.dft.tests.test_grids import find_matching_index_between_two_grids from gpu4pyscf.pbc.grad.rhf import _finite_diff_cells -def numerical_gradient_exc_becke(cell, xc, kmesh, auxbasis, atom_grid, dm, unrestricted=False): +def numerical_gradient_exc_becke(cell, xc, kmesh, auxbasis, atom_grid, dm, unrestricted=False, dx = 1e-4): assert np.array(kmesh).shape == (3,) def get_energy(cell): kpts = cell.make_kpts(kmesh) @@ -50,7 +50,6 @@ def get_energy(cell): return exc numerical_gradient = np.zeros((cell.natm + 3, 3)) - dx = 1e-4 cell_copy = cell.copy() for i_atom in range(cell.natm): for i_xyz in range(3): @@ -443,7 +442,7 @@ def test_xc_gradient_unrestricted_no_k_without_response(self): ) kpts = cell.make_kpts((1,1,1)) - mf = kuks.KUKS(cell, xc="HSE06", kpts=kpts).density_fit(auxbasis='def2-universal-jkfit') + mf = uks.UKS(cell, xc="HSE06").density_fit(auxbasis='def2-universal-jkfit') mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (99,590) mf.conv_tol = 1e-10 @@ -455,16 +454,17 @@ def test_xc_gradient_unrestricted_no_k_without_response(self): dm = dm[:,None,:,:] test_gradient = unrestricted_get_vxc(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) - # ref_gradient = numerical_gradient_exc_becke(cell, "HSE06", (1,1,1), 'def2-universal-jkfit', (99,590), dm, unrestricted=True) + # ref_gradient = numerical_gradient_exc_becke(cell, "HSE06", (1,1,1), 'def2-universal-jkfit', (99,590), dm, unrestricted=True, dx=1e-5) ref_gradient = np.array([ - [ 0.0000210218686902, -0.0175452375472673, 0.0000210222950159], - [-0.0000210218686902, 0.0175452375383855, -0.0000210222772523], - [-0.5849730212226234, -0.0000565129454344, -0.0047818453818849], - [-0.0000557453816441, -0.5846170845913434, -0.0000557508350596], - [-0.0047818449822046, -0.0000565183100321, -0.5849730417928356], + [ 0.0000210217621088, -0.0175452376183216, 0.000021022206198 ], + [-0.0000210217621088, 0.0175452377071394, -0.000021022206198 ], + [-0.5849730015938803, -0.0000565129276708, -0.0047818454085302], + [-0.0000557453638805, -0.5846170649803639, -0.0000557508705867], + [-0.0047818451420767, -0.0000565183455592, -0.5849730222884375], ]) - assert np.max(np.abs(test_gradient - ref_gradient)) < 2e-4 + assert np.max(np.abs(test_gradient[:-3] - ref_gradient[:-3])) < 2e-4 + assert np.max(np.abs(test_gradient[-3:] - ref_gradient[-3:])) < 5e-4 def test_xc_gradient_unrestricted_no_k_with_response(self): cell = pyscf.M( @@ -527,17 +527,13 @@ def test_xc_gradient_unrestricted_k_with_response(self): verbose = 0, ) - kmesh = (1,3,1) + kmesh = (1,4,1) kpts = cell.make_kpts(kmesh) mf = kuks.KUKS(cell, xc="lda", kpts=kpts) mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (40,194) mf.conv_tol = 1e-10 - # TODO: This is a hack to avoid a OOM issue in get_hcore() function (using multigrid_v3 internally) - hcore = mf.to_cpu().get_hcore() - mf.get_hcore = lambda: cp.asarray(hcore) - mf.kernel() dm = mf.make_rdm1() @@ -546,9 +542,10 @@ def test_xc_gradient_unrestricted_k_with_response(self): test_gradient = unrestricted_get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) # dm is not very stable, and numerical gradient is super fast - ref_gradient = numerical_gradient_exc_becke(cell, "lda", kmesh, None, (40,194), dm, unrestricted=True) + ref_gradient = numerical_gradient_exc_becke(cell, "lda", kmesh, None, (40,194), dm, unrestricted=True, dx=1e-5) - assert np.max(np.abs(test_gradient - ref_gradient)) < 1e-9 + assert np.max(np.abs(test_gradient[:-3] - ref_gradient[:-3])) < 1e-9 + assert np.max(np.abs(test_gradient[-3:] - ref_gradient[-3:])) < 4e-9 if __name__ == '__main__': print("Full Tests for PBC Becke grids") diff --git a/gpu4pyscf/pbc/grad/kuks.py b/gpu4pyscf/pbc/grad/kuks.py index a4d1dbb28..3e2d0268a 100644 --- a/gpu4pyscf/pbc/grad/kuks.py +++ b/gpu4pyscf/pbc/grad/kuks.py @@ -28,6 +28,9 @@ from gpu4pyscf.lib.cupy_helper import contract from gpu4pyscf.pbc.dft import multigrid, multigrid_v3, BeckeGrids from gpu4pyscf.pbc.dft.gen_grid import get_becke_weight_derivative +from gpu4pyscf.pbc.dft.numint import _GTOvalOpt +from gpu4pyscf.pbc.grad.krks_stress import _eval_ao_strain_derivatives +from gpu4pyscf.pbc.grad.krks import get_d2mu_dr2 __all__ = ['Gradients'] @@ -36,15 +39,28 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): xctype = ni._xc_type(xc_code) nao = cell.nao nkpts = len(kpts) - vmat = cp.zeros((2,nkpts,3,nao,nao), dtype=dm_kpts.dtype) + if xctype == 'LDA': + ao_deriv = 0 + elif xctype == 'GGA': + ao_deriv = 1 + elif xctype == 'MGGA': ao_deriv = 1 - for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts, - sort_grids=True): + else: + raise NotImplementedError(f"Unrecognized xctype = {xctype}") + eval_gto_opt = _GTOvalOpt(cell, kpts, deriv=ao_deriv) + + vmat = cp.zeros((2,nkpts,3,nao,nao), dtype=dm_kpts.dtype) + de_stress_rho = cp.zeros((3,3)) + exc_sum = 0 + + if xctype == 'LDA': + for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho_a = ni.eval_rho(cell, ao_ks[:,0], dm_kpts[0], xctype=xctype, hermi=hermi) rho_b = ni.eval_rho(cell, ao_ks[:,0], dm_kpts[1], xctype=xctype, hermi=hermi) rho = cp.stack([rho_a, rho_b], axis=0) - vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=1)[1] + del rho_a, rho_b + exc, vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=1)[:2] wv = weight * vxc[:,0] aowa = cp.einsum('xpi,p->xpi', ao_ks[:,0], wv[0]) aowb = cp.einsum('xpi,p->xpi', ao_ks[:,0], wv[1]) @@ -52,28 +68,65 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): vmat[0,kn] += krks_grad._d1_dot_(ao_ks[kn,1:4], aowa[kn]) vmat[1,kn] += krks_grad._d1_dot_(ao_ks[kn,1:4], aowb[kn]) + del aowa, aowb, vxc + + exc_sum += cp.sum(weight * ((rho[0] + rho[1]) * exc)) + + del rho, exc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain = ao_ks_strain[:,:,:,0] + ao_ks_strain += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) + dm_nu = contract('ukpq,kgq->ukpg', dm_kpts, ao_ks[:,0].conj()) + drho_stress = contract('kxypg,ukpg->uxyg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('uxyg,ug->xy', drho_stress, wv).real + + del ao_ks_strain, dm_nu, drho_stress, wv + elif xctype == 'GGA': - ao_deriv = 2 - for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts, - sort_grids=True): + for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho_a = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts[0], xctype=xctype, hermi=hermi) rho_b = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts[1], xctype=xctype, hermi=hermi) rho = cp.stack([rho_a, rho_b], axis=0) - vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=1)[1] + del rho_a, rho_b + exc, vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=1)[:2] wv = weight * vxc wv[:,0] *= .5 for kn in range(nkpts): vmat[0,kn] += krks_grad._gga_grad_sum_(ao_ks[kn], wv[0]) vmat[1,kn] += krks_grad._gga_grad_sum_(ao_ks[kn], wv[1]) + del vxc + + exc_sum += cp.sum(weight * ((rho[0,0] + rho[1,0]) * exc)) + + del rho, exc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, coords.T) + del d2ao + + wv[:,0] *= 2 + dm_nu = contract('ukpq,kgq->ukpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,ukpg->uxydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dmu_stress, wv).real + del dmu_stress, dm_nu + dm_dmu = contract('ukpq,kdgp->ukdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,ukdqg->uxydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dnu_stress, wv[:,1:4]).real + del dnu_stress, dm_dmu + + del ao_ks_strain, wv + elif xctype == 'MGGA': - ao_deriv = 2 - for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts, - sort_grids=True): + for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts, sort_grids=True): rho_a = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts[0], xctype=xctype, hermi=hermi) rho_b = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts[1], xctype=xctype, hermi=hermi) rho = cp.stack([rho_a, rho_b], axis=0) - vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=1)[1] + del rho_a, rho_b + exc, vxc = ni.eval_xc_eff(xc_code, rho, deriv=1, xctype=xctype, spin=1)[:2] wv = weight * vxc wv[:,0] *= .5 wv[:,4] *= .5 # for the factor 1/2 in tau @@ -83,15 +136,47 @@ def get_vxc(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): vmat[0,kn] += krks_grad._tau_grad_dot_(ao_ks[kn], wv[0,4]) vmat[1,kn] += krks_grad._tau_grad_dot_(ao_ks[kn], wv[1,4]) + del vxc + + exc_sum += cp.sum(weight * ((rho[0,0] + rho[1,0]) * exc)) + + del rho, exc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, coords.T) + del d2ao + + wv[:,0] *= 2 + dm_nu = contract('ukpq,kgq->ukpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,ukpg->uxydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dmu_stress, wv[:,0:4]).real + del dmu_stress, dm_nu + dm_dmu = contract('ukpq,kdgp->ukdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,ukdqg->uxydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dnu_stress, wv[:,1:4]).real + del dnu_stress, dm_dmu + dm_dnu = contract('ukpq,kdgq->ukdpg', dm_kpts, ao_ks[:,1:4].conj()) + tau_stress = contract('kxydpg,ukdpg->uxyg', ao_ks_strain[:,:,:,1:4], dm_dnu) + de_stress_rho += 2 * contract('uxyg,ug->xy', tau_stress, wv[:,4]).real + del tau_stress, dm_dnu + + del ao_ks_strain, wv + elif xctype == 'HF': pass elif xctype == 'NLC': raise NotImplementedError("NLC") else: - raise NotImplementedError(xc_code) + raise NotImplementedError(f"Unrecognized xctype = {xctype}") - exc = krhf_grad.contract_h1e_dm(cell, vmat, dm_kpts, hermi=1) - exc *= -1.0 / nkpts + de_stress_weight = exc_sum * cp.eye(3) + + exc = np.zeros((cell.natm + 3, 3)) + exc[:-3] = -krhf_grad.contract_h1e_dm(cell, vmat, dm_kpts, hermi=1) + exc[:-3] *= 1.0 / nkpts + exc[-3:] = (de_stress_rho / nkpts + de_stress_weight).get() return exc def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): @@ -115,8 +200,9 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): ao_deriv = 1 else: raise NotImplementedError(f"Unrecognized xctype = {xctype}") + eval_gto_opt = _GTOvalOpt(cell, kpts, deriv=ao_deriv) - de_grid_response_weight = cp.zeros((natm, 3), dtype=cp.float64) + de_grid_response_weight = cp.zeros((natm + 3, 3), dtype=cp.float64) g1 = 0 for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv, kpts): g0, g1 = g1, g1 + weight.size @@ -131,13 +217,13 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): else: assert rho.ndim == 1 dweight_dA = get_becke_weight_derivative(grids, natm, (g0,g1)) - dweight_dA = dweight_dA[:-3] de_grid_response_weight += cp.einsum("Adg->Ad", dweight_dA * (rho * exc)) del dweight_dA, rho, exc assert g1 == ngrids dvmat_orbital_response = cp.zeros((2,nkpts,3,nao,nao), dtype=dm_kpts.dtype) de_grid_response_rho = cp.zeros((natm, 3), dtype=dm_kpts.dtype) + de_stress_rho = cp.zeros((3,3), dtype=cp.float64) g1 = 0 for ao_ks, weight, coords in ni.block_loop(cell, grids, ao_deriv + 1, kpts): @@ -163,7 +249,19 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): dvmat_orbital_response[1,kn] += vtmp_b de_grid_response_rho[i_atom] += cp.einsum('xij,ji->x', vtmp_b, dm_kpts[1,kn]) * 2 del vtmp_a, vtmp_b - del wv, rho, aowa, aowb, vxc + del rho, aowa, aowb, vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + ao_ks_strain = ao_ks_strain[:,:,:,0] + associated_supatm_coords = grids.supatm_coords[grids.supatm_idx[g0:g1]] + ao_ks_strain += contract('kxgp,yg->kxypg', ao_ks[:,1:4], associated_supatm_coords.T) + del associated_supatm_coords + + dm_nu = contract('ukpq,kgq->ukpg', dm_kpts, ao_ks[:,0].conj()) + drho_stress = contract('kxypg,ukpg->uxyg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('uxyg,ug->xy', drho_stress, wv).real + + del ao_ks_strain, dm_nu, drho_stress, wv elif xctype == 'GGA': rho_a = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts[0], xctype=xctype, hermi=hermi) @@ -181,7 +279,26 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): dvmat_orbital_response[1,kn] += vtmp_b de_grid_response_rho[i_atom] += cp.einsum('xij,ji->x', vtmp_b, dm_kpts[1,kn]) * 2 del vtmp_a, vtmp_b - del wv, rho, vxc + del rho, vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + associated_supatm_coords = grids.supatm_coords[grids.supatm_idx[g0:g1]] + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], associated_supatm_coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, associated_supatm_coords.T) + del d2ao, associated_supatm_coords + + wv[:,0] *= 2 + dm_nu = contract('ukpq,kgq->ukpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,ukpg->uxydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dmu_stress, wv).real + del dmu_stress, dm_nu + dm_dmu = contract('ukpq,kdgp->ukdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,ukdqg->uxydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dnu_stress, wv[:,1:4]).real + del dnu_stress, dm_dmu + + del ao_ks_strain, wv elif xctype == 'MGGA': rho_a = ni.eval_rho(cell, ao_ks[:,:4], dm_kpts[0], xctype=xctype, hermi=hermi) @@ -200,14 +317,39 @@ def get_vxc_full_response(ni, cell, grids, xc_code, dm_kpts, kpts, hermi=1): dvmat_orbital_response[1,kn] += vtmp_b de_grid_response_rho[i_atom] += cp.einsum('xij,ji->x', vtmp_b, dm_kpts[1,kn]) * 2 del vtmp_a, vtmp_b - del wv, rho, vxc + del rho, vxc + + ao_ks_strain = _eval_ao_strain_derivatives(cell, coords, kpts, deriv=ao_deriv, opt=eval_gto_opt) + associated_supatm_coords = grids.supatm_coords[grids.supatm_idx[g0:g1]] + ao_ks_strain[:,:,:,0] += contract('kxgp,yg->kxypg', ao_ks[:,1:4], associated_supatm_coords.T) + d2ao = get_d2mu_dr2(ao_ks) + ao_ks_strain[:,:,:,1:4] += contract('kxdgp,yg->kxydpg', d2ao, associated_supatm_coords.T) + del d2ao, associated_supatm_coords + + wv[:,0] *= 2 + dm_nu = contract('ukpq,kgq->ukpg', dm_kpts, ao_ks[:,0].conj()) + dmu_stress = contract('kxydpg,ukpg->uxydg', ao_ks_strain, dm_nu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dmu_stress, wv[:,0:4]).real + del dmu_stress, dm_nu + dm_dmu = contract('ukpq,kdgp->ukdqg', dm_kpts, ao_ks[:,1:4]) + dnu_stress = contract('kxyqg,ukdqg->uxydg', ao_ks_strain[:,:,:,0].conj(), dm_dmu) + de_stress_rho += 2 * contract('uxydg,udg->xy', dnu_stress, wv[:,1:4]).real + del dnu_stress, dm_dmu + dm_dnu = contract('ukpq,kdgq->ukdpg', dm_kpts, ao_ks[:,1:4].conj()) + tau_stress = contract('kxydpg,ukdpg->uxyg', ao_ks_strain[:,:,:,1:4], dm_dnu) + de_stress_rho += 2 * contract('uxyg,ug->xy', tau_stress, wv[:,4]).real + del tau_stress, dm_dnu + + del ao_ks_strain, wv else: raise NotImplementedError(f"Unrecognized xctype = {xctype}") assert g1 == ngrids - exc = de_grid_response_rho.get().real - exc -= krhf_grad.contract_h1e_dm(cell, dvmat_orbital_response, dm_kpts, hermi=1) + exc = np.zeros((cell.natm + 3, 3), dtype=np.float64) + exc[:-3] = de_grid_response_rho.get().real + exc[:-3] -= krhf_grad.contract_h1e_dm(cell, dvmat_orbital_response, dm_kpts, hermi=1) + exc[-3:] = de_stress_rho.get() exc *= 1.0 / nkpts exc += de_grid_response_weight.get() return exc @@ -253,12 +395,7 @@ def energy_ee(self, dm, kpts): else: fn = get_vxc cell = self.cell - de[:-3] = fn(ni, cell, grids, xc, dm, kpts) - if isinstance(grids, BeckeGrids): - de[-3:] = np.nan - else: - de[-3:] = multigrid_v3.MultiGridNumInt(cell).energy_strain_gradient( - xc, dm, kpts=kpts, spin=1, with_j=False, with_nuc=False) + de = fn(ni, cell, grids, xc, dm, kpts) t0 = log.timer_debug1('vxc', *t0) if j_factor != 0 or k_sr != 0 or k_lr != 0: diff --git a/gpu4pyscf/pbc/grad/uks.py b/gpu4pyscf/pbc/grad/uks.py index 296271b7d..b3d311e03 100644 --- a/gpu4pyscf/pbc/grad/uks.py +++ b/gpu4pyscf/pbc/grad/uks.py @@ -60,13 +60,7 @@ def energy_ee(self, dm): fn = get_vxc cell = self.cell # UKS densities use (spin, k-point, AO, AO) ordering. - de[:-3] = fn(ni, cell, grids, xc, dm[:,None], np.zeros((1, 3))) - if isinstance(grids, BeckeGrids): - # Atomic grid response is supported; Becke-grid stress is not. - de[-3:] = np.nan - else: - de[-3:] = multigrid_v3.MultiGridNumInt(cell).energy_strain_gradient( - xc, dm, spin=1, with_j=False, with_nuc=False) + de = fn(ni, cell, grids, xc, dm[:,None], np.zeros((1, 3))) if j_factor != 0 or k_sr != 0 or k_lr != 0: de += rhf._get_ejk_derivatives(mf, dm, None, j_factor, omega, k_lr, k_sr) From ba13799fc1f33a420ee45401def51040620d9297 Mon Sep 17 00:00:00 2001 From: "henry.wang1" Date: Fri, 25 Sep 2026 08:48:18 +0800 Subject: [PATCH 5/6] Fix unstable tests --- gpu4pyscf/pbc/dft/tests/test_pbc_grids.py | 123 +++++++++++++++++++++- 1 file changed, 121 insertions(+), 2 deletions(-) diff --git a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py index 7c2fe6650..7850d055f 100644 --- a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py +++ b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py @@ -266,8 +266,126 @@ def test_xc_gradient_lda_with_response(self): mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (50,194) mf.conv_tol = 1e-10 + mf.with_df.linear_dep_threshold = 1e-10 # The default in pyscf==2.8.0 is 1e-9, in pyscf==2.14.0 is 1e-10 + + mocc = cp.array([ + [-0.85521099, -0.00000003, 0.00000002, 0.00000008, 0. , 0. , -0.3492905 , 0.34152485, 0.09633393, 0.21895687, 0.00000215, 0. , -0. , 0. , -0. , 0. , -0. , 0. , -0. , 0. , -0. , -0. , -0. , -0. ], + [-0.08831735, -0. , 0. , 0.00000002, 0. , 0. , -0.01693935, 0.02506856, -0.16650864, -0.44926685, -0.0000044 , -0. , 0. , -0.00000001, 0. , 0. , 0. , -0. , 0. , -0. , 0. , -0. , -0. , -0. ], + [ 0.2305187 , 0.00000001, -0.00000001, -0.00000017, -0. , -0. , 0.00496981, -0.11457228, -0.04837164, -0.02759234, -0.00000028, -0. , 0. , -0. , 0. , -0. , -0. , 0. , -0. , 0. , -0. , 0.00000002, 0. , -0. ], + [ 0. , 0. , 0. , -0.00045206, 0.00055216, -0.00032234, 0. , 0. , -0. , -0.00000098, 0.09957528, 0.14032745, -0.01239044, -0.00001928, -0.00002749, -0.1885912 , -0.2508752 , -0.09070458, 0.00000683, -0.0000023 , 0.00000003, -0.16268662, -0.13183595, 0.18857193], + [ 0. , -0. , 0. , -0.00045213, 0.00000304, 0.0006393 , 0. , 0. , -0. , -0.00000098, 0.09961574, -0.08086849, -0.11531499, 0.00003348, -0.00000296, -0.18855876, 0.04678102, 0.26265845, -0.00000541, -0.00000477, -0.00000003, -0.16300158, 0.2290356 , 0.01949894], + [ 0. , 0. , -0. , -0.00045206, -0.0005552 , -0.00031707, 0. , 0. , -0. , -0.00000098, 0.09960935, -0.05940578, 0.12770859, -0.00001419, 0.00003048, -0.18870952, 0.20397426, -0.17180089, -0.00000148, 0.00000708, 0. , -0.16239021, -0.09782128, -0.20848848], + [-0.00000004, -0.00000001, -0.00000001, -0.01733036, 0.02116779, -0.0123573 , 0.00000002, -0.00000001, 0. , -0.0000005 , 0.05499091, 0.07749642, -0.00684267, -0.00001065, -0.00001518, -0.09871982, -0.13132295, -0.04748016, 0.00000354, -0.00000119, 0.00000002, -0.15002503, -0.12157574, 0.17389622], + [-0.00000004, 0.00000002, -0. , -0.01733309, 0.0001167 , 0.02450857, 0.00000002, -0.00000001, 0. , -0.0000005 , 0.05501325, -0.04465996, -0.06368318, 0.00001849, -0.00000164, -0.09870284, 0.02448796, 0.137491 , -0.0000028 , -0.00000247, -0.00000001, -0.15031548, 0.21121078, 0.01798142], + [-0.00000004, -0.00000001, 0.00000002, -0.01733034, -0.02128453, -0.01215514, 0.00000002, -0.00000001, 0. , -0.00000049, 0.05500972, -0.03280709, 0.0705276 , -0.00000783, 0.00001683, -0.09878175, 0.10677222, -0.08993077, -0.00000076, 0.00000366, 0. , -0.14975169, -0.09020828, -0.19226274], + [-0. , -0. , 0. , 0.00057538, 0.00070666, 0.00040356, -0. , -0. , 0. , -0.00000006, 0.00584017, -0.00348301, 0.00748765, -0.00000083, 0.00000179, -0.00382084, 0.00412991, -0.00347849, -0.00000003, 0.00000014, 0. , -0.01384961, -0.00834279, -0.01778115], + [-0. , -0. , -0. , 0.00057538, -0.00070278, 0.00041027, -0. , -0. , 0. , -0.00000006, 0.00583817, 0.0082275 , -0.00072646, -0.00000113, -0.00000161, -0.00381845, -0.00507953, -0.00183652, 0.00000014, -0.00000005, 0. , -0.01387489, -0.01124376, 0.01608255], + [ 0. , -0.00002756, 0.00008737, 0. , -0. , -0. , 0. , -0. , 0. , 0. , -0. , 0.00000111, -0.00000238, -0.00464484, 0.00997621, -0. , -0. , 0. , 0. , -0. , 0. , -0. , -0. , -0. ], + [-0. , 0. , -0. , 0.00057547, -0.00000387, -0.0008137 , -0. , -0. , 0. , -0.00000006, 0.00584055, -0.00474138, -0.006761 , 0.00000196, -0.00000017, -0.00381779, 0.00094719, 0.0053181 , -0.00000011, -0.0000001 , -0. , -0.01390175, 0.01953353, 0.00166299], + [ 0. , -0.00008737, -0.00002756, -0. , 0. , -0. , 0. , -0. , -0. , 0. , 0. , -0.00000238, -0.00000111, -0.00997621, -0.00464484, -0. , 0. , 0. , -0. , -0. , -0. , 0. , -0. , 0. ], + [-0.00036546, 0. , 0.00000001, -0.85254518, -0.00000079, -0.00006349, -0.34537452, -0.35201109, 0.09633286, 0.00000069, -0.21895957, -0.00003897, -0.00000232, -0.00000001, -0.00000002, -0.00000579, -0. , 0. , -0. , 0. , 0. , 0.0000004 , 0. , 0. ], + [ 0.00001055, 0. , 0. , -0.08815982, -0.00000008, -0.00000657, -0.01676465, -0.02573301, -0.16652394, -0.00000477, 0.44928141, 0.00007997, 0.00000476, 0.00000001, 0.00000003, 0.00002493, 0.00000001, -0.00000001, 0. , 0. , -0. , 0.00000172, 0. , 0. ], + [-0.0003627 , -0. , -0. , 0.23011912, 0.00000021, 0.00001714, 0.00344315, 0.11543023, -0.04833963, 0.00008864, 0.02753048, 0.0000049 , 0.00000029, 0. , 0. , -0.00001898, -0. , 0.00000001, 0. , -0. , 0. , -0.00002307, -0.00000003, -0.00000002], + [ 0.00045208, -0.0004319 , -0.00047141, 0.00000001, 0. , -0. , -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000085, -0.0000332 , 0.00000294, -0.08085838, -0.11533339, 0.18863828, -0.12538688, -0.04543211, -0.21908948, 0.0732591 , 0.0000114 , 0.1626865 , 0.13183761, -0.18857432], + [ 0.00045208, 0.0006242 , -0.00013833, 0.00000001, 0. , 0. , -0.00000199, 0.00000148, 0.00000002, 0.09959752, 0.00000083, 0.00001915, 0.00002731, 0.14031084, -0.01235872, 0.18865447, 0.02345219, 0.13126238, 0.17298903, 0.15310751, 0.00004989, 0.16300146, -0.2290385 , -0.01949919], + [ 0.00045208, -0.0001923 , 0.00060974, 0.00000001, -0. , -0. , -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000082, 0.00001406, -0.00003026, -0.05945245, 0.12769211, 0.18857909, 0.10205461, -0.08598325, 0.04610051, -0.22636662, -0.0000613 , 0.16239009, 0.09782251, 0.20849111], + [ 0.01733071, -0.01655722, -0.01807177, -0.00000025, -0.00000006, 0.00000004, -0.00007611, 0.0000568 , -0.00000929, 0.05502888, -0.00002725, 0.00000581, -0.00000051, -0.04467114, -0.06371726, 0.0985865 , -0.06553283, -0.02374486, -0.1145189 , 0.03829281, 0.00000596, 0.15007133, 0.12158971, -0.17391612], + [ 0.01733071, 0.02392922, -0.00530308, -0.00000025, -0. , -0.00000007, -0.00007611, 0.0000568 , -0.00000929, 0.05502888, -0.00002726, -0.00000335, -0.00000477, 0.07751634, -0.00682772, 0.09859496, 0.01225717, 0.06860363, 0.09042202, 0.08002988, 0.00002608, 0.15036181, -0.21123494, -0.01798346], + [ 0.01733071, -0.007372 , 0.02337486, -0.00000025, 0.00000006, 0.00000004, -0.00007611, 0.0000568 , -0.00000929, 0.05502888, -0.00002726, -0.00000247, 0.00000528, -0.03284519, 0.07054497, 0.09855556, 0.05333833, -0.04493872, 0.02409691, -0.11832269, -0.00003204, 0.14979796, 0.09021866, 0.1922848 ], + [ 0.0005754 , -0.00024476, 0.00077607, 0. , -0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, 0.00000013, -0.00000071, 0.00000154, 0.00348591, -0.00748706, -0.00381672, -0.00206623, 0.00174084, -0.00093321, 0.00458235, 0.00000124, -0.01384984, -0.00834291, -0.0177814 ], + [ 0.0005754 , -0.00054972, -0.0006 , 0. , 0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, 0.00000013, 0.00000169, -0.00000015, 0.00474102, 0.00676242, -0.00381792, 0.00253862, 0.00091983, 0.00443504, -0.00148299, -0.00000023, -0.01387512, -0.01124392, 0.01608278], + [ 0. , 0. , -0. , -0. , 0.00007956, 0.00004543, -0. , -0. , -0. , -0. , -0.00000072, 0.00464234, -0.0099769 , 0.00000124, -0.00000265, -0. , 0.00000025, -0.00000021, -0.00000005, 0.00000027, 0. , -0. , 0.00000001, 0.00000003], + [ 0.0005754 , 0.00079447, -0.00017607, 0. , 0. , 0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, 0.00000013, -0.00000097, -0.00000139, -0.00822694, 0.00072464, -0.00381825, -0.00047482, -0.00265758, -0.00350182, -0.00309936, -0.00000101, -0.01390199, 0.01953381, 0.00166301], + [-0. , 0. , 0. , -0.00000001, -0.00004543, 0.00007956, -0. , -0. , 0. , 0. , 0.00000182, -0.00997689, -0.00464234, 0.00000265, 0.00000123, -0. , -0.00000021, -0.00000025, 0.00000027, 0.00000005, 0. , -0. , 0.00000003, -0.00000001], + [ 0.28149775, 0.24176118, -0.76656657, -0.00000019, 0.000001 , 0.00000057, -0.35428689, 0.34524753, 0.09633267, -0.0729876 , -0.00000073, 0.00002076, -0.00004468, -0.08713396, 0.18714652, -0. , -0. , 0. , 0. , -0. , -0. , 0. , -0. , -0. ], + [ 0.02922689, 0.025 , -0.07926899, -0. , 0.0000001 , 0.00000006, -0.01745596, 0.02545375, -0.16653584, 0.14976831, 0.0000015 , -0.0000426 , 0.00009168, 0.1787895 , -0.3840045 , -0. , -0. , 0. , -0. , 0. , 0. , 0. , -0. , -0. ], + [-0.07630111, -0.06525625, 0.20691189, -0.0000001 , -0.00000027, -0.00000015, 0.00631767, -0.11557807, -0.04821647, 0.00907678, 0.00000009, -0.00000261, 0.00000561, 0.01095183, -0.02352237, 0. , 0. , -0. , 0. , -0. , -0. , -0.00000001, 0.00000001, 0.00000001], + [ 0. , -0. , -0. , 0.00045206, -0.00055216, 0.00032234, -0. , -0. , 0. , 0.00000098, -0.09957426, -0.14032601, 0.01239032, 0.00001928, 0.00002749, -0.00000216, -0.00000041, -0.0000078 , -0.08460097, -0.2530379 , 0.1885492 , -0.16268774, -0.13183686, 0.18857324], + [ 0. , 0. , -0. , 0.00045213, -0.00000304, -0.0006393 , -0. , -0. , 0. , 0.00000098, -0.09961472, 0.08086766, 0.1153138 , -0.00003348, 0.00000296, -0.00000219, 0.00000764, 0.00000173, 0.17680287, -0.19969978, -0.1886718 , -0.1630027 , 0.22903719, 0.01949908], + [-0. , 0. , -0. , -0.00045206, -0.0005552 , -0.00031706, 0. , 0. , 0. , -0.00000097, 0.09960744, -0.05940465, 0.12770615, -0.00001419, 0.00003048, 0.18871586, -0.20398112, 0.17180667, 0.00000148, -0.00000708, -0. , -0.16239242, -0.09782261, -0.20849132], + [ 0. , -0.00000002, -0.00000003, 0.01733016, -0.0211675 , 0.01235718, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05500797, -0.07752052, 0.0068448 , 0.00001065, 0.00001519, -0.00001467, -0.00001823, -0.00001055, -0.04426999, -0.13240966, 0.09866402, -0.15004811, -0.12159439, 0.17392231], + [ 0. , 0.00000003, -0.00000001, 0.01733289, -0.00011665, -0.02450826, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05503032, 0.04467385, 0.06370299, -0.0000185 , 0.00000164, -0.00001468, 0.00000731, 0.00001976, 0.0925174 , -0.10449889, -0.09872817, -0.1503386 , 0.21124272, 0.01798381], + [ 0.00000004, 0.00000003, -0.00000008, -0.01732993, -0.02128398, -0.01215483, 0.00000002, 0.00000004, -0. , -0.00000058, 0.05504905, -0.0328306 , 0.07057814, -0.00000784, 0.00001684, 0.09874795, -0.10673569, 0.08990001, 0.00000076, -0.00000366, -0. , -0.14978715, -0.09022915, -0.19230721], + [ 0. , -0. , 0. , 0.00057538, 0.00070666, 0.00040356, -0. , -0. , -0. , -0.00000006, 0.0058406 , -0.00348326, 0.0074882 , -0.00000083, 0.00000179, 0.003821 , -0.00413009, 0.00347864, 0.00000003, -0.00000014, -0. , -0.01384994, -0.00834298, -0.01778157], + [-0. , 0. , 0. , -0.00057538, 0.00070278, -0.00041027, 0. , 0. , 0. , 0.00000006, -0.00583821, -0.00822755, 0.00072647, 0.00000113, 0.00000161, -0.00000048, -0.00000058, -0.00000037, -0.00171288, -0.00512316, 0.00381749, -0.01387487, -0.01124375, 0.01608252], + [ 0.00008638, -0.00000918, 0.00002912, -0. , -0. , -0. , -0.00000039, 0.00000028, 0.00000066, 0.01037336, 0.0000001 , 0.00000037, -0.00000079, -0.00154829, 0.00332542, 0. , -0. , 0. , 0. , -0. , -0. , -0. , -0. , -0. ], + [-0. , -0. , 0. , -0.00057547, 0.00000387, 0.0008137 , 0. , 0. , 0. , 0.00000006, -0.00584058, 0.00474141, 0.00676104, -0.00000196, 0.00000017, -0.00000048, 0.00000026, 0.00000064, 0.00357966, -0.00404325, -0.00381997, -0.01390173, 0.0195335 , 0.00166298], + [-0. , 0.00008737, 0.00002756, 0. , -0. , 0. , -0. , 0. , 0. , -0. , -0. , 0.00000238, 0.00000111, 0.00997538, 0.00464446, 0. , -0. , -0. , 0. , 0. , 0. , -0. , 0. , -0. ], + [-0.00036519, 0.00000035, -0.0000011 , 0.28421093, -0.69800517, -0.39855274, -0.34537609, -0.3520123 , 0.09633286, -0.00000217, 0.072973 , 0.08710356, -0.18716625, 0.0000208 , -0.00004466, 0.00000193, 0.00000417, -0.00000352, -0. , 0. , 0. , -0.00000013, 0.00000016, 0.00000034], + [ 0.00001058, 0.00000004, -0.00000011, 0.02938968, -0.07217921, -0.04121348, -0.01676482, -0.02573314, -0.16652394, 0.0000011 , -0.14973273, -0.1787271 , 0.38404494, -0.00004268, 0.00009163, -0.0000083 , -0.00001798, 0.00001515, 0. , -0. , 0. , -0.00000058, 0.00000069, 0.00000147], + [-0.00036278, -0.00000009, 0.0000003 , -0.07671439, 0.18840574, 0.10757746, 0.00344358, 0.11543056, -0.04833963, 0.00008899, -0.00917512, -0.01095181, 0.02353301, -0.00000262, 0.00000562, 0.00000632, 0.00001368, -0.00001153, -0. , 0. , -0. , 0.00000772, -0.00000922, -0.00001969], + [-0.00045208, 0.0004319 , 0.00047141, 0. , -0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000092, 0.00003358, -0.00000271, 0.08085838, 0.11533339, -0.00007569, 0.07859684, -0.21723151, 0.04229099, 0.12644661, 0.18866153, 0.16268795, 0.13183849, -0.18857247], + [-0.00045208, -0.0006242 , 0.00013833, 0. , -0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959752, -0.0000009 , -0.00001957, -0.00002745, -0.14031084, 0.01235872, -0.00005946, 0.22742656, -0.04054806, -0.0883944 , 0.09993416, -0.18860023, 0.16300291, -0.22903763, -0.01949734], + [ 0.00045208, -0.0001923 , 0.00060974, -0.00000001, 0.00000001, 0.00000001, -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000125, 0.00001419, -0.00003053, -0.05945245, 0.12769211, -0.1885797 , -0.10205394, 0.08598269, -0.04610051, 0.22636662, 0.0000613 , 0.16239299, 0.09782165, 0.20848927], + [-0.01733075, 0.01655727, 0.01807176, -0.00000001, 0.00000016, 0.00000029, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000167, 0.00002485, 0.00003476, 0.04467113, 0.06371725, -0.00004099, 0.04107512, -0.11353227, 0.02210568, 0.06609417, 0.09861411, 0.15005196, 0.12157799, -0.17394032], + [-0.01733075, -0.02392926, 0.00530305, -0.00000001, 0.00000033, -0.00000001, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000169, -0.0000426 , 0.00000337, -0.07751633, 0.00682773, -0.00003251, 0.11886017, -0.02118961, -0.04620412, 0.052236 , -0.09858207, 0.15034244, -0.21124606, -0.01800794], + [ 0.01733079, -0.00737202, 0.02337491, 0.0000002 , -0.00000017, -0.0000001 , -0.00007611, 0.00005685, -0.00000929, 0.05502879, 0.00003277, 0.00001547, -0.00003327, -0.03284518, 0.07054496, -0.09855875, -0.05333488, 0.04493581, -0.02409691, 0.11832269, 0.00003204, 0.14975912, 0.09022976, 0.19230857], + [ 0.0005754 , -0.00024476, 0.00077607, -0. , -0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, -0.00000038, -0.00000087, 0.00000187, 0.00348591, -0.00748706, 0.00381761, 0.00206527, -0.00174004, 0.00093321, -0.00458235, -0.00000124, -0.0138498 , -0.00834292, -0.01778143], + [-0.0005754 , 0.00054972, 0.0006 , -0. , -0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, -0.00000208, -0.00000013, -0.00474102, -0.00676242, 0.00000197, -0.00159035, 0.00439734, -0.0008561 , -0.00255966, -0.00381908, -0.0138751 , -0.01124391, 0.0160828 ], + [ 0. , 0. , -0. , 0.00008638, 0.00002652, 0.00001515, -0. , -0. , -0. , 0.0000001 , -0.01037499, 0.0015456 , -0.00332574, 0.0000005 , -0.00000107, -0.00000031, 0.00000008, -0.00000007, 0.00000005, -0.00000027, -0. , 0.00000003, 0. , 0.00000001], + [-0.0005754 , -0.00079447, 0.00017607, -0. , -0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000144, 0.00000151, 0.00822694, -0.00072464, 0.00000164, -0.0046036 , 0.00082015, 0.00178937, -0.00202297, 0.00381784, -0.01390196, 0.01953382, 0.00166304], + [ 0. , 0. , 0. , 0.00000001, 0.00004543, -0.00007956, 0. , 0. , -0. , -0. , -0.00000182, 0.0099769 , 0.00464234, -0.00000211, -0.00000098, 0. , 0.00000021, 0.00000025, 0.00000009, 0.00000002, 0.00000026, 0. , -0.00000003, 0.00000001], + [ 0.2814978 , -0.78474669, 0.17391196, -0.00000018, -0.00000001, -0.00000116, -0.35428689, 0.34524753, 0.09633267, -0.0729876 , -0.00000073, 0.00002828, 0.00004033, 0.20564062, -0.01811304, -0. , -0. , -0. , 0. , 0. , 0. , 0. , 0. , 0. ], + [ 0.0292269 , -0.08114896, 0.01798386, -0. , -0. , -0.00000012, -0.01745596, 0.02545375, -0.16653584, 0.14976832, 0.00000149, -0.00005802, -0.00008276, -0.4219524 , 0.037166 , -0. , -0. , -0. , -0. , -0. , -0. , 0. , 0. , 0. ], + [-0.07630112, 0.21181907, -0.04694237, -0.0000001 , 0. , 0.00000031, 0.00631767, -0.11557807, -0.04821647, 0.00907678, 0.00000009, -0.00000355, -0.00000507, -0.02584689, 0.00227662, 0. , 0. , 0. , 0. , 0. , -0. , -0.00000001, -0.00000002, -0. ], + [ 0. , -0. , -0. , 0.00045206, -0.00055216, 0.00032234, -0. , -0. , 0. , 0.00000098, -0.09957426, -0.14032601, 0.01239032, 0.00001928, 0.00002749, -0.00000218, -0.00000537, 0.00000572, 0.08460097, 0.2530379 , -0.1885492 , -0.16268774, -0.13183687, 0.18857324], + [-0. , -0. , 0. , -0.00045213, 0.00000304, 0.0006393 , 0. , 0. , 0. , -0.00000098, 0.09961384, -0.08086694, -0.11531278, 0.00003348, -0.00000296, 0.1885651 , -0.04678259, -0.26266728, 0.00000541, 0.00000477, 0.00000003, -0.16300379, 0.22903872, 0.0194992 ], + [ 0. , -0. , 0. , 0.00045206, 0.0005552 , 0.00031707, -0. , -0. , 0. , 0.00000098, -0.09960833, 0.05940517, -0.12770728, 0.00001419, -0.00003048, -0.00000219, 0.00000702, 0.00000355, 0.26138042, 0.05318012, 0.18864902, -0.16239133, -0.09782196, -0.20848992], + [ 0. , -0.00000002, -0.00000003, 0.01733016, -0.02116754, 0.0123571 , -0.00000002, -0.00000002, -0. , 0.00000054, -0.05500797, -0.07752052, 0.0068448 , 0.00001065, 0.00001519, -0.00001467, -0.00002079, -0.00000356, 0.04426999, 0.13240966, -0.09866402, -0.15004811, -0.12159384, 0.17392269], + [ 0.00000004, -0.00000009, 0.00000002, -0.01733268, 0.0001167 , 0.02450793, 0.00000002, 0.00000004, -0. , -0.00000059, 0.05505259, -0.04469196, -0.06372882, 0.0000185 , -0.00000164, 0.09866906, -0.02447958, -0.13744397, 0.0000028 , 0.00000247, 0.00000001, -0.15035106, 0.21125963, 0.01798558], + [ 0. , -0.00000001, 0.00000003, 0.01733014, 0.02128429, 0.01215495, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05502679, 0.03281729, -0.07054953, 0.00000784, -0.00001684, -0.00001469, 0.00001828, -0.0000105 , 0.13677514, 0.0278281 , 0.09871625, -0.14977473, -0.09022162, -0.19229194], + [-0. , 0. , -0. , -0.00057538, -0.00070666, -0.00040356, 0. , 0. , 0. , 0.00000006, -0.0058402 , 0.00348303, -0.00748769, 0.00000083, -0.00000179, -0.00000048, 0.00000061, -0.00000032, 0.00529207, 0.00107672, 0.00381951, -0.01384959, -0.00834277, -0.01778112], + [-0. , 0. , 0. , -0.00057538, 0.00070278, -0.00041027, 0. , 0. , 0. , 0.00000006, -0.00583821, -0.00822755, 0.00072647, 0.00000113, 0.00000161, -0.00000048, -0.00000068, -0.00000009, 0.00171288, 0.00512316, -0.00381749, -0.01387487, -0.01124374, 0.01608252], + [-0.00004319, -0.00003207, -0.00007416, 0. , 0. , -0. , 0.0000002 , -0.00000014, -0.00000033, -0.00518668, -0.00000005, -0.00000211, 0.00000095, -0.00266314, -0.00933172, -0. , 0. , 0. , -0. , 0. , 0. , 0. , -0. , 0. ], + [ 0. , 0. , -0. , 0.00057547, -0.00000387, -0.0008137 , -0. , -0. , -0. , -0.00000006, 0.00584098, -0.00474173, -0.0067615 , 0.00000196, -0.00000017, 0.00381795, -0.00094723, -0.00531833, 0.00000011, 0.0000001 , 0. , -0.01390208, 0.01953399, 0.00166303], + [-0.0000748 , -0.00001591, 0.00005044, 0. , -0. , -0. , 0.00000034, -0.00000024, -0.00000057, -0.00898359, -0.00000009, 0.00000064, -0.00000137, -0.00268176, 0.00575931, -0. , -0. , 0. , 0. , -0. , -0. , 0. , -0. , -0. ], + [-0.00036519, -0.00000113, 0.00000024, 0.28412074, 0.00382786, 0.80379856, -0.34537609, -0.3520123 , 0.09633286, -0.00000217, 0.07296363, 0.1185591 , 0.16900692, -0.00004906, 0.00000436, 0.00000193, 0.00000096, 0.00000537, -0. , -0. , -0. , -0.00000013, -0.00000037, -0.00000003], + [ 0.00001058, -0.00000012, 0.00000003, 0.02938035, 0.00039583, 0.08311908, -0.01676482, -0.02573314, -0.16652394, 0.0000011 , -0.1497135 , -0.24327048, -0.34678397, 0.00010067, -0.00000895, -0.00000832, -0.00000413, -0.00002314, 0. , 0. , 0. , -0.00000057, -0.00000162, -0.00000014], + [-0.00036278, 0.0000003 , -0.00000007, -0.07669004, -0.00103322, -0.21696152, 0.00344358, 0.11543056, -0.04833963, 0.00008899, -0.00917395, -0.01490681, -0.02124978, 0.00000617, -0.00000055, 0.00000633, 0.00000314, 0.00001761, -0. , -0. , 0. , 0.00000766, 0.00002168, 0.00000186], + [-0.00045208, 0.0004319 , 0.00047141, 0. , 0. , 0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000092, 0.00003354, -0.00000324, 0.08085838, 0.11533339, 0.00007508, -0.07859766, 0.21723122, -0.04229099, -0.12644661, -0.18866153, 0.16268795, 0.13183558, -0.1885745 ], + [ 0.00045208, 0.0006242 , -0.00013833, -0.00000001, -0. , -0.00000001, -0.00000199, 0.00000148, 0.00000002, 0.09959752, 0.00000126, 0.00001932, 0.00002756, 0.14031084, -0.01235872, -0.18865508, -0.02345204, -0.13126153, -0.17298903, -0.15310751, -0.00004989, 0.16300436, -0.22903647, -0.01949901], + [-0.00045208, 0.0001923 , -0.00060974, 0. , -0. , 0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000089, -0.00001445, 0.00003047, 0.05945245, -0.12769211, 0.00001593, 0.14882962, 0.17668265, -0.1307088 , -0.02667045, 0.18861157, 0.16239154, 0.09782048, 0.20849094], + [-0.01733075, 0.0165572 , 0.01807182, -0.00000001, -0.00000017, -0.00000028, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000167, 0.00001837, -0.00003859, 0.04467115, 0.06371724, 0.0000378 , -0.04107936, 0.11353074, -0.02210568, -0.06609417, -0.09861411, 0.15005189, 0.12161643, -0.17391351], + [ 0.01733079, 0.02392927, -0.0053031 , 0.0000002 , 0. , 0.0000002 , -0.00007611, 0.00005685, -0.00000929, 0.05502879, 0.00003277, 0.00002107, 0.00003004, 0.07751632, -0.00682772, -0.09859815, -0.01225638, -0.06859919, -0.09042202, -0.08002988, -0.00002608, 0.15032303, -0.21126119, -0.01798573], + [-0.01733075, 0.00737197, -0.0233749 , -0.00000001, 0.00000017, -0.00000028, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000169, -0.00004253, 0.00000416, 0.03284519, -0.07054496, 0.00000689, 0.07778429, 0.09233825, -0.06832204, -0.01394076, 0.098588 , 0.14977852, 0.09024543, 0.19228674], + [-0.0005754 , 0.00024476, -0.00077607, -0. , -0. , 0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000116, -0.00000174, -0.00348591, 0.00748706, 0.00000011, -0.00301303, -0.00357596, 0.00264594, 0.00053989, -0.00381807, -0.01384982, -0.00834294, -0.0177814 ], + [-0.0005754 , 0.00054972, 0.0006 , -0. , 0. , 0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, -0.00000203, 0.0000005 , -0.00474102, -0.00676242, -0.00000108, 0.00159153, -0.00439692, 0.0008561 , 0.00255966, 0.00381908, -0.0138751 , -0.01124395, 0.01608277], + [-0. , 0. , -0. , -0.00004319, -0.00007927, 0.00001564, 0. , 0. , 0. , -0.00000005, 0.00518891, -0.00885413, 0.00397106, 0.00000076, 0.00000195, 0.00000016, -0.00000029, -0. , -0.00000015, -0.00000003, 0.00000022, -0.00000001, 0.00000001, -0.00000002], + [ 0.0005754 , 0.00079447, -0.00017607, -0. , 0. , 0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, -0.00000038, -0.00000118, -0.00000169, -0.00822694, 0.00072464, 0.00381913, 0.0004746 , 0.00265635, 0.00350182, 0.00309936, 0.00000101, -0.01390194, 0.01953384, 0.00166301], + [-0. , 0. , 0. , -0.00007481, 0.00004593, 0.00002622, 0. , 0. , 0. , -0.00000009, 0.00898438, 0.00268186, -0.00576007, 0.00000091, -0.00000125, 0.00000027, 0.00000015, -0.00000012, -0.00000015, -0.00000019, -0.00000013, -0.00000002, 0.00000001, 0.00000001], + [ 0.28149777, 0.54298554, 0.59265459, -0.00000018, -0.000001 , 0.00000059, -0.35428689, 0.34524753, 0.09633267, -0.0729876 , -0.00000069, -0.00004904, 0.00000435, -0.11850666, -0.16903348, -0. , 0. , 0. , -0. , 0. , -0. , 0. , -0. , 0. ], + [ 0.02922689, 0.05614896, 0.06128513, -0. , -0.0000001 , 0.00000006, -0.01745596, 0.02545375, -0.16653584, 0.14976831, 0.00000142, 0.00010062, -0.00000892, 0.24316291, 0.3468385 , -0. , 0. , 0. , 0. , -0. , 0. , 0. , -0. , 0. ], + [-0.07630111, -0.14656283, -0.15996951, -0.0000001 , 0.00000027, -0.00000016, 0.00631767, -0.11557807, -0.04821647, 0.00907678, 0.00000009, 0.00000616, -0.00000055, 0.01489506, 0.02124575, 0. , -0. , -0. , -0. , 0. , -0. , -0.00000001, 0.00000001, -0.00000001], + [-0. , 0. , 0. , -0.00045206, 0.00055216, -0.00032234, 0. , 0. , 0. , -0.00000098, 0.09957338, 0.14032477, -0.01239021, -0.00001928, -0.00002749, 0.18859754, 0.25088363, 0.09070763, -0.00000683, 0.0000023 , -0.00000003, -0.16268883, -0.13183775, 0.18857451], + [ 0. , 0. , -0. , 0.00045213, -0.00000304, -0.0006393 , -0. , -0. , 0. , 0.00000098, -0.09961472, 0.08086766, 0.1153138 , -0.00003348, 0.00000296, -0.00000215, -0.00000657, 0.00000431, -0.17680287, 0.19969978, 0.1886718 , -0.1630027 , 0.22903719, 0.01949907], + [ 0. , -0. , 0. , 0.00045206, 0.0005552 , 0.00031706, -0. , -0. , 0. , 0.00000098, -0.09960833, 0.05940517, -0.12770728, 0.00001419, -0.00003048, -0.00000215, -0.00000233, -0.0000075 , -0.26138042, -0.05318012, -0.18864902, -0.16239133, -0.09782195, -0.20848992], + [ 0.00000004, 0.00000006, 0.00000006, -0.01732996, 0.02116724, -0.01235698, 0.00000002, 0.00000004, -0. , -0.00000059, 0.05503023, 0.07755195, -0.00684758, -0.00001066, -0.00001519, 0.09868604, 0.13127803, 0.04746391, -0.00000353, 0.00000119, -0.00000002, -0.15006055, -0.12160387, 0.17393644], + [ 0. , 0.00000003, -0.00000001, 0.01733289, -0.00011674, -0.02450826, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05503032, 0.04467385, 0.06370299, -0.0000185 , 0.00000164, -0.00001465, -0.00000004, 0.00002109, -0.0925174 , 0.10449889, 0.09872817, -0.1503386 , 0.21124267, 0.01798447], + [ 0. , -0.00000001, 0.00000003, 0.01733014, 0.02128424, 0.01215503, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05502679, 0.03281729, -0.07054953, 0.00000784, -0.00001684, -0.00001467, 0.00001344, -0.00001622, -0.13677514, -0.02782809, -0.09871625, -0.14977473, -0.09022222, -0.19229166], + [-0. , 0. , -0. , -0.00057538, -0.00070666, -0.00040356, 0. , 0. , 0. , 0.00000006, -0.0058402 , 0.00348303, -0.00748769, 0.00000083, -0.00000179, -0.00000048, 0.00000042, -0.00000055, -0.00529207, -0.00107672, -0.00381951, -0.01384959, -0.00834278, -0.01778112], + [ 0. , -0. , -0. , 0.00057538, -0.00070278, 0.00041027, -0. , -0. , -0. , -0.00000006, 0.00583861, 0.00822811, -0.00072651, -0.00000113, -0.00000161, 0.00381861, 0.00507974, 0.00183659, -0.00000014, 0.00000005, -0. , -0.01387522, -0.01124403, 0.01608293], + [-0.00004319, 0.00006881, -0.00004234, 0. , 0. , 0. , 0.0000002 , -0.00000014, -0.00000033, -0.00518668, -0.00000005, 0.00000064, 0.00000223, 0.00885569, -0.00396864, -0. , -0. , -0. , 0. , 0. , -0. , 0. , 0. , 0. ], + [-0. , -0. , 0. , -0.00057547, 0.00000387, 0.0008137 , 0. , 0. , 0. , 0.00000006, -0.00584058, 0.00474141, 0.00676104, -0.00000196, 0.00000017, -0.00000048, -0.00000003, 0.00000069, -0.00357966, 0.00404325, 0.00381997, -0.01390173, 0.0195335 , 0.00166299], + [ 0.0000748 , 0.00001591, -0.00005044, -0. , 0. , 0. , -0.00000034, 0.00000024, 0.00000057, 0.00898359, 0.00000009, -0.00000064, 0.00000137, 0.00268132, -0.00575952, 0. , 0. , -0. , 0. , 0. , -0. , -0. , 0. , 0. ], + [-0.00036519, 0.00000078, 0.00000085, 0.28421012, 0.69417811, -0.40518233, -0.34537609, -0.3520123 , 0.09633286, -0.00000217, 0.07302293, -0.20562368, 0.01816164, 0.00002827, 0.00004031, 0.00000193, -0.00000513, -0.00000186, 0. , -0. , 0. , -0.00000013, 0.00000021, -0.00000031], + [ 0.00001058, 0.00000008, 0.00000009, 0.02938959, 0.07178346, -0.04189903, -0.01676482, -0.02573314, -0.16652394, 0.0000011 , -0.14983518, 0.4219176 , -0.03726573, -0.000058 , -0.00008272, -0.00000831, 0.0000221 , 0.000008 , -0. , 0. , -0. , -0.00000057, 0.00000093, -0.00000133], + [-0.00036278, -0.00000021, -0.00000023, -0.07671417, -0.18737274, 0.10936692, 0.00344358, 0.11543056, -0.04833963, 0.00008899, -0.0091814 , 0.02585372, -0.00228352, -0.00000356, -0.00000507, 0.00000633, -0.00001682, -0.00000609, 0. , -0. , 0. , 0.00000769, -0.00001243, 0.00001785], + [ 0.00045208, -0.0004319 , -0.00047141, -0.00000001, -0.00000001, 0.00000001, -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000128, -0.00003351, 0.00000297, -0.08085838, -0.11533339, -0.18863889, 0.12538606, 0.04543181, 0.21908948, -0.0732591 , -0.0000114 , 0.1626894 , 0.13183645, -0.18857265], + [-0.00045208, -0.0006242 , 0.00013833, 0. , 0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959752, -0.0000009 , -0.00001914, -0.00002775, -0.14031084, 0.01235872, 0.00005884, -0.22742641, 0.04054892, 0.0883944 , -0.09993416, 0.18860023, 0.16300291, -0.22903733, -0.01950086], + [-0.00045208, 0.0001923 , -0.00060974, 0. , 0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000089, -0.00001397, 0.00003069, 0.05945245, -0.12769211, -0.00001655, -0.14882896, -0.1766832 , 0.1307088 , 0.02667045, -0.18861157, 0.16239154, 0.09782369, 0.20848944], + [ 0.01733079, -0.01655725, -0.01807181, 0.0000002 , 0.00000017, -0.0000001 , -0.00007611, 0.00005685, -0.00000929, 0.05502879, 0.00003279, -0.00003652, 0.00000323, -0.04467114, -0.06371725, -0.09858969, 0.06552859, 0.02374333, 0.1145189 , -0.03829281, -0.00000596, 0.15003252, 0.1216047 , -0.17393772], + [-0.01733075, -0.02392924, 0.00530313, -0.00000001, -0.00000033, -0.00000001, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000168, 0.00001768, -0.00003892, -0.07751633, 0.00682771, 0.00002932, -0.11885938, 0.02119405, 0.04620412, -0.052236 , 0.09858207, 0.1503424 , -0.21125006, -0.01796125], + [-0.01733075, 0.00737206, -0.02337487, -0.00000002, -0.00000015, 0.00000029, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000168, 0.00002423, 0.00003522, 0.03284518, -0.07054497, -0.00001008, -0.07778084, -0.09234115, 0.06832204, 0.01394076, -0.098588 , 0.14977855, 0.09020299, 0.19230662], + [-0.0005754 , 0.00024476, -0.00077607, -0. , 0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000058, -0.00000201, -0.00348591, 0.00748706, 0.00000077, 0.00301208, 0.00357676, -0.00264594, -0.00053989, 0.00381807, -0.01384982, -0.00834289, -0.01778143], + [ 0.0005754 , -0.00054972, -0.0006 , -0. , 0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, -0.00000038, 0.00000205, -0.00000018, 0.00474102, 0.00676242, 0.0038188 , -0.00253745, -0.00091941, -0.00443504, 0.00148299, 0.00000023, -0.01387508, -0.01124394, 0.0160828 ], + [-0. , 0. , -0. , -0.00004318, -0.00002681, -0.00007622, 0. , 0. , 0. , -0.00000005, 0.0051868 , 0.00266619, 0.00933158, -0.00000199, 0.00000068, 0.00000016, -0.00000005, 0.00000029, 0.00000015, 0.00000003, -0.00000022, -0.00000001, -0.00000002, -0.00000001], + [-0.0005754 , -0.00079447, 0.00017607, -0. , 0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000093, 0.00000187, 0.00822694, -0.00072464, -0.00000076, 0.00460338, -0.00082138, -0.00178937, 0.00202297, -0.00381784, -0.01390196, 0.01953383, 0.00166298], + [ 0. , 0. , 0. , 0.00007481, -0.00004593, -0.00002622, -0. , -0. , -0. , 0.00000009, -0.00898438, -0.00268186, 0.00576007, -0.00000037, 0.0000015 , -0.00000027, -0.00000015, 0.00000012, -0.00000021, 0.00000012, -0.00000013, 0.00000002, -0.00000001, -0.00000001], + ]) + dm0 = 2 * mocc @ mocc.T + dm0 = dm0[None,:,:] - mf.kernel() + mf.kernel(dm0 = dm0) dm = mf.make_rdm1() if dm.ndim == 2: @@ -291,7 +409,7 @@ def test_xc_gradient_lda_with_response(self): # It can match down to 1e-9, if the finite difference is computed using the same dm from SCF. # However if we save the finite difference result, it suffers from the numerical instability of dm, and the a 3e-7 error is observed. - assert np.max(np.abs(test_gradient - ref_gradient)) < 5e-7 + assert np.max(np.abs(test_gradient - ref_gradient)) < 5e-8 def test_xc_gradient_gga_with_response(self): cell = pyscf.M( @@ -489,6 +607,7 @@ def test_xc_gradient_unrestricted_no_k_with_response(self): mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (50,194) mf.conv_tol = 1e-10 + mf.with_df.linear_dep_threshold = 1e-10 # The default in pyscf==2.8.0 is 1e-9, in pyscf==2.14.0 is 1e-10 mf.kernel() From 39c37e52ba1eae6057c464020e3ff8a5d7e9c3be Mon Sep 17 00:00:00 2001 From: "henry.wang1" Date: Fri, 25 Sep 2026 17:16:41 +0800 Subject: [PATCH 6/6] Change the test to a much more stable one --- gpu4pyscf/pbc/dft/tests/test_pbc_grids.py | 174 ++++------------------ 1 file changed, 26 insertions(+), 148 deletions(-) diff --git a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py index 7850d055f..00c783135 100644 --- a/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py +++ b/gpu4pyscf/pbc/dft/tests/test_pbc_grids.py @@ -246,170 +246,48 @@ def test_becke_weight_derivative(self): def test_xc_gradient_lda_with_response(self): cell = pyscf.M( - a = np.eye(3) * 3.5668, - atom = ''' - C 0. 0. 0. - C 0.8917 0.8917 0.8917 - C 1.7834 1.7834 0. - C 2.6751 2.6751 0.8917 - C 1.7834 0. 1.7834 - C 2.6751 0.8917 2.6751 - C 0. 1.7834 1.7834 - C 0.8917 2.6751 2.6751 - ''', - basis = 'def2-svp', + a = np.array([ + [15.9069652593, 0, 0], + [0, 15.9069652593, 0], + [0, 6, 15.9069652593], + ]), + atom = """ + O 15.43509000 9.59549000 8.94968000 + H 15.05724000 9.21878000 9.73314000 + H 0.51550474 9.33856000 9.01857000 + He 2.51550474 9.33856000 9.01857000 + """, + basis = {'default': 'def2-svp', 'He': """ + #BASIS SET: (2s) -> [1s] + He S + 0.2432879285E+01 0.4301284983E+00 + 0.4330512863E+00 0.6789135305E+00 + He F + 1.0 1.0 + He G + 1.2 1.0 + END + """}, verbose = 0, ) kpts = np.array([[0,0,0]]) - mf = krks.KRKS(cell, xc="LDA0", kpts=kpts).density_fit(auxbasis='def2-universal-jkfit') + mf = rks.RKS(cell, xc="LDA0").density_fit(auxbasis='def2-universal-jkfit') mf.grids = gen_grid.BeckeGrids(cell) mf.grids.atom_grid = (50,194) mf.conv_tol = 1e-10 mf.with_df.linear_dep_threshold = 1e-10 # The default in pyscf==2.8.0 is 1e-9, in pyscf==2.14.0 is 1e-10 - mocc = cp.array([ - [-0.85521099, -0.00000003, 0.00000002, 0.00000008, 0. , 0. , -0.3492905 , 0.34152485, 0.09633393, 0.21895687, 0.00000215, 0. , -0. , 0. , -0. , 0. , -0. , 0. , -0. , 0. , -0. , -0. , -0. , -0. ], - [-0.08831735, -0. , 0. , 0.00000002, 0. , 0. , -0.01693935, 0.02506856, -0.16650864, -0.44926685, -0.0000044 , -0. , 0. , -0.00000001, 0. , 0. , 0. , -0. , 0. , -0. , 0. , -0. , -0. , -0. ], - [ 0.2305187 , 0.00000001, -0.00000001, -0.00000017, -0. , -0. , 0.00496981, -0.11457228, -0.04837164, -0.02759234, -0.00000028, -0. , 0. , -0. , 0. , -0. , -0. , 0. , -0. , 0. , -0. , 0.00000002, 0. , -0. ], - [ 0. , 0. , 0. , -0.00045206, 0.00055216, -0.00032234, 0. , 0. , -0. , -0.00000098, 0.09957528, 0.14032745, -0.01239044, -0.00001928, -0.00002749, -0.1885912 , -0.2508752 , -0.09070458, 0.00000683, -0.0000023 , 0.00000003, -0.16268662, -0.13183595, 0.18857193], - [ 0. , -0. , 0. , -0.00045213, 0.00000304, 0.0006393 , 0. , 0. , -0. , -0.00000098, 0.09961574, -0.08086849, -0.11531499, 0.00003348, -0.00000296, -0.18855876, 0.04678102, 0.26265845, -0.00000541, -0.00000477, -0.00000003, -0.16300158, 0.2290356 , 0.01949894], - [ 0. , 0. , -0. , -0.00045206, -0.0005552 , -0.00031707, 0. , 0. , -0. , -0.00000098, 0.09960935, -0.05940578, 0.12770859, -0.00001419, 0.00003048, -0.18870952, 0.20397426, -0.17180089, -0.00000148, 0.00000708, 0. , -0.16239021, -0.09782128, -0.20848848], - [-0.00000004, -0.00000001, -0.00000001, -0.01733036, 0.02116779, -0.0123573 , 0.00000002, -0.00000001, 0. , -0.0000005 , 0.05499091, 0.07749642, -0.00684267, -0.00001065, -0.00001518, -0.09871982, -0.13132295, -0.04748016, 0.00000354, -0.00000119, 0.00000002, -0.15002503, -0.12157574, 0.17389622], - [-0.00000004, 0.00000002, -0. , -0.01733309, 0.0001167 , 0.02450857, 0.00000002, -0.00000001, 0. , -0.0000005 , 0.05501325, -0.04465996, -0.06368318, 0.00001849, -0.00000164, -0.09870284, 0.02448796, 0.137491 , -0.0000028 , -0.00000247, -0.00000001, -0.15031548, 0.21121078, 0.01798142], - [-0.00000004, -0.00000001, 0.00000002, -0.01733034, -0.02128453, -0.01215514, 0.00000002, -0.00000001, 0. , -0.00000049, 0.05500972, -0.03280709, 0.0705276 , -0.00000783, 0.00001683, -0.09878175, 0.10677222, -0.08993077, -0.00000076, 0.00000366, 0. , -0.14975169, -0.09020828, -0.19226274], - [-0. , -0. , 0. , 0.00057538, 0.00070666, 0.00040356, -0. , -0. , 0. , -0.00000006, 0.00584017, -0.00348301, 0.00748765, -0.00000083, 0.00000179, -0.00382084, 0.00412991, -0.00347849, -0.00000003, 0.00000014, 0. , -0.01384961, -0.00834279, -0.01778115], - [-0. , -0. , -0. , 0.00057538, -0.00070278, 0.00041027, -0. , -0. , 0. , -0.00000006, 0.00583817, 0.0082275 , -0.00072646, -0.00000113, -0.00000161, -0.00381845, -0.00507953, -0.00183652, 0.00000014, -0.00000005, 0. , -0.01387489, -0.01124376, 0.01608255], - [ 0. , -0.00002756, 0.00008737, 0. , -0. , -0. , 0. , -0. , 0. , 0. , -0. , 0.00000111, -0.00000238, -0.00464484, 0.00997621, -0. , -0. , 0. , 0. , -0. , 0. , -0. , -0. , -0. ], - [-0. , 0. , -0. , 0.00057547, -0.00000387, -0.0008137 , -0. , -0. , 0. , -0.00000006, 0.00584055, -0.00474138, -0.006761 , 0.00000196, -0.00000017, -0.00381779, 0.00094719, 0.0053181 , -0.00000011, -0.0000001 , -0. , -0.01390175, 0.01953353, 0.00166299], - [ 0. , -0.00008737, -0.00002756, -0. , 0. , -0. , 0. , -0. , -0. , 0. , 0. , -0.00000238, -0.00000111, -0.00997621, -0.00464484, -0. , 0. , 0. , -0. , -0. , -0. , 0. , -0. , 0. ], - [-0.00036546, 0. , 0.00000001, -0.85254518, -0.00000079, -0.00006349, -0.34537452, -0.35201109, 0.09633286, 0.00000069, -0.21895957, -0.00003897, -0.00000232, -0.00000001, -0.00000002, -0.00000579, -0. , 0. , -0. , 0. , 0. , 0.0000004 , 0. , 0. ], - [ 0.00001055, 0. , 0. , -0.08815982, -0.00000008, -0.00000657, -0.01676465, -0.02573301, -0.16652394, -0.00000477, 0.44928141, 0.00007997, 0.00000476, 0.00000001, 0.00000003, 0.00002493, 0.00000001, -0.00000001, 0. , 0. , -0. , 0.00000172, 0. , 0. ], - [-0.0003627 , -0. , -0. , 0.23011912, 0.00000021, 0.00001714, 0.00344315, 0.11543023, -0.04833963, 0.00008864, 0.02753048, 0.0000049 , 0.00000029, 0. , 0. , -0.00001898, -0. , 0.00000001, 0. , -0. , 0. , -0.00002307, -0.00000003, -0.00000002], - [ 0.00045208, -0.0004319 , -0.00047141, 0.00000001, 0. , -0. , -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000085, -0.0000332 , 0.00000294, -0.08085838, -0.11533339, 0.18863828, -0.12538688, -0.04543211, -0.21908948, 0.0732591 , 0.0000114 , 0.1626865 , 0.13183761, -0.18857432], - [ 0.00045208, 0.0006242 , -0.00013833, 0.00000001, 0. , 0. , -0.00000199, 0.00000148, 0.00000002, 0.09959752, 0.00000083, 0.00001915, 0.00002731, 0.14031084, -0.01235872, 0.18865447, 0.02345219, 0.13126238, 0.17298903, 0.15310751, 0.00004989, 0.16300146, -0.2290385 , -0.01949919], - [ 0.00045208, -0.0001923 , 0.00060974, 0.00000001, -0. , -0. , -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000082, 0.00001406, -0.00003026, -0.05945245, 0.12769211, 0.18857909, 0.10205461, -0.08598325, 0.04610051, -0.22636662, -0.0000613 , 0.16239009, 0.09782251, 0.20849111], - [ 0.01733071, -0.01655722, -0.01807177, -0.00000025, -0.00000006, 0.00000004, -0.00007611, 0.0000568 , -0.00000929, 0.05502888, -0.00002725, 0.00000581, -0.00000051, -0.04467114, -0.06371726, 0.0985865 , -0.06553283, -0.02374486, -0.1145189 , 0.03829281, 0.00000596, 0.15007133, 0.12158971, -0.17391612], - [ 0.01733071, 0.02392922, -0.00530308, -0.00000025, -0. , -0.00000007, -0.00007611, 0.0000568 , -0.00000929, 0.05502888, -0.00002726, -0.00000335, -0.00000477, 0.07751634, -0.00682772, 0.09859496, 0.01225717, 0.06860363, 0.09042202, 0.08002988, 0.00002608, 0.15036181, -0.21123494, -0.01798346], - [ 0.01733071, -0.007372 , 0.02337486, -0.00000025, 0.00000006, 0.00000004, -0.00007611, 0.0000568 , -0.00000929, 0.05502888, -0.00002726, -0.00000247, 0.00000528, -0.03284519, 0.07054497, 0.09855556, 0.05333833, -0.04493872, 0.02409691, -0.11832269, -0.00003204, 0.14979796, 0.09021866, 0.1922848 ], - [ 0.0005754 , -0.00024476, 0.00077607, 0. , -0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, 0.00000013, -0.00000071, 0.00000154, 0.00348591, -0.00748706, -0.00381672, -0.00206623, 0.00174084, -0.00093321, 0.00458235, 0.00000124, -0.01384984, -0.00834291, -0.0177814 ], - [ 0.0005754 , -0.00054972, -0.0006 , 0. , 0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, 0.00000013, 0.00000169, -0.00000015, 0.00474102, 0.00676242, -0.00381792, 0.00253862, 0.00091983, 0.00443504, -0.00148299, -0.00000023, -0.01387512, -0.01124392, 0.01608278], - [ 0. , 0. , -0. , -0. , 0.00007956, 0.00004543, -0. , -0. , -0. , -0. , -0.00000072, 0.00464234, -0.0099769 , 0.00000124, -0.00000265, -0. , 0.00000025, -0.00000021, -0.00000005, 0.00000027, 0. , -0. , 0.00000001, 0.00000003], - [ 0.0005754 , 0.00079447, -0.00017607, 0. , 0. , 0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, 0.00000013, -0.00000097, -0.00000139, -0.00822694, 0.00072464, -0.00381825, -0.00047482, -0.00265758, -0.00350182, -0.00309936, -0.00000101, -0.01390199, 0.01953381, 0.00166301], - [-0. , 0. , 0. , -0.00000001, -0.00004543, 0.00007956, -0. , -0. , 0. , 0. , 0.00000182, -0.00997689, -0.00464234, 0.00000265, 0.00000123, -0. , -0.00000021, -0.00000025, 0.00000027, 0.00000005, 0. , -0. , 0.00000003, -0.00000001], - [ 0.28149775, 0.24176118, -0.76656657, -0.00000019, 0.000001 , 0.00000057, -0.35428689, 0.34524753, 0.09633267, -0.0729876 , -0.00000073, 0.00002076, -0.00004468, -0.08713396, 0.18714652, -0. , -0. , 0. , 0. , -0. , -0. , 0. , -0. , -0. ], - [ 0.02922689, 0.025 , -0.07926899, -0. , 0.0000001 , 0.00000006, -0.01745596, 0.02545375, -0.16653584, 0.14976831, 0.0000015 , -0.0000426 , 0.00009168, 0.1787895 , -0.3840045 , -0. , -0. , 0. , -0. , 0. , 0. , 0. , -0. , -0. ], - [-0.07630111, -0.06525625, 0.20691189, -0.0000001 , -0.00000027, -0.00000015, 0.00631767, -0.11557807, -0.04821647, 0.00907678, 0.00000009, -0.00000261, 0.00000561, 0.01095183, -0.02352237, 0. , 0. , -0. , 0. , -0. , -0. , -0.00000001, 0.00000001, 0.00000001], - [ 0. , -0. , -0. , 0.00045206, -0.00055216, 0.00032234, -0. , -0. , 0. , 0.00000098, -0.09957426, -0.14032601, 0.01239032, 0.00001928, 0.00002749, -0.00000216, -0.00000041, -0.0000078 , -0.08460097, -0.2530379 , 0.1885492 , -0.16268774, -0.13183686, 0.18857324], - [ 0. , 0. , -0. , 0.00045213, -0.00000304, -0.0006393 , -0. , -0. , 0. , 0.00000098, -0.09961472, 0.08086766, 0.1153138 , -0.00003348, 0.00000296, -0.00000219, 0.00000764, 0.00000173, 0.17680287, -0.19969978, -0.1886718 , -0.1630027 , 0.22903719, 0.01949908], - [-0. , 0. , -0. , -0.00045206, -0.0005552 , -0.00031706, 0. , 0. , 0. , -0.00000097, 0.09960744, -0.05940465, 0.12770615, -0.00001419, 0.00003048, 0.18871586, -0.20398112, 0.17180667, 0.00000148, -0.00000708, -0. , -0.16239242, -0.09782261, -0.20849132], - [ 0. , -0.00000002, -0.00000003, 0.01733016, -0.0211675 , 0.01235718, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05500797, -0.07752052, 0.0068448 , 0.00001065, 0.00001519, -0.00001467, -0.00001823, -0.00001055, -0.04426999, -0.13240966, 0.09866402, -0.15004811, -0.12159439, 0.17392231], - [ 0. , 0.00000003, -0.00000001, 0.01733289, -0.00011665, -0.02450826, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05503032, 0.04467385, 0.06370299, -0.0000185 , 0.00000164, -0.00001468, 0.00000731, 0.00001976, 0.0925174 , -0.10449889, -0.09872817, -0.1503386 , 0.21124272, 0.01798381], - [ 0.00000004, 0.00000003, -0.00000008, -0.01732993, -0.02128398, -0.01215483, 0.00000002, 0.00000004, -0. , -0.00000058, 0.05504905, -0.0328306 , 0.07057814, -0.00000784, 0.00001684, 0.09874795, -0.10673569, 0.08990001, 0.00000076, -0.00000366, -0. , -0.14978715, -0.09022915, -0.19230721], - [ 0. , -0. , 0. , 0.00057538, 0.00070666, 0.00040356, -0. , -0. , -0. , -0.00000006, 0.0058406 , -0.00348326, 0.0074882 , -0.00000083, 0.00000179, 0.003821 , -0.00413009, 0.00347864, 0.00000003, -0.00000014, -0. , -0.01384994, -0.00834298, -0.01778157], - [-0. , 0. , 0. , -0.00057538, 0.00070278, -0.00041027, 0. , 0. , 0. , 0.00000006, -0.00583821, -0.00822755, 0.00072647, 0.00000113, 0.00000161, -0.00000048, -0.00000058, -0.00000037, -0.00171288, -0.00512316, 0.00381749, -0.01387487, -0.01124375, 0.01608252], - [ 0.00008638, -0.00000918, 0.00002912, -0. , -0. , -0. , -0.00000039, 0.00000028, 0.00000066, 0.01037336, 0.0000001 , 0.00000037, -0.00000079, -0.00154829, 0.00332542, 0. , -0. , 0. , 0. , -0. , -0. , -0. , -0. , -0. ], - [-0. , -0. , 0. , -0.00057547, 0.00000387, 0.0008137 , 0. , 0. , 0. , 0.00000006, -0.00584058, 0.00474141, 0.00676104, -0.00000196, 0.00000017, -0.00000048, 0.00000026, 0.00000064, 0.00357966, -0.00404325, -0.00381997, -0.01390173, 0.0195335 , 0.00166298], - [-0. , 0.00008737, 0.00002756, 0. , -0. , 0. , -0. , 0. , 0. , -0. , -0. , 0.00000238, 0.00000111, 0.00997538, 0.00464446, 0. , -0. , -0. , 0. , 0. , 0. , -0. , 0. , -0. ], - [-0.00036519, 0.00000035, -0.0000011 , 0.28421093, -0.69800517, -0.39855274, -0.34537609, -0.3520123 , 0.09633286, -0.00000217, 0.072973 , 0.08710356, -0.18716625, 0.0000208 , -0.00004466, 0.00000193, 0.00000417, -0.00000352, -0. , 0. , 0. , -0.00000013, 0.00000016, 0.00000034], - [ 0.00001058, 0.00000004, -0.00000011, 0.02938968, -0.07217921, -0.04121348, -0.01676482, -0.02573314, -0.16652394, 0.0000011 , -0.14973273, -0.1787271 , 0.38404494, -0.00004268, 0.00009163, -0.0000083 , -0.00001798, 0.00001515, 0. , -0. , 0. , -0.00000058, 0.00000069, 0.00000147], - [-0.00036278, -0.00000009, 0.0000003 , -0.07671439, 0.18840574, 0.10757746, 0.00344358, 0.11543056, -0.04833963, 0.00008899, -0.00917512, -0.01095181, 0.02353301, -0.00000262, 0.00000562, 0.00000632, 0.00001368, -0.00001153, -0. , 0. , -0. , 0.00000772, -0.00000922, -0.00001969], - [-0.00045208, 0.0004319 , 0.00047141, 0. , -0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000092, 0.00003358, -0.00000271, 0.08085838, 0.11533339, -0.00007569, 0.07859684, -0.21723151, 0.04229099, 0.12644661, 0.18866153, 0.16268795, 0.13183849, -0.18857247], - [-0.00045208, -0.0006242 , 0.00013833, 0. , -0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959752, -0.0000009 , -0.00001957, -0.00002745, -0.14031084, 0.01235872, -0.00005946, 0.22742656, -0.04054806, -0.0883944 , 0.09993416, -0.18860023, 0.16300291, -0.22903763, -0.01949734], - [ 0.00045208, -0.0001923 , 0.00060974, -0.00000001, 0.00000001, 0.00000001, -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000125, 0.00001419, -0.00003053, -0.05945245, 0.12769211, -0.1885797 , -0.10205394, 0.08598269, -0.04610051, 0.22636662, 0.0000613 , 0.16239299, 0.09782165, 0.20848927], - [-0.01733075, 0.01655727, 0.01807176, -0.00000001, 0.00000016, 0.00000029, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000167, 0.00002485, 0.00003476, 0.04467113, 0.06371725, -0.00004099, 0.04107512, -0.11353227, 0.02210568, 0.06609417, 0.09861411, 0.15005196, 0.12157799, -0.17394032], - [-0.01733075, -0.02392926, 0.00530305, -0.00000001, 0.00000033, -0.00000001, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000169, -0.0000426 , 0.00000337, -0.07751633, 0.00682773, -0.00003251, 0.11886017, -0.02118961, -0.04620412, 0.052236 , -0.09858207, 0.15034244, -0.21124606, -0.01800794], - [ 0.01733079, -0.00737202, 0.02337491, 0.0000002 , -0.00000017, -0.0000001 , -0.00007611, 0.00005685, -0.00000929, 0.05502879, 0.00003277, 0.00001547, -0.00003327, -0.03284518, 0.07054496, -0.09855875, -0.05333488, 0.04493581, -0.02409691, 0.11832269, 0.00003204, 0.14975912, 0.09022976, 0.19230857], - [ 0.0005754 , -0.00024476, 0.00077607, -0. , -0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, -0.00000038, -0.00000087, 0.00000187, 0.00348591, -0.00748706, 0.00381761, 0.00206527, -0.00174004, 0.00093321, -0.00458235, -0.00000124, -0.0138498 , -0.00834292, -0.01778143], - [-0.0005754 , 0.00054972, 0.0006 , -0. , -0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, -0.00000208, -0.00000013, -0.00474102, -0.00676242, 0.00000197, -0.00159035, 0.00439734, -0.0008561 , -0.00255966, -0.00381908, -0.0138751 , -0.01124391, 0.0160828 ], - [ 0. , 0. , -0. , 0.00008638, 0.00002652, 0.00001515, -0. , -0. , -0. , 0.0000001 , -0.01037499, 0.0015456 , -0.00332574, 0.0000005 , -0.00000107, -0.00000031, 0.00000008, -0.00000007, 0.00000005, -0.00000027, -0. , 0.00000003, 0. , 0.00000001], - [-0.0005754 , -0.00079447, 0.00017607, -0. , -0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000144, 0.00000151, 0.00822694, -0.00072464, 0.00000164, -0.0046036 , 0.00082015, 0.00178937, -0.00202297, 0.00381784, -0.01390196, 0.01953382, 0.00166304], - [ 0. , 0. , 0. , 0.00000001, 0.00004543, -0.00007956, 0. , 0. , -0. , -0. , -0.00000182, 0.0099769 , 0.00464234, -0.00000211, -0.00000098, 0. , 0.00000021, 0.00000025, 0.00000009, 0.00000002, 0.00000026, 0. , -0.00000003, 0.00000001], - [ 0.2814978 , -0.78474669, 0.17391196, -0.00000018, -0.00000001, -0.00000116, -0.35428689, 0.34524753, 0.09633267, -0.0729876 , -0.00000073, 0.00002828, 0.00004033, 0.20564062, -0.01811304, -0. , -0. , -0. , 0. , 0. , 0. , 0. , 0. , 0. ], - [ 0.0292269 , -0.08114896, 0.01798386, -0. , -0. , -0.00000012, -0.01745596, 0.02545375, -0.16653584, 0.14976832, 0.00000149, -0.00005802, -0.00008276, -0.4219524 , 0.037166 , -0. , -0. , -0. , -0. , -0. , -0. , 0. , 0. , 0. ], - [-0.07630112, 0.21181907, -0.04694237, -0.0000001 , 0. , 0.00000031, 0.00631767, -0.11557807, -0.04821647, 0.00907678, 0.00000009, -0.00000355, -0.00000507, -0.02584689, 0.00227662, 0. , 0. , 0. , 0. , 0. , -0. , -0.00000001, -0.00000002, -0. ], - [ 0. , -0. , -0. , 0.00045206, -0.00055216, 0.00032234, -0. , -0. , 0. , 0.00000098, -0.09957426, -0.14032601, 0.01239032, 0.00001928, 0.00002749, -0.00000218, -0.00000537, 0.00000572, 0.08460097, 0.2530379 , -0.1885492 , -0.16268774, -0.13183687, 0.18857324], - [-0. , -0. , 0. , -0.00045213, 0.00000304, 0.0006393 , 0. , 0. , 0. , -0.00000098, 0.09961384, -0.08086694, -0.11531278, 0.00003348, -0.00000296, 0.1885651 , -0.04678259, -0.26266728, 0.00000541, 0.00000477, 0.00000003, -0.16300379, 0.22903872, 0.0194992 ], - [ 0. , -0. , 0. , 0.00045206, 0.0005552 , 0.00031707, -0. , -0. , 0. , 0.00000098, -0.09960833, 0.05940517, -0.12770728, 0.00001419, -0.00003048, -0.00000219, 0.00000702, 0.00000355, 0.26138042, 0.05318012, 0.18864902, -0.16239133, -0.09782196, -0.20848992], - [ 0. , -0.00000002, -0.00000003, 0.01733016, -0.02116754, 0.0123571 , -0.00000002, -0.00000002, -0. , 0.00000054, -0.05500797, -0.07752052, 0.0068448 , 0.00001065, 0.00001519, -0.00001467, -0.00002079, -0.00000356, 0.04426999, 0.13240966, -0.09866402, -0.15004811, -0.12159384, 0.17392269], - [ 0.00000004, -0.00000009, 0.00000002, -0.01733268, 0.0001167 , 0.02450793, 0.00000002, 0.00000004, -0. , -0.00000059, 0.05505259, -0.04469196, -0.06372882, 0.0000185 , -0.00000164, 0.09866906, -0.02447958, -0.13744397, 0.0000028 , 0.00000247, 0.00000001, -0.15035106, 0.21125963, 0.01798558], - [ 0. , -0.00000001, 0.00000003, 0.01733014, 0.02128429, 0.01215495, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05502679, 0.03281729, -0.07054953, 0.00000784, -0.00001684, -0.00001469, 0.00001828, -0.0000105 , 0.13677514, 0.0278281 , 0.09871625, -0.14977473, -0.09022162, -0.19229194], - [-0. , 0. , -0. , -0.00057538, -0.00070666, -0.00040356, 0. , 0. , 0. , 0.00000006, -0.0058402 , 0.00348303, -0.00748769, 0.00000083, -0.00000179, -0.00000048, 0.00000061, -0.00000032, 0.00529207, 0.00107672, 0.00381951, -0.01384959, -0.00834277, -0.01778112], - [-0. , 0. , 0. , -0.00057538, 0.00070278, -0.00041027, 0. , 0. , 0. , 0.00000006, -0.00583821, -0.00822755, 0.00072647, 0.00000113, 0.00000161, -0.00000048, -0.00000068, -0.00000009, 0.00171288, 0.00512316, -0.00381749, -0.01387487, -0.01124374, 0.01608252], - [-0.00004319, -0.00003207, -0.00007416, 0. , 0. , -0. , 0.0000002 , -0.00000014, -0.00000033, -0.00518668, -0.00000005, -0.00000211, 0.00000095, -0.00266314, -0.00933172, -0. , 0. , 0. , -0. , 0. , 0. , 0. , -0. , 0. ], - [ 0. , 0. , -0. , 0.00057547, -0.00000387, -0.0008137 , -0. , -0. , -0. , -0.00000006, 0.00584098, -0.00474173, -0.0067615 , 0.00000196, -0.00000017, 0.00381795, -0.00094723, -0.00531833, 0.00000011, 0.0000001 , 0. , -0.01390208, 0.01953399, 0.00166303], - [-0.0000748 , -0.00001591, 0.00005044, 0. , -0. , -0. , 0.00000034, -0.00000024, -0.00000057, -0.00898359, -0.00000009, 0.00000064, -0.00000137, -0.00268176, 0.00575931, -0. , -0. , 0. , 0. , -0. , -0. , 0. , -0. , -0. ], - [-0.00036519, -0.00000113, 0.00000024, 0.28412074, 0.00382786, 0.80379856, -0.34537609, -0.3520123 , 0.09633286, -0.00000217, 0.07296363, 0.1185591 , 0.16900692, -0.00004906, 0.00000436, 0.00000193, 0.00000096, 0.00000537, -0. , -0. , -0. , -0.00000013, -0.00000037, -0.00000003], - [ 0.00001058, -0.00000012, 0.00000003, 0.02938035, 0.00039583, 0.08311908, -0.01676482, -0.02573314, -0.16652394, 0.0000011 , -0.1497135 , -0.24327048, -0.34678397, 0.00010067, -0.00000895, -0.00000832, -0.00000413, -0.00002314, 0. , 0. , 0. , -0.00000057, -0.00000162, -0.00000014], - [-0.00036278, 0.0000003 , -0.00000007, -0.07669004, -0.00103322, -0.21696152, 0.00344358, 0.11543056, -0.04833963, 0.00008899, -0.00917395, -0.01490681, -0.02124978, 0.00000617, -0.00000055, 0.00000633, 0.00000314, 0.00001761, -0. , -0. , 0. , 0.00000766, 0.00002168, 0.00000186], - [-0.00045208, 0.0004319 , 0.00047141, 0. , 0. , 0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000092, 0.00003354, -0.00000324, 0.08085838, 0.11533339, 0.00007508, -0.07859766, 0.21723122, -0.04229099, -0.12644661, -0.18866153, 0.16268795, 0.13183558, -0.1885745 ], - [ 0.00045208, 0.0006242 , -0.00013833, -0.00000001, -0. , -0.00000001, -0.00000199, 0.00000148, 0.00000002, 0.09959752, 0.00000126, 0.00001932, 0.00002756, 0.14031084, -0.01235872, -0.18865508, -0.02345204, -0.13126153, -0.17298903, -0.15310751, -0.00004989, 0.16300436, -0.22903647, -0.01949901], - [-0.00045208, 0.0001923 , -0.00060974, 0. , -0. , 0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000089, -0.00001445, 0.00003047, 0.05945245, -0.12769211, 0.00001593, 0.14882962, 0.17668265, -0.1307088 , -0.02667045, 0.18861157, 0.16239154, 0.09782048, 0.20849094], - [-0.01733075, 0.0165572 , 0.01807182, -0.00000001, -0.00000017, -0.00000028, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000167, 0.00001837, -0.00003859, 0.04467115, 0.06371724, 0.0000378 , -0.04107936, 0.11353074, -0.02210568, -0.06609417, -0.09861411, 0.15005189, 0.12161643, -0.17391351], - [ 0.01733079, 0.02392927, -0.0053031 , 0.0000002 , 0. , 0.0000002 , -0.00007611, 0.00005685, -0.00000929, 0.05502879, 0.00003277, 0.00002107, 0.00003004, 0.07751632, -0.00682772, -0.09859815, -0.01225638, -0.06859919, -0.09042202, -0.08002988, -0.00002608, 0.15032303, -0.21126119, -0.01798573], - [-0.01733075, 0.00737197, -0.0233749 , -0.00000001, 0.00000017, -0.00000028, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000169, -0.00004253, 0.00000416, 0.03284519, -0.07054496, 0.00000689, 0.07778429, 0.09233825, -0.06832204, -0.01394076, 0.098588 , 0.14977852, 0.09024543, 0.19228674], - [-0.0005754 , 0.00024476, -0.00077607, -0. , -0. , 0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000116, -0.00000174, -0.00348591, 0.00748706, 0.00000011, -0.00301303, -0.00357596, 0.00264594, 0.00053989, -0.00381807, -0.01384982, -0.00834294, -0.0177814 ], - [-0.0005754 , 0.00054972, 0.0006 , -0. , 0. , 0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, -0.00000203, 0.0000005 , -0.00474102, -0.00676242, -0.00000108, 0.00159153, -0.00439692, 0.0008561 , 0.00255966, 0.00381908, -0.0138751 , -0.01124395, 0.01608277], - [-0. , 0. , -0. , -0.00004319, -0.00007927, 0.00001564, 0. , 0. , 0. , -0.00000005, 0.00518891, -0.00885413, 0.00397106, 0.00000076, 0.00000195, 0.00000016, -0.00000029, -0. , -0.00000015, -0.00000003, 0.00000022, -0.00000001, 0.00000001, -0.00000002], - [ 0.0005754 , 0.00079447, -0.00017607, -0. , 0. , 0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, -0.00000038, -0.00000118, -0.00000169, -0.00822694, 0.00072464, 0.00381913, 0.0004746 , 0.00265635, 0.00350182, 0.00309936, 0.00000101, -0.01390194, 0.01953384, 0.00166301], - [-0. , 0. , 0. , -0.00007481, 0.00004593, 0.00002622, 0. , 0. , 0. , -0.00000009, 0.00898438, 0.00268186, -0.00576007, 0.00000091, -0.00000125, 0.00000027, 0.00000015, -0.00000012, -0.00000015, -0.00000019, -0.00000013, -0.00000002, 0.00000001, 0.00000001], - [ 0.28149777, 0.54298554, 0.59265459, -0.00000018, -0.000001 , 0.00000059, -0.35428689, 0.34524753, 0.09633267, -0.0729876 , -0.00000069, -0.00004904, 0.00000435, -0.11850666, -0.16903348, -0. , 0. , 0. , -0. , 0. , -0. , 0. , -0. , 0. ], - [ 0.02922689, 0.05614896, 0.06128513, -0. , -0.0000001 , 0.00000006, -0.01745596, 0.02545375, -0.16653584, 0.14976831, 0.00000142, 0.00010062, -0.00000892, 0.24316291, 0.3468385 , -0. , 0. , 0. , 0. , -0. , 0. , 0. , -0. , 0. ], - [-0.07630111, -0.14656283, -0.15996951, -0.0000001 , 0.00000027, -0.00000016, 0.00631767, -0.11557807, -0.04821647, 0.00907678, 0.00000009, 0.00000616, -0.00000055, 0.01489506, 0.02124575, 0. , -0. , -0. , -0. , 0. , -0. , -0.00000001, 0.00000001, -0.00000001], - [-0. , 0. , 0. , -0.00045206, 0.00055216, -0.00032234, 0. , 0. , 0. , -0.00000098, 0.09957338, 0.14032477, -0.01239021, -0.00001928, -0.00002749, 0.18859754, 0.25088363, 0.09070763, -0.00000683, 0.0000023 , -0.00000003, -0.16268883, -0.13183775, 0.18857451], - [ 0. , 0. , -0. , 0.00045213, -0.00000304, -0.0006393 , -0. , -0. , 0. , 0.00000098, -0.09961472, 0.08086766, 0.1153138 , -0.00003348, 0.00000296, -0.00000215, -0.00000657, 0.00000431, -0.17680287, 0.19969978, 0.1886718 , -0.1630027 , 0.22903719, 0.01949907], - [ 0. , -0. , 0. , 0.00045206, 0.0005552 , 0.00031706, -0. , -0. , 0. , 0.00000098, -0.09960833, 0.05940517, -0.12770728, 0.00001419, -0.00003048, -0.00000215, -0.00000233, -0.0000075 , -0.26138042, -0.05318012, -0.18864902, -0.16239133, -0.09782195, -0.20848992], - [ 0.00000004, 0.00000006, 0.00000006, -0.01732996, 0.02116724, -0.01235698, 0.00000002, 0.00000004, -0. , -0.00000059, 0.05503023, 0.07755195, -0.00684758, -0.00001066, -0.00001519, 0.09868604, 0.13127803, 0.04746391, -0.00000353, 0.00000119, -0.00000002, -0.15006055, -0.12160387, 0.17393644], - [ 0. , 0.00000003, -0.00000001, 0.01733289, -0.00011674, -0.02450826, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05503032, 0.04467385, 0.06370299, -0.0000185 , 0.00000164, -0.00001465, -0.00000004, 0.00002109, -0.0925174 , 0.10449889, 0.09872817, -0.1503386 , 0.21124267, 0.01798447], - [ 0. , -0.00000001, 0.00000003, 0.01733014, 0.02128424, 0.01215503, -0.00000002, -0.00000002, -0. , 0.00000054, -0.05502679, 0.03281729, -0.07054953, 0.00000784, -0.00001684, -0.00001467, 0.00001344, -0.00001622, -0.13677514, -0.02782809, -0.09871625, -0.14977473, -0.09022222, -0.19229166], - [-0. , 0. , -0. , -0.00057538, -0.00070666, -0.00040356, 0. , 0. , 0. , 0.00000006, -0.0058402 , 0.00348303, -0.00748769, 0.00000083, -0.00000179, -0.00000048, 0.00000042, -0.00000055, -0.00529207, -0.00107672, -0.00381951, -0.01384959, -0.00834278, -0.01778112], - [ 0. , -0. , -0. , 0.00057538, -0.00070278, 0.00041027, -0. , -0. , -0. , -0.00000006, 0.00583861, 0.00822811, -0.00072651, -0.00000113, -0.00000161, 0.00381861, 0.00507974, 0.00183659, -0.00000014, 0.00000005, -0. , -0.01387522, -0.01124403, 0.01608293], - [-0.00004319, 0.00006881, -0.00004234, 0. , 0. , 0. , 0.0000002 , -0.00000014, -0.00000033, -0.00518668, -0.00000005, 0.00000064, 0.00000223, 0.00885569, -0.00396864, -0. , -0. , -0. , 0. , 0. , -0. , 0. , 0. , 0. ], - [-0. , -0. , 0. , -0.00057547, 0.00000387, 0.0008137 , 0. , 0. , 0. , 0.00000006, -0.00584058, 0.00474141, 0.00676104, -0.00000196, 0.00000017, -0.00000048, -0.00000003, 0.00000069, -0.00357966, 0.00404325, 0.00381997, -0.01390173, 0.0195335 , 0.00166299], - [ 0.0000748 , 0.00001591, -0.00005044, -0. , 0. , 0. , -0.00000034, 0.00000024, 0.00000057, 0.00898359, 0.00000009, -0.00000064, 0.00000137, 0.00268132, -0.00575952, 0. , 0. , -0. , 0. , 0. , -0. , -0. , 0. , 0. ], - [-0.00036519, 0.00000078, 0.00000085, 0.28421012, 0.69417811, -0.40518233, -0.34537609, -0.3520123 , 0.09633286, -0.00000217, 0.07302293, -0.20562368, 0.01816164, 0.00002827, 0.00004031, 0.00000193, -0.00000513, -0.00000186, 0. , -0. , 0. , -0.00000013, 0.00000021, -0.00000031], - [ 0.00001058, 0.00000008, 0.00000009, 0.02938959, 0.07178346, -0.04189903, -0.01676482, -0.02573314, -0.16652394, 0.0000011 , -0.14983518, 0.4219176 , -0.03726573, -0.000058 , -0.00008272, -0.00000831, 0.0000221 , 0.000008 , -0. , 0. , -0. , -0.00000057, 0.00000093, -0.00000133], - [-0.00036278, -0.00000021, -0.00000023, -0.07671417, -0.18737274, 0.10936692, 0.00344358, 0.11543056, -0.04833963, 0.00008899, -0.0091814 , 0.02585372, -0.00228352, -0.00000356, -0.00000507, 0.00000633, -0.00001682, -0.00000609, 0. , -0. , 0. , 0.00000769, -0.00001243, 0.00001785], - [ 0.00045208, -0.0004319 , -0.00047141, -0.00000001, -0.00000001, 0.00000001, -0.00000199, 0.00000148, 0.00000002, 0.09959753, 0.00000128, -0.00003351, 0.00000297, -0.08085838, -0.11533339, -0.18863889, 0.12538606, 0.04543181, 0.21908948, -0.0732591 , -0.0000114 , 0.1626894 , 0.13183645, -0.18857265], - [-0.00045208, -0.0006242 , 0.00013833, 0. , 0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959752, -0.0000009 , -0.00001914, -0.00002775, -0.14031084, 0.01235872, 0.00005884, -0.22742641, 0.04054892, 0.0883944 , -0.09993416, 0.18860023, 0.16300291, -0.22903733, -0.01950086], - [-0.00045208, 0.0001923 , -0.00060974, 0. , 0.00000001, -0.00000001, 0.00000199, -0.00000148, -0.00000002, -0.09959753, -0.00000089, -0.00001397, 0.00003069, 0.05945245, -0.12769211, -0.00001655, -0.14882896, -0.1766832 , 0.1307088 , 0.02667045, -0.18861157, 0.16239154, 0.09782369, 0.20848944], - [ 0.01733079, -0.01655725, -0.01807181, 0.0000002 , 0.00000017, -0.0000001 , -0.00007611, 0.00005685, -0.00000929, 0.05502879, 0.00003279, -0.00003652, 0.00000323, -0.04467114, -0.06371725, -0.09858969, 0.06552859, 0.02374333, 0.1145189 , -0.03829281, -0.00000596, 0.15003252, 0.1216047 , -0.17393772], - [-0.01733075, -0.02392924, 0.00530313, -0.00000001, -0.00000033, -0.00000001, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000168, 0.00001768, -0.00003892, -0.07751633, 0.00682771, 0.00002932, -0.11885938, 0.02119405, 0.04620412, -0.052236 , 0.09858207, 0.1503424 , -0.21125006, -0.01796125], - [-0.01733075, 0.00737206, -0.02337487, -0.00000002, -0.00000015, 0.00000029, 0.00007611, -0.00005682, 0.00000929, -0.05502883, 0.00000168, 0.00002423, 0.00003522, 0.03284518, -0.07054497, -0.00001008, -0.07778084, -0.09234115, 0.06832204, 0.01394076, -0.098588 , 0.14977855, 0.09020299, 0.19230662], - [-0.0005754 , 0.00024476, -0.00077607, -0. , 0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000058, -0.00000201, -0.00348591, 0.00748706, 0.00000077, 0.00301208, 0.00357676, -0.00264594, -0.00053989, 0.00381807, -0.01384982, -0.00834289, -0.01778143], - [ 0.0005754 , -0.00054972, -0.0006 , -0. , 0. , -0. , -0.00000254, 0.00000188, 0.00000057, -0.00584035, -0.00000038, 0.00000205, -0.00000018, 0.00474102, 0.00676242, 0.0038188 , -0.00253745, -0.00091941, -0.00443504, 0.00148299, 0.00000023, -0.01387508, -0.01124394, 0.0160828 ], - [-0. , 0. , -0. , -0.00004318, -0.00002681, -0.00007622, 0. , 0. , 0. , -0.00000005, 0.0051868 , 0.00266619, 0.00933158, -0.00000199, 0.00000068, 0.00000016, -0.00000005, 0.00000029, 0.00000015, 0.00000003, -0.00000022, -0.00000001, -0.00000002, -0.00000001], - [-0.0005754 , -0.00079447, 0.00017607, -0. , 0. , -0. , 0.00000254, -0.00000188, -0.00000057, 0.00584035, -0.00000001, 0.00000093, 0.00000187, 0.00822694, -0.00072464, -0.00000076, 0.00460338, -0.00082138, -0.00178937, 0.00202297, -0.00381784, -0.01390196, 0.01953383, 0.00166298], - [ 0. , 0. , 0. , 0.00007481, -0.00004593, -0.00002622, -0. , -0. , -0. , 0.00000009, -0.00898438, -0.00268186, 0.00576007, -0.00000037, 0.0000015 , -0.00000027, -0.00000015, 0.00000012, -0.00000021, 0.00000012, -0.00000013, 0.00000002, -0.00000001, -0.00000001], - ]) - dm0 = 2 * mocc @ mocc.T - dm0 = dm0[None,:,:] - - mf.kernel(dm0 = dm0) + mf.kernel() dm = mf.make_rdm1() if dm.ndim == 2: dm = dm[None,:,:] test_gradient = get_vxc_full_response(mf._numint, cell, mf.grids, mf.xc, dm, kpts, hermi=1) - # ref_gradient = numerical_gradient_exc_becke(cell, "LDA0", [1,1,1], 'def2-universal-jkfit', (50,194), dm) - ref_gradient = np.array([ - [ 0.0000000095035091, 0.0000000094857455, 0.0000000095390362], - [-0.0002934405785027, -0.0002934406140298, -0.0002934405607391], - [-0.0000000263788991, -0.0000000263611355, 0.0000000984989867], - [ 0.0002934578446911, 0.0002934579157454, -0.0002935301246509], - [-0.0000000261657362, 0.0000000985878046, -0.0000000263788991], - [ 0.0002934578979819, -0.0002935301424145, 0.0002934578802183], - [ 0.0000000986233317, -0.0000000261834998, -0.0000000262190269], - [-0.0002935301779416, 0.0002934578802183, 0.0002934578979819], - [ 2.1212388574554097, 0.0000000038546943, 0.0000000030730973], - [ 0.0000000035527137, 2.1212388573843555, 0.000000003126388 ], - [ 0.0000000037658765, 0.0000000031974423, 2.121238857029084 ], - ]) + ref_gradient = numerical_gradient_exc_becke(cell, "LDA0", [1,1,1], 'def2-universal-jkfit', (50,194), dm) - # It can match down to 1e-9, if the finite difference is computed using the same dm from SCF. - # However if we save the finite difference result, it suffers from the numerical instability of dm, and the a 3e-7 error is observed. - assert np.max(np.abs(test_gradient - ref_gradient)) < 5e-8 + assert np.max(np.abs(test_gradient - ref_gradient)) < 2e-8 def test_xc_gradient_gga_with_response(self): cell = pyscf.M(