We report the results of molecular dynamics simulations of the properties of a pseudo-atom model of dodecane thiol ligated 5-nm diameter gold nanoparticles (AuNP) in vacuum as a function of ligand coverage and particle separation in three state of aggregation: the isolated AuNP, an isolated pair of AuNPs and a square assembly of AuNPs. Our calculations show that for all values of the coverage the ligand density along a radius emanating from the core of an isolated AuNP oscillates along the chain up to the fourth pseudo-atom, then smoothly decays to zero. We examine the ligand envelope as a function of the coverage and demonstrate that the deformation of that envelope generated by interaction between the NPs is coverage-dependent, so that the shape, depth and position of the minimum of the potential of mean force displays a systematic dependence on the coverage. We propose an accurate analytical description of the calculated potential of mean force with parameters that scale linearly with the ligand coverage. We define and calculate an effective pair potential of mean force for a square configuration of particles; our definition contains, implicitly, both the three- and four-particle contributions to deviation from additivity. We find that the effective pair potential of mean force in this configuration has a different minimum and a different well depth than the isolated pair potential of mean force. Previous work has found that the three-particle contribution to deviation from additivity is monotone repulsive, whereas we find that the combined three- and four-particle contributions have an attractive well, implying that the three- and four-particle contributions are of comparable magnitude but opposite sign, thereby suggesting that even higher order correction terms likely play a significant role in the behavior of assemblies of many nanoparticles.