# Determining the torsion energy profile of butadiene by molecular dynamics

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?

• I would change the tag thermodynamics for statistical mechanics. – B. Kelly May 6 '20 at 7:09
• @CharlieCrown you can ping editors directly (i.e. @orthocresol), so that I get notified ;) – orthocresol May 6 '20 at 10:22
• @orthocresol yes but sometimes when you play with the bull, you get the horns – B. Kelly May 6 '20 at 17:20

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 $$$$P_i = \frac{\exp ^{-E_i/RT}}{Z}$$$$

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

$$$$n_i N = \exp ^{-E_i/RT}$$$$

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

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

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

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

• You are right, I made a leap at that part, thanks for pointing it out. However, it does not appear to cause the mismatch. Since N and Z and in the logarithm, they will subtract/add a constant to the energies. If I calculate the relative energies afterwards, the result is the same. – Raphaël May 6 '20 at 11:11
• @Raphaël The probability of one conformer to another would not need the partition function, but that is a special case. Overall, just by making things relative, you still have a Z, just a simpler Z made up of a sum of $\Delta E$'s rather than a sum of $E_i$'s. i.e, if $n_i = N \frac{\exp^{-E_i/RT}}{Z}$ then you could say $\frac{n_i}{n_j}=\exp^{-\Delta E_{\rm ij}/RT}$, but that cancellation of both N and Z is because it is a ratio of 2 specific configurations. – B. Kelly May 6 '20 at 17:15