The exact expressions for the dipole, quadrupole, and octupoles of a collection of N point charges involve summations of corresponding tensors over the N sites weighted by their charge magnitudes. When the point charges are atoms (in a molecule) the N‐site formula is an approximation, and one must integrate over the electron density to recover the exact multipoles. In the present work we revisit the N(N + 1)/2‐site point charge density model of Hall (Chem. Phys. Lett. 6, 501, 1973) for the purpose of fitting ab initio derived multipole moment hypersurfaces using permutationally invariant polynomials (PIP). We examine new approaches in PIP‐fitting procedures for the dipole, quadrupole, octupole moments, and polarizability tensor surfaces (DMS, QMS, OMS and PTS, respectively) for a non‐polar CCl4 and a polar CHCl3 and show that compared to the primitive N‐site model the N(N + 1)/2‐site model appreciably improves the relative RMSE of the DMS and does much more substantially so, by an order of magnitude, for the corresponding ones of QMS and OMS. Training datasets are obtained by sampling potential energies up to 18 000 cm−1 above the global minima, generated by molecular dynamics simulations at the DFT B3LYP/aug‐cc‐pVDZ level of theory.