I report on the derivation, development and computer implementation of methods for computing the energies and forces between small rigid polarizable molecules, that are defined by the center-of-mass moments of their electronic and nuclear charge distributions and their linear response moments. The formalism is based on compact and efficient storage and manipulation of symmetric Cartesian tensors of arbitrary rank, and a general formula for the Cartesian gradients of one-dimensional interaction (kernel) potentials. The theory is applied to many-body interactions among water molecules. Permanent moments of the water molecule are computed up to the 9th order with quantum-chemistry software and their basis-set dependence is investigated. Response moments up to the 5th order are similarly investigated. Kernel potentials for electronic, nuclear and polarized interactions are suggested and compared to interaction energies from symmetry-adapted perturbation-theory. I discuss vibrational degrees of freedom and report on a novel method for fitting high rank moment tensors to a exible geometry. The method is based on decomposition of the tensor into a sum of outer products of vectors, which are defined in the lab-frame by the molecular geometry. I show that the formalism, which is based on an asymptotic expansion, can
give good results at all ranges.