Boltzmann Distribution in Molecular Dynamics Simulation?

Question

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?

Answer 1

Within statistical mechanics (SM) a molecular property $X$ is computed by $$\left<{X}\right>_{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<{X}\right>_{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<{X_{SM}}\right>$ 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<X\right>_{SM} = \lim_{t \to \infty} \left<{X}\right>_{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).

Answer 2

This question is actually very broad as the answer can change quite a bit based on what kind of thing you are doing with your MD simulation.

The first thing to note is that the Boltzmann distribution only applies at thermal equilibrium. So if you start an MD simulation from a non-equilibrium state and then try to calculate the average of some property as you suggest, you will notice that the value behaves very erratically and then eventually converges to the expected value. This convergence basically means the system equilibrated so that the Boltzmann distribution actually is obeyed.

Another thing is that it is not true that a system will be immediately obeying a Boltzmann distribution just because you run an MD simulation. Perhaps you are thinking of the fact that many MD codes can begin a simulation by seeding the velocities of the particles according to a Maxwell-Boltzmann distribution. This is by no means required, but is useful because it generally gets you closer to equilibrium than if you just started with velocities from zero. The reason you can't just start from a Boltzmann distribution is because the Boltzmann distribution assumes a knowledge of the states of the system. Generally these states and their properties are what you are trying to find. All the MB distribution asks for, however, is a temperature and the masses of the particles. These are user-defined so no problem there.

A third note is that this answer somewhat depends on if by MD you mean classical molecular dynamics for quantum molecular dynamics. The Boltzmann distribution is somewhat more useful for quantum systems because quantum systems have states of discrete energy which works very naturally with the Boltzmann distribution. That's not to say that you can't use the Boltzmann distribution for classical, but it is somewhat less natural because in order to get meaningful statistics on the dynamics of the higher-energy state you will have to use a constant energy ensemble. And even then, how constant of energy is needed as too tight a constraint might force very non-physical dynamics?

To be extra clear about the main part of your question, you say:

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).

The misunderstanding is that an MD simulation does not guarantee higher energy states are less likely to be observed for short simulations times. If you start very far from equilibrium (because how are you going to guess the "equilibrium structure" of a liquid for instance) then you might wander around in high energy states for quite some time before thermal equilibrium is found. This is most definitely not a Boltzmann distribution and the states will not be Boltzmann averaged (and won't follow MB statistics). Everything you say after this is fine, but you implicitly assumed that we started from equilibrium when in fact we do not.

If there is anything more specific I can try to clarify, or if I am missing the heart of your question, please comment and I will edit.