Identifying stable speciation in multicomponent liquid solutions is of fundamental importance to areas ranging from electrochemistry to organic chemistry and biomolecular systems. However, elucidating this complex solvation environment is a daunting task even when using advanced experimental and computational techniques. Here, we introduce a fully automated, high-throughput computational framework for the accurate and robust prediction of stable species present in liquid solutions by computing the nuclear magnetic resonance (NMR) chemical shifts of molecules. The framework automatically extracts and categorizes hundreds of thousands of atomic clusters from classical molecular dynamics (CMD) simulations to identify the most stable speciation in the solution and calculate their NMR chemical shifts via DFT calculations. Additionally, the framework creates an output database of computed chemical shifts for liquid solutions across a wide chemical and parameter space. This task can be infeasible experimentally and challenging using conventional computational methods. To demonstrate the capabilities of our framework, we compare our computational results to experimental measurements for a complex test case of magnesium bis(trifluoromethanesulfonyl)imide Mg(TFSI)2 salt in dimethoxyethane (DME) solvent, which is a common electrolyte system for Mg-based batteries. Our extensive benchmarking and analysis of the Mg2+ solvation structural evolutions reveal critical factors such as the effect of force field parameters that influence the accuracy of NMR chemical shift predictions in liquid solutions. Furthermore, we show how the framework reduces the efforts of performing and managing over 300 13C and 600 1H DFT chemical shift predictions to a single submission procedure. By enabling more efficient and accurate high-throughput computations of NMR chemical shifts, our approach can accelerate theory-guided design of liquid solutions for various applications.