The self-assembly and stability of DNA G-quadruplexes (GQs) are affected by the intrinsic stability of different GpG base steps embedded in their G-quartet sterns. We have carried out MD simulations followed by MM-PBSA (molecular mechanics Poisson-Boltzmann surface area) free energy calculations on all the experimentally observed three-quartet intramolecular human telomeric GQ topologies. We also studied antiparallel GQ models with alternative syn-anti patterns of the G-quartets. We tested different ions, dihedral variants of the DNA force field, water models, and simulation lengths. In total, similar to 35 mu s of simulations have been carried out. The systems studied here are considerably more complete than the previously analyzed two-quartet stems. Among other effects, our computations included the stem loop coupling and ion ion interactions inside the stem. The calculations showed a broad agreement with the earlier predictions. However, the increase in the completeness of the system was associated with increased noise of the free energy data which could be related, for example, to-the presence of long-lived loop substates and rather complex dynamics-for the. two bound ions inside the G-stem. As a result, the MM-PBSA data, were noisy and we could not improve their quantitative convergence even by expanding the simulations to 2.5 mu s long trajectories. We also suggest that the quality of MM-based free energy computations based on MD simulations of complete GQs is more sensitive to the neglect of explicit polarization effects, which, in real systems, are associated with the presence of multiple closely spaced-ions inside the GQs. Thus, although the MM-PBSA procedure provides very useful insights that complement the structural-dynamics data from MD trajectories of GQs, the method is far from reaching quantitative accuracy. Our conclusions are in agreement with critical assessments of the MM-PBSA methodology available in contemporary literature for other types of problems.