Determining the torsion energy profile of butadiene by molecular dynamics

Question

Disclaimer: this question is related to my previous question, although it is much more detailed

Consider the following rotation of butadiene:

I want to determine the relative energy of each dihedral angle by molecular dynamics (MD).

My initial idea was to simulate butadiene in the NVT ensemble and use the relative frequencies of each dihedral angle to determine the relative energy. In the NVT ensemble, the probability of each microstate $P_i$ is:

$$P_i = \frac{\mathrm e^{-E/RT}}{Z}$$

Thus, I performed MD simulations in the NVT ensemble and recorded the dihedral angle periodically. With $N$ total dihedral angles, I made a histogram of the dihedral angles (e.g. a bin every 5 degrees). If we define the number of angles contained in a given bin as $n_i$, I am using the following equation:

$$n_i = \frac{\mathrm e^{-E/RT}}{N}$$

This equation can be rearranged to isolate the approximate energy associated with that dihedral angle:

$$-RT \ln(n_i N)= E$$

However, the curve obtained does not match the distribution I obtained by a simple relaxed scan of the dihedral:

Where is the problem in my method?

Answer 1

I really only want to make a comment, but it is too long...

Your equation $n_i = \exp^{-E_i/RT}/N$ does not make sense to me.

I agree with \begin{equation} P_i = \frac{\exp ^{-E_i/RT}}{Z} \end{equation}

where $Z = \sum \exp ^{-E_i/RT}$ is the configurational partition function.

I can also accept that using a histogram, you could say $P_i = n_i/N$ where $n_i$ is number of samples in that bin, and N is total samples.

I do not see however how

\begin{equation} n_i N = \exp ^{-E_i/RT} \end{equation}

This does not make sense to me. Perhaps you could swap $P_i$ for $n_i/N$ so that you have \begin{equation} P_i =\frac{n_i}{N} = \frac{\exp ^{-E_i/RT}}{Z} \end{equation}

which could lead to \begin{equation} n_i = \frac{N}{Z}\exp ^{-E_i/RT} \end{equation}

Which could be rearranged to get \begin{equation} E_i =-RT \ln \left(\frac{n_i }{N}Z\right) \end{equation}

If you are comparing to a reference energy then you could get

\begin{eqnarray} E_i - E_{\rm ref} = \Delta E_{\rm i,ref} &=&-RT \ln \left(\frac{n_i }{N}Z\right) +RT \ln \left(\frac{n_{\rm ref} }{N}Z\right) \\ &=& -RT\ln \left(\frac{n_i }{n_{\rm ref}}\right) \end{eqnarray}

It would be interesting to see if this works? I don't have the raw data, so I cannot try it (I also don't have the time, only the interest).