A general numerical method for solving the Wertheim association theory in the case of fully anisotropic polyatomic rigid molecules is proposed. In order to handle the nonspherical shapes of the molecules, the Wertheim association theory is combined with the rotational invariant approach of the molecular Ornstein-Zernike (MOZ) method. Therefore, this approach is called the WMOZ method. It is carried out in the association hypernetted chain (AHNC) approximation which is the only approximation and which consists of the neglect of the bridge functions. The method of solution is implemented by translating the set operations appearing in the Wertheim association theory into an algebraic formalism developed through a one-to-one correspondence between subsets and binary numbers. The Wertheim association theory is explicitly solved for hydrogen fluoride using two interaction site models. For these models, site-site distribution functions, internal excess energies, and dielectric constants obtained by the MOZ and WMOZ methods are compared to the exact molecular dynamics results. The WMOZ method strongly improves the structural description of the H bonding in comparison with the MOZ approximation. The quality of the dielectric and thermodynamic results obtained by the WMOZ approach is discussed in comparison with the simulation values. To our knowledge, this is the first solution of the Wertheim association theory for a liquid of particles of nonspherical shapes.