def METHOD_NAME(self):
    axis = np.array([1., 0., 0.])
    dh_k = 2.
    N = 100  # even number to get singularities at phi=0 and phi=pi
    d_phi = 2 * np.pi / N
    for dh_n, dh_phi0_div in [(2, 3), (3, 6)]:
        with self.subTest(multiplicity=dh_n, phi_0=f"pi / {dh_phi0_div}"):
            dh_phi0 = np.pi / dh_phi0_div
            dihedral = espressomd.interactions.Dihedral(
                bend=dh_k, mult=dh_n, phase=dh_phi0)
            self.system.bonded_inter.add(dihedral)
            self.system.part.clear()
            p0, p1, p2, p3 = self.system.part.add(pos=4 * [(0., 0., 0.)])
            p1.add_bond((dihedral, p0, p2, p3))
            p1.pos = [5., 5., 5.]
            p2.pos = p1.pos + [1., 0., 0.]
            p0.pos = p1.pos + [0., 1., 0.]
            for i in range(N):
                phi = i * d_phi
                p3.pos = p2.pos + rotate_vector([0., 1., 0.], axis, phi)
                self.system.integrator.run(recalc_forces=True, steps=0)
                # Calculate expected forces and energies
                E_ref, forces_ref = dihedral_potential_and_forces(
                    dh_k, dh_n, dh_phi0, p0.pos, p1.pos, p2.pos, p3.pos)
                self.check_values(E_ref, forces_ref)
    self.check_undefined_angle()