A molecular modeling approach is introduced as a way to treat multibody (more than two molecules) contributions to the intermolecular potential. There are two key features to the method. First, it employs polarizable electrostatics on the molecules, but converges the charges and fields for only three molecules at a time, taken separately for all trimers (three molecules falling within a cutoff distance) in the system. This feature introduces significant computational savings when applied in Monte Carlo simulation (in comparison to a full N-body polarization treatment), as movement of a single molecule does not require re-converging of the polarization of all molecules, and it achieves this without approximations that cause the value of the energy to depend on the history of the simulation. Second, the approach defines the polarization energy in excess of the pairwise contribution, meaning that the trimer energy has subtracted from it the sum of the energies obtained by converging the polarization of each molecule pair in the trimer. This feature is advantageous because it removes the need (often found in polarizable models) to stiffen inappropriately the repulsive part of the pair potential. The polarization contribution is thus a purely three-body potential. The approach is applied to model hydrogen fluoride, which in experiments exhibits unusual properties that have proven difficult to capture well by molecular models. The new HF model is shown to be much more successful than previous modeling efforts in obtaining agreement with a broad range of experimental data (volumetric properties, heat effects, molecular structure, and vapor-liquid equilibria).