Change DF gradient routine full/short-range j2c inversion algorithm to CD

[henryw7]

Similar to PR #783, for a big system (gpu4pyscf/benchmarks/molecules/organic/095_Azadirachtin.xyz) with def2-universal-jkfit auxbasis, the smallest eigenvalue of the j2c matrix is close to threshold (1e-7 by default), and it introduces numerical instability to gradient calculation as well. Here we noticed that the general matrix solver cp.linalg.solve is not numerically stable enough (when changing GPU type from A100 to A30, a 1e-4 error on final gradient is introduced). Surprisingly, eigenvalue solver is also not stable enough (same magnitude of error as cp.linalg.solve). Only the Cholesky solver provides a reasonable stability (1e-5 error across GPUs). This conclusion might be coincidence on the particular system, and any suggestion is welcome.

This PR also fix the CD part in

gpu4pyscf/gpu4pyscf/df/grad/rhf.py

Line 32 in 3edad9f

def _gen_metric_solver(int2c, decompose_j2c='CD', lindep=df.LINEAR_DEP_THR):