With the aid of the Heisenberg operators for the electromagnetic field in the neighbourhood of several molecules, the many-body interaction potentials are calculated by finding the energy of a test molecule in this field. In the first instance the calculation is made for the case where the intermolecular separations are very large. In this far-zone range, the energies depend on the molecules through their static polarizabilities. The three-body non-additive dispersion energy is calculated for an arbitrary configuration and the dependence on the geometry is investigated. The four-body non-additive potential is found for the regular tetrahedral configuration. The theory is extended to cover all separations outside regions of overlap. It is shown that the interaction energy in this case depends on the dynamic polarizabilities. The full Casimir-Polder potential follows directly from the general expression. The near-zone limit of the complete formula tends to the form obtained by using electrostatic couplings only: as with, for example, the Axilrod-Teller potential for N = 3. The work described here presents an alternative viewpoint to the conventional perturbation theory and lends further insight into the nature of intermolecular forces.