Hessian failure for heavy element (Zn)

[cvsik]

Hi,

I have problems computing the Hessian for a system containing Zn (dummy system in the example given here):

from pyscf import scf, gto

mol = gto.M(atom="Zn 0 0 0; H 1 0 0; H 0 1 0; H 1 1 0; H 0 1 1", basis="def2-svp")

mf = scf.RKS(mol).density_fit()
mf.xc = "b3lyp"

mf.kernel()

hess = mf.Hessian()
hess.auxbasis_response = 2
result = hess.kernel()
print(result)

mf_gpu = mf.to_gpu()
mf_gpu.kernel()

hess_gpu = mf_gpu.Hessian()
hess_gpu.auxbasis_response = 2
result_gpu = hess_gpu.kernel()
print(result_gpu)

The Hessian evaluates fine on the CPU, but fails on the GPU with the following error:

Traceback (most recent call last):
  File "/fs/home/cvsik/Projects/transition_states/debug/example.py", line 20, in <module>
    result_gpu = hess_gpu.kernel()
                 ^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py", line 803, in kernel
    de = hessobj.hess_elec(mo_energy, mo_coeff, mo_occ, atmlst=atmlst)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py", line 80, in hess_elec
    mo1, mo_e1 = hessobj.solve_mo1(mo_energy, mo_coeff, mo_occ, h1mo,
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py", line 945, in solve_mo1
    return solve_mo1(self.base, mo_energy, mo_coeff, mo_occ, h1mo,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py", line 707, in solve_mo1
    raw_mo1 = krylov(fvind_vo, mo1.reshape(-1,nmo*nocc),
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/lib/cupy_helper.py", line 632, in krylov
    axt = aop(x1)
          ^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py", line 661, in fvind_vo
    v = fx(mo1).reshape(-1,nmo, nocc)
        ^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rhf.py", line 739, in fx
    return hessobj.get_veff_resp_mo(mol, dm1, mo_coeff, mo_occ, hermi=1)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/hessian/rks.py", line 832, in get_veff_resp_mo
    grids.build(mol=mol, with_non0tab=False, sort_grids=True)
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/dft/gen_grid.py", line 502, in build
    atom_grids_tab = self.gen_atomic_grids(
                     ^^^^^^^^^^^^^^^^^^^^^^
  File "<string>", line 6, in gen_atomic_grids
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/dft/gen_grid.py", line 243, in gen_atomic_grids
    angs = prune(chg, rad, n_ang)
           ^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/shared/conda/envs/tstoolkit-dev3/lib/python3.11/site-packages/gpu4pyscf/dft/gen_grid.py", line 80, in sg1_prune
    r_atom = radii[nuc] + 1e-200
             ~~~~~^^^^^
IndexError: index 30 is out of bounds for axis 0 with size 19

I'm using gpu4pyscf==1.3.0 and pyscf==2.7.0
Looks like grid pruning fails because there are not enough radii available 😀

Any help on this would be much appreciated! Thanks!

Best,
Max

[wxj6000]

@cvsik Thank you for the feedback. The issue is caused by switching the default grids of CPSCF to SG1 grids. The original SG1 grids do not support heavy atoms after Ar. The issue has been fixed by turning off the pruning scheme for heavy atoms. (#330)

Please note different grids are used in PySCF and GPU4PySCF, which will result in a small discrepancy between the two. If you would like to use the finer grids, please follow this example.
https://github.com/pyscf/gpu4pyscf/blob/master/examples/00-h2o.py#L63