The harmonic oscillator model and the pure WS result in ⟨r⟩ h = r e because the harmonic potential is symmetric. However, MD simulations with initial conditions from WS (WS+MD) produce an average value ⟨r⟩ WS+MD very close to that computed with the numerical ground state anharmonic wavefunction ⟨r⟩ anh . In the case of the diatomic molecule, the SWS distribution can be exactly mapped to WS. Therefore we can get the same solution ⟨r⟩ SWS+MD = ⟨r⟩ WS+MD by scanning the parameter τ . The resulting curve has a single minimum at τ = 2 fs where ⟨r⟩ SWS+MD = ⟨r⟩ WS+MD = ⟨r⟩ anh . Figure also shows the results from MBS+MD simulations at various temperatures, where the temperature T is mapped to τ through the Debye temperature relation τ = ν -1 = h/(k B T ). The resulting value coincides with the expected one at T around 3000 K (τ = 16 fs). The trend for MBS shows only the growth of ⟨r⟩ with the increase in temperature. Thus MBS requires an a priori knowledge of the appropriate temperature, in contrast to SWS+MD, where the best τ can be found in the variational procedure by minimizing the ⟨r⟩. The results for the amplitude l show that the harmonic value is close to those from anharmonic and WS+MD calculations. This is expected since the second moment of the internuclear distribution is dominated by the harmonic potential. The single minimum of ⟨r⟩ and l in the SWS+MD simulations suggest a similar behavior for other observables, which can be used for finding optimal parameter τ for larger molecules. The obvious choice is an energy-related parameter, for example, the average total energy (electronic and nuclear) ⟨E tot ⟩. The other option is the kinetic energy of nuclei (⟨E kin ⟩), which is often expressed in terms of temperature through the relation