def METHOD_NAME(mol_h2o):
    psi4.set_options({'basis': 'cc-pvdz'})
    _, wfn = psi4.energy('HF', molecule=mol_h2o, return_wfn=True)
    order = 8
    ints_grad = fdiff_multipole_integral(mol_h2o, 'cc-pvdz', [1.0, 2.0, 3.0], order=order)
    nmul = cumulative_cart_dim(order) - 1
    grad_fdiff = np.einsum('ncmij,ij->ncm', ints_grad, wfn.Da().np).reshape(-1, nmul)
    mints = psi4.core.MintsHelper(wfn.basisset())
    # test against finite differences
    grad = mints.multipole_grad(D=wfn.Da(), order=order, origin=[1.0, 2.0, 3.0])
    np.testing.assert_allclose(grad_fdiff, grad.np, atol=1e-7)
    # test that we get the same result from the 'hard-wired' dipole_grad
    grad_dip = mints.dipole_grad(wfn.Da())
    grad = mints.multipole_grad(D=wfn.Da(), order=1, origin=[0.0, 0.0, 0.0])
    np.testing.assert_allclose(grad_dip.np, grad.np, atol=1e-14)