Boltzmann Distribution in Molecular Dynamics Simulation? The Boltzmann distribution is often used to describe the population of energy states at different temperatures, i.e. something like

$$\dfrac{g_i}{g_j} = \mathrm{exp}\left(-(E_i-E_j)\times \frac{1}{k_B\cdot T}\right).$$

For two given states $i,j$ with energies $E_i,E_j$, we thus have the population given by $g_i$ and $g_j.$

So if I have some observable $\sigma$ calculated for each state $i$, the natural way to calculate the observable of the ensemble would be

$$ \sigma_\mathrm{ensemble} = \left(\sum_i g_i \times \sigma_i\right) / ~\sum_i g_i~.$$

So, if I calculate $n$ ground states, let's say by geometry optimization of several guess structures, I employ the above weighting scheme to calculate the observable for the ensemble.

But here is my problem: If I perform an MD simulation instead, then the states should already be automatically populated according to the Boltzmann population (e.g. higher energy states are less likely to be observed in the MD simulation). So the population of the MD simulation already converges to the Boltzmann population for increasing simulation times $t$

$$ \lim_{t\rightarrow\infty} \frac{g_i^{MD}}{g_j^{MD}} = \frac{g_i}{g_j} $$ 

Now I am confused, because that would basically mean that there is a continuous transition between e.g. two states that clearly have to be boltzmann averaged up to the "long MD simulation", that is automatically populated accordingly. Where in the latter case, I would accidentally apply the Boltzmann population weighting twice.

Are there any resources on this weird population dichotomy or am I just overlooking something?
Within statistical mechanics (SM) a molecular property $X$ is computed by $$\left&lt;{X}\right&gt;_{SM}=\sum^{states}_i X_i p_i$$where $X_i$ is the value of $X$ for energy state $i$ and $p_i$ is the probability of being in energy state $i$ with energy $E_i$:$$p_i=\frac{e^{-E_i/kT}}{\sum_i e^{-E_i/kT}}$$Within molecular dynamics (MD) the corresponding property is computed by$$\left&lt;{X}\right&gt;_{MD}=\frac{1}{M}\sum^M_i X(t_i)$$where $M$ is the number of time-steps and $X(t_i)$ is the value of $X$ at time $t_i$.

You can also compute $\left&lt;{X_{SM}}\right&gt;$ from an MD by finding all the "unique" states (e.g. conformations) and computing their energy.  "Unique" means that even if a state is observed more than once in the MD you only include it once in the sum over states. 

In principle $$\left&lt;X\right&gt;_{SM} = \lim_{t \to \infty} \left&lt;{X}\right&gt;_{MD}$$ but this assumes that the method used to calculate the energy and forces are identical in both calculations.

For example, for MD explicit solvent calculations with periodic boundary conditions it not really possible to calculate the energy of a state ($E_i$).  Also, vibrational effects must be included in $E_i$, which is typically done using the harmonic oscillator approximation, while the MD will include anharmonic effects. Similarly, the harmonic oscillator approximation is quantized (QM) while the vibrations in the MD simulation are not (classical).