def METHOD_NAME(mol, basis_name, origin, order, step=1e-4):
    prefactors_5p = np.array([1.0, -8.0, 8.0, -1.0]) / 12.0
    multipliers_5p = [-2, -1, 1, 2]
    natoms = mol.natom()
    basis = psi4.core.BasisSet.build(mol, 'orbital', basis_name)
    nbf = basis.nbf()
    nmul = cumulative_cart_dim(order) - 1
    int_grad = np.zeros((natoms, 3, nmul, nbf, nbf))
    for i in range(natoms):
        for c in range(3):
            for f, p in zip(multipliers_5p, prefactors_5p):
                mol_p = mol.clone()
                coords_p = mol_p.geometry()
                coords_p.np[i, c] += f * step
                mol_p.set_geometry(coords_p)
                basis = psi4.core.BasisSet.build(mol_p, 'orbital', basis_name)
                mints = psi4.core.MintsHelper(basis)
                ints_pert = matlist_to_ndarray(mints.ao_multipoles(order, origin))
                int_grad[i, c, :, :, :] += p * ints_pert / step
    return int_grad