def METHOD_NAME(co2, molecule):
    kt, theta0 = 1, np.pi / 4
    calc = Manybody({1: ZeroPair()}, {1: HarmonicAngle(kt, theta0)}, molecule)
    co2.calc = calc
    angle = np.radians(co2.get_angle(0, 1, 2))
    # Testing potential energy
    epot = co2.get_potential_energy()
    epot_ref = 0.5 * kt * (angle - theta0)**2
    nt.assert_allclose(epot, epot_ref, rtol=1e-14,
                       err_msg="Wrong energy")
    # Testing forces
    f = co2.get_forces()
    f_ref = numerical_forces(co2, d=1e-6)
    nt.assert_allclose(f, f_ref, rtol=1e-6, atol=2e-9,
                       err_msg="Wrong forces")
    # Checking zeros
    nt.assert_allclose(np.abs(f.sum()), 0, atol=1e-13)
    # Testing stress
    s = co2.get_stress()
    s_ref = numerical_stress(co2, d=1e-6)
    nt.assert_allclose(s, s_ref, rtol=1e-6, atol=2e-9)
    # Testing nonaffine forces with finite differences
    nf_ref = numerical_nonaffine_forces_reference(co2, d=1e-6)
    nf = co2.calc.get_property('nonaffine_forces', co2)
    nt.assert_allclose(nf, nf_ref, rtol=1e-8, atol=1e-6)
    # Testing hessian
    h = co2.calc.get_property('hessian', co2)
    h_ref = numerical_hessian(co2, d=1e-4)
    # print(f"{calc.get_block_sparse_hessian(co2)[0]}")
    # print(f"{h.todense()}\n\n{h_ref.todense()}")
    nt.assert_allclose(h.todense(), h_ref.todense(), atol=1e-5)