Multivalent ions (Mg2+) in RNA tertiary structure folding can be strongly correlated and thus cannot be treated by mean-field theories such as the Poisson-Boltzmann equation. We recently developed a statistical mechanical model (TBI) to account for ion correlation by considering ensemble of discrete ion distributions. Experimental tests show that the TBI model gives improved predictions for nucleic folding folding stability over the Poisson-Boltzmann equation, which generally underestimates the (multivalent) ion-dependent folding stability due to ignoring the ion correlation. Using the TBI theory, we investigate the folding energy landscape for a simple system with loop-tethered short DNA helices and find that Na+ and Mg2+ play contrasting roles in helix–helix assembly. High [Na+] (>0.3 M) causes a reduced helix–helix electrostatic repulsion and a subsequent disordered packing of helices, while Mg2+ of concentration > 1 mM is predicted to induce a more compact and ordered helix–helix packing. Mg2+ is much more efficient in causing nucleic acid compaction and is predicted to induce a collapse transition around 1mM of [Mg2+].