Based on the thermodynamic perturbation theory (TPT) and the random phase approximation (RPA), we present a statistical theory of solutions of electrically neutral soft molecules, every of which is modelled as a set of sites that interact with each other through the potentials, presented as the sum of the Coulomb potential and arbitrary soft-core potential. As an application of our formalism, we formulate a general statistical theory of solution of the soft-core dipolar particles. For the latter, we obtain a new analytical relation for the screening function. As a special case, we apply this theory to describing the phase behavior of a solution of the dipolar particles interacting with each other in addition to the electrostatic potential through the repulsive Gaussian potential - Gaussian core dipolar model (GCDM). Using the obtained analytic expression for the total free energy of the GCDM, we obtain the liquid-liquid phase separation with an upper critical point. The developed formalism could be used as a general framework for the coarse-grained description of thermodynamic properties of solutions of macromolecules, such as proteins, betaines, polypeptides, etc. (C) 2019 Elsevier B.V. All rights reserved.