The generalized energy-based fragmentation (GEBF) method is extended to allow calculations of nuclear magnetic resonance (NMR) chemical shifts of macromolecular and condensed-phase systems feasible at a low computational cost. In this approach, NMR shielding constants in a large system are evaluated as a linear combination of the corresponding quantities from a series of small "electrostatically embedded" subsystems. Comparison of NMR shielding constants from the GEBF-X method [where X is an electronic structure method, such as Hartree-Fock (HF), density functional theory (DFT), ...] with those from the conventional quantum chemistry method for two representative systems verifies that the GEBF approach can reproduce the results of the conventional quantum chemistry method very well. This procedure has further been applied to compute NMR shielding constants of a large foldamer and a supramolecular aggregate, and the 15N shielding constant for CH3CN in the CHCl3 solvent. For the former two systems, the predicted 1H chemical shifts are in good agreement with the experimental data. For the CH3CN/CHCl3 solution, the 15N shielding constant of CH3CN is evaluated as the ensemble average of up to 200 sufficiently large CH3CN/CHCl3 clusters from either classical or QM/MM (quantum mechanics/molecular mechanics) molecular dynamics (MD) simulations. Our results reveal that the gas-to-solution shift of 15N (from an isolated CH3CN to the CH3CN/CHCl3 solution) based on PM6-DH+/MM MD simulation is in good accord with the experimental value, outperforming those based on classical MD simulation and the previous polarizable continuum model using integral equation formalism (IEF-PCM) study. This study unravels that the generation of representative liquid structures is critical in evaluating the NMR shielding constants of condensed-phase systems.