Spatial structures of molecular clusters modeling a solvate shell around phosphorus-containing methyl-and butyl-derivatives of phosphine and betaine molecules dissolved in different solvents (acetone, toluene, formamide) have been calculated by using different variants of density functional theory (unrestricted Becke three-parameter Lee-Yang-Parr [UB3LYP], Perdew-Burke-Ernzerhof [PBE], optimized exchange functional [OPTX] developed by Handy and Cohen in conjunction with Lee-Yang-Parr [LYP] correlational functional [OLYP]) with 6-31G(d,p) and 6-31G++(d,p) basis sets. The P-31 magnetic shielding constants for the structures are calculated with the usage of gauge-including atomic orbitals in UB3LYP/6-31G(d, p) and 6-31G++(d, p) methods. The modeling of molecular clusters is done by using the supermolecular model, the molecular mechanics method and the combination of quantum chemistry and molecular mechanics methods (QM/MM). The own N-layered integrated molecular orbital method (ONIOM) has been applied for modeling and calculating of isotropic P-31 nucleus magnetic shielding of clusters of trimethylphosphine and trimethylbetaine molecules dissolved in acetone using combinations of UB3LYP/6-31G(d, p) (higher level) and unrestricted Hartree-Fock (UHF)/6-31G(d, p) (lower level) methods. Applicability of the ONIOM approach and different ways of modeling to the calculation of P-31 nucleus magnetic shielding constants is studied. A comparison of the results obtained by the density functional theory, ONIOM and MM methods is given.