# Why are MD simulations necessary for obtaining Boltzmann populations?

Given that MD simulations converge to the Boltzmann distribution $$\rho \sim \exp(-\beta \epsilon)$$ after sufficiently long times, and all the macroscopic quantities can be computed from the Boltzmann distribution itself, what is the need for MD simulations?

Specifically, I am asking in the context of short peptides (tripeptides or tetrapeptides).

For instance, in Beck, D. A. C.; Alonso, D. O. V.; Inoyama, D.; Daggett, V. Proc. Natl. Acad. Sci. U. S. A. 2008, 105 (34), 12259–12264, the authors use MD to generate Ramachandran distributions of pentapeptide conformations at a constant temperature. This distribution should obey statistical mechanics, and specifically the Boltzmann distribution. So in theory one should be able to write down the distributions using the Boltzmann weights as follows,

$$\rho(\{\phi_{i},\psi_{i}\}) \sim \exp(-\beta V(\{ \phi_{i},\psi_{i}\})).$$

Here, $$\{ \phi_{i},\psi_{i}\}$$ denotes the set of Ramachandran angle coordinates.

Why should I run MD to get the same distributions?

• In the paper you reference, it is noted that the conformations do not converge and are dynamic in range. Most important and more broadly, experimentally-determined structures are static and do not represent full state space of the system. Commented Oct 11, 2021 at 1:58
• I still have a question. Even if they get dynamic pentapeptide when they use MD to get the distributions. They can simply use boltzman weight to express that right? why should I use MD to get this result. Is boltzmann weight incorrect in predicting conformational distributions? @ToddMinehardt Commented Oct 11, 2021 at 2:18

In principle you can, assuming you are given $$V(\{\bf{r_i}\})$$ where $$\{\bf {r_i} \}$$ is the set of variables that define a configuration in the system - this is typically the coordinates of the atoms, but it might include things like external fields. The problem is in practice it will be at best horribly, horribly inefficient, and in fact almost certainly impossible.

To make progress let's try and define your problem at little more carefully. You have written the potential energy of the system as a function of just two variables, $$\{ \phi_{i},\psi_{i}\}$$. This is not the case - it is a function of the coordinates of all the atoms in the system, assuming no external fields, and so to get a function of just these two variables we are going to have to average over all configurations which contain the atoms in positions satisfying the constraint that $$\{ \phi_{i},\psi_{i}\}$$ are the values we desire. Thus in fact we should write

$$\rho(\{\phi_{i},\psi_{i}\}) \sim \int d{\bf {r_1}}d{\bf {r_2}}... d{\bf {r_N}}\exp(-\beta V(\{ \phi_{i},\psi_{i}\};\{\bf{r_i}\}))$$

where the integral has limits such that the atomic configurations considered have the required angles, $$\{\phi_{i},\psi_{i}\}$$.

Now in principle this complicated multi-dimensional integral might be able to be solved analytically. If this is the case just calculate the value for all $$\{\phi_{i},\psi_{i}\}$$ you are interested in, normalise the probabilities appropriately, and you are done. In practice for any realistic forcefield you will never be able to perform this integral analytically. Thus you have to resort to numerical integration on the computer.

But now you have a secondary problem. To do this multidimensional integral on a computer naïvely is going to be impossible for anything but the smallest system. Let's say we just have 35 atoms, thus approximately 100 degrees of freedom, and we want to have just 10 points along each direction. That means we have approximately $$10^{100}$$ configurations to sample for each $$\{\phi_{i},\psi_{i}\}$$. This is impossible even on the most powerful computers we have today. And it will be impossible for a very long time.

But there's another problem, even if we could do this it will likely be extremely inaccurate to the point of uselessness. The problem is the vast, vast, vast majority of the atomic configurations will contribute almost nothing to the integral; their energy will be just too high - remember these are integrals over all space. As an example consider the case where the hydrogen bonded to a $$\alpha$$ carbon is stretched to almost infinity. We know this is very high energy, and therefore through the exponential will make a negligible contribution to the integral, but our above naïve sampling scheme makes no guarantee that we won't include such configurations and so pointlessly calculate their contribution to the probability, which will be zero. But worse is the flip side of this: We have absolutely no guarantee that we will actually include configurations that do have large contributions to the probability! Thus unless we have many more than 10 points along each axis most likely our result will be utterly, utterly wrong.

This is where Molecular Dynamics (or, alternatively, Monte Carlo) comes in. It is an "intelligent" way of sampling the above integral such that we only consider atomic configurations that have a sizeable contribution to the configuration integral. Thus we can calculate accurate estimates to the quantities of interest by using a relatively small sample of important configurations, and thus making the evaluation of those properties via the configuration integral viable on a computer in a reasonable amount of time and memory.

Finally note this is not specific to the case here - the same argument applies to virtually any quantity of interest, and the conclusion remains the same: In statistical physics the prime uses of Molecular Dynamics and Monte Carlo are as intelligent sampling techniques to overcome the Curse of Dimensionality