Simulation of the solid solutions in the system of double sodium–gadolinium and sodium–europium molybdates, which are promising matrices for solid state lasers and phosphors has been carried out by the method of interatomic potentials. Two types of solid solutions have been studied, one of which contains finite components corresponding to the stoichiometric NaGd(MoO4)2–NaEu(MoO4)2 compositions with statistical distribution of cations in the crystal lattice. Another object is a cation-deficient Na2Gd4(MoO4)7–Na2Eu4(MoO4)7 system, in which we have examined the variants of statistical distribution and partial ordering of cations over structural positions. Atomistic simulation has been performed using the GULP 4.0.1 software package (General Utility Lattice Program). It is shown that when we pass from sodium-gadolinium molybdate to sodium-europium molybdate, both of stoichiometric and cation-deficient compositions, an increase in the unit cell volume is observed, while the density of the crystal, the energy of interatomic interactions in the structure, the vibrational entropy and the heat capacity decrease along with increasing europium content. The energy of interatomic interactions in the structure for cation-deficient solid solutions is less than for stoichiometric ones. Other aforementioned characteristics for cation-deficient solid solutions have greater values than for stoichiometric ones. The role of cluster europium centers in concentration quenching in NaGd(MoO4)2–NaEu(MoO4)2 solid solutions has been examined.