def METHOD_NAME(mf):
    mol = mf.mol
    mfs = mf.Gradients().set(grid_response=True).as_scanner()
    def grad_full(ia, inc):
        coord = mol.atom_coord(ia).copy()
        ptr = mol._atm[ia,gto.PTR_COORD]
        de = []
        for i in range(3):
            mol._env[ptr+i] = coord[i] + inc
            e1a = mfs(mol.copy())[1]
            mol._env[ptr+i] = coord[i] - inc
            e1b = mfs(mol.copy())[1]
            mol._env[ptr+i] = coord[i]
            de.append((e1a-e1b)/(2*inc))
        return de
    natm = mol.natm
    e2ref = [grad_full(ia, .5e-3) for ia in range(mol.natm)]
    e2ref = numpy.asarray(e2ref).reshape(natm,3,natm,3).transpose(0,2,1,3)
    return e2ref