diff --git a/pyscf/lib/pbc/nr_direct.c b/pyscf/lib/pbc/nr_direct.c index 0aa088a369..e19f1019a7 100644 --- a/pyscf/lib/pbc/nr_direct.c +++ b/pyscf/lib/pbc/nr_direct.c @@ -1174,7 +1174,7 @@ void PBCVHFnr_sindex(int16_t *sindex, int *atm, int natm, { float fac_guess = .5f - logf(omega2)/4; int ijb, ib, jb, i0, j0, i1, j1, i, j, li, lj; - float dx, dy, dz, ai, aj, ci, cj, aij, a1, fi, fj, rr, rij, dri, drj; + float dx, dy, dz, ai, aj, ci, cj, aij, a1, rr; float log_fac, theta, theta_r, r_guess, v; #pragma omp for schedule(dynamic, 1) for (ijb = 0; ijb < ngroups*(ngroups+1)/2; ijb++) { @@ -1193,8 +1193,6 @@ void PBCVHFnr_sindex(int16_t *sindex, int *atm, int natm, cj = cs[jb]; aij = ai + aj; - fi = ai / aij; - fj = aj / aij; a1 = ai * aj / aij; theta = omega2/(omega2+aij); r_guess = sqrtf(-logf(1e-9f) / (aij * theta)); @@ -1209,10 +1207,7 @@ void PBCVHFnr_sindex(int16_t *sindex, int *atm, int natm, dy = ry[i] - ry[j]; dz = rz[i] - rz[j]; rr = dx * dx + dy * dy + dz * dz; - rij = sqrtf(rr); - dri = fj * rij + theta_r; - drj = fi * rij + theta_r; - v = li*logf(dri) + lj*logf(drj) - a1*rr + log_fac; + v = (li+lj)*logf(MAX(theta_r, 1.f)) - a1*rr + log_fac; sindex[i*Nbas+j] = v * LOG_ADJUST; } } if (ib > jb) { diff --git a/pyscf/lib/vhf/nr_sr_vhf.c b/pyscf/lib/vhf/nr_sr_vhf.c index e0e5e08838..bd4a93fda1 100644 --- a/pyscf/lib/vhf/nr_sr_vhf.c +++ b/pyscf/lib/vhf/nr_sr_vhf.c @@ -888,9 +888,9 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond, float fac_guess = .5f - logf(omega2)/4; int ish, jsh, li, lj; int ij, i, j, di, dj, dij, di2, dj2; - float ai, aj, fi, fj, aij, a1, ci, cj; + float ai, aj, aij, ai_aij, a1, ci, cj; float xi, yi, zi, xj, yj, zj, xij, yij, zij; - float dx, dy, dz, r2, r, dri, drj, v, log_fac, r_guess, theta, theta_r; + float dx, dy, dz, r2, v, log_fac, r_guess, theta, theta_r; double qtmp, tmp; float log_qmax; int shls[4]; @@ -924,12 +924,11 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond, dy = yj - yi; dz = zj - zi; aij = ai + aj; - fi = ai / aij; - fj = aj / aij; - a1 = fi * aj; - xij = xi + fj * dx; - yij = yi + fj * dy; - zij = zi + fj * dz; + ai_aij = ai / aij; + a1 = ai_aij * aj; + xij = xi + ai_aij * dx; + yij = yi + ai_aij * dy; + zij = zi + ai_aij * dz; theta = omega2/(omega2+aij); r_guess = R_GUESS_FAC / sqrtf(aij * theta); @@ -938,10 +937,7 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond, log_fac = logf(ci*cj * sqrtf((2*li+1.f)*(2*lj+1.f))/(4*M_PI)) + 1.5f*logf(M_PI/aij) + fac_guess; r2 = dx * dx + dy * dy + dz * dz; - r = sqrtf(r2); - dri = fj * r + theta_r; - drj = fi * r + theta_r; - v = li*logf(dri) + lj*logf(drj) - a1*r2 + log_fac; + v = (li+lj)*logf(MAX(theta_r, 1.f)) - a1*r2 + log_fac; s_index[ish*Nbas+jsh] = v; s_index[jsh*Nbas+ish] = v; xij_cond[ish*Nbas+jsh] = xij; diff --git a/pyscf/lib/vhf/test/test_nr_direct.py b/pyscf/lib/vhf/test/test_nr_direct.py index ffa208469d..0e5894a25e 100644 --- a/pyscf/lib/vhf/test/test_nr_direct.py +++ b/pyscf/lib/vhf/test/test_nr_direct.py @@ -285,6 +285,40 @@ def test_direct_jk_s2(self): self.assertTrue(numpy.allclose(vj0,vj1)) self.assertTrue(numpy.allclose(vk0,vk1)) + def test_sr_vhf_q_cond(self): + for omega in [.1, .2, .3]: + for l in [0, 1, 2, 3]: + for a in [.15, .5, 2.5]: + rs = numpy.arange(1, 10) * 2. + mol = gto.M(atom=['H 0 0 0'] + [f'H 0 0 {r}' for r in rs], + basis=[[l,[a,1]]], unit='B') + nbas = mol.nbas + q_cond = numpy.empty((6,nbas,nbas), dtype=numpy.float32) + ao_loc = mol.ao_loc + cintopt = lib.c_null_ptr() + with mol.with_short_range_coulomb(omega): + with mol.with_integral_screen(1e-26): + libcvhf2.CVHFnr_sr_int2e_q_cond( + libcvhf2.int2e_sph, cintopt, + q_cond.ctypes.data_as(ctypes.c_void_p), + ao_loc.ctypes.data_as(ctypes.c_void_p), + mol._atm.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.natm), + mol._bas.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(mol.nbas), + mol._env.ctypes.data_as(ctypes.c_void_p)) + + s_index = q_cond[2] + si_0 = s_index[0,0] + si_others = s_index.diagonal()[1:] + with mol.with_short_range_coulomb(omega): + ints = [abs(mol.intor_by_shell('int2e', (0,0,i,i))).max() + for i in range(1, mol.nbas)] + + aij = akl = a * 2 + omega2 = mol.omega**2 + theta = 1/(2/aij+1/omega2) + rr = rs**2 + estimator = rr * numpy.exp(si_0 + si_others - theta*rr) + assert all(estimator / ints > 1) if __name__ == '__main__':