We investigate the possibility of generating initial conditions for cosmological N-body simulations by simulating a system whose correlations at thermal equilibrium approximate well those of cosmological density perturbations. The system is an appropriately modified version of the standard ``one component plasma'' (OCP). We show first how a well-known semi-analytic method can be used to determine the potential required to produce the desired correlations, and then verify our results for some cosmological type spectra with simulations of the full molecular dynamics. The advantage of the method, compared to the standard one, is that it gives by construction an accurate representation of both the real and reciprocal space correlation properties of the theoretical model. Furthermore the distributions are also statistically homogeneous and isotropic. We discuss briefly the modifications needed to implement the method to produce configurations appropriate for large N-body simulations in cosmology, and also the generation of initial velocities in this context.