I have one question - I've done simulations using Car Parrinello Molecular Dynamics together with Thermodynamic Integration method (liquid phase), hence I have calculated Helmholtz energies for my molecule at given values of its torsion angle. My molecule has two global energy minima at radically different values of torsion angle, and I'd like to calculate % abundance of each conformer (population).

Now, if I'd've done statical calculations using Gaussian I'd have Gibbs energies, and for those two conformers I'd calculate it like that:

conformer 1: $\frac{e^{\frac{-\Delta G_{conf_{1}}}{RT}}}{e^{\frac{-\Delta G_{conf_{1}}}{RT}} + e^{\frac{-\Delta G_{conf_{2}}}{RT}}}$ conformer 2: $\frac{e^{\frac{-\Delta G_{conf_{2}}}{RT}}}{e^{\frac{-\Delta G_{conf_{1}}}{RT}} + e^{\frac{-\Delta G_{conf_{2}}}{RT}}}$

But don't know how to do it having $\Delta F$

Can I just use Boltzmann distribution in form like this:

conformer 1: $\frac{e^{\frac{-\Delta F_{conf_{1}}}{kT}}}{e^{\frac{-\Delta F_{conf_{1}}}{kT}} + e^{\frac{-\Delta F_{conf_{2}}}{kT}}}~?$


I'm not sure if this is exactly the point at issue, but here goes anyway.

You have a simulated Helmholtz function $F(\theta)$ tabulated as a function of the torsional angle $\theta$. From this, it is possible to define a normalized probability density, as a function of $\theta$, like this: $$ \mathcal{P}(\theta) = \frac{1}{\mathcal{N}}\exp[-F(\theta)/k_BT] \quad\text{where}\quad \mathcal{N} = \int_0^{2\pi} d\theta \exp[-F(\theta)/k_BT] . $$ This is similar to your suggested formula, but the normalizing factor is an integral, rather than a sum of two possibilities.

Couple of technical points. Firstly, the reference value from which you measure $F(\theta)$ doesn't matter. In other words, you can add or subtract any constant $F_0$ from $F(\theta)$, without affecting the result. One sensible choice is to redefine $F(\theta)\rightarrow F(\theta)-F_0$ where $F_0$ is the minimum value of $F(\theta)$ that you have determined. So now, $F(\theta)\ge 0$ for all $\theta$. Secondly, you'll almost certainly have to calculate the integral, to get $\mathcal{N}$, numerically.

Now, your two low-energy conformers should correspond to two reasonably sharp peaks in $\mathcal{P}(\theta)$. If you integrate $\mathcal{P}(\theta)$ over each peak separately, you will get the contributions of the populations of each conformer to the total. As long as the two peaks are reasonably well separated, and the probability decreases to a very low value in between them, the result should not be over-sensitive to the precise limits of integration. You just need to capture the weight in each peak separately. Assuming that there are no other low-lying conformations, of course, the two probabilities (relative populations) should add up to 1.

You mentioned that a static calculation using the Gaussian package would give you two Gibbs free energies. I don't know exactly what the basis is, for this package to give you $G$ rather than $F$ (or, indeed, how it estimates free energies rather than just energies). But for a torsional degree of freedom, I doubt that it makes much difference whether you consider your system to be at constant volume, or constant pressure, so I would not be inclined to worry about the difference between $G$ and $F$ in this case.

[EDIT following OP comment]

It seems that the Gaussian calculation is for isolated gas-phase molecules, in well-defined conformational minima. The fact that $G$ is calculated in this situation rather than $F$ is not important here.

Integrating $\mathcal{P}(\theta)$ over each peak separately is illustrated below. I've shown $F(\theta)$ and $\mathcal{P}(\theta)$ schematically.

enter image description here

The integrals are over the first peak, giving the blue area $$ p_1 = \int_{\theta_1}^{\theta_2} \mathcal{P}(\theta) d\theta $$ and over the second peak, giving the red area $$ p_2 = \int_{\theta_3}^{\theta_4} \mathcal{P}(\theta) d\theta . $$ The choice of the integration limits $\theta_1\ldots\theta_4$ should not be crucial, provided the probability density is very small in the regions between the peaks. If $\mathcal{P}(\theta)$ is properly normalized, $p_1+p_2=1$, and $p_1$ and $p_2$ are the desired populations.

| cite | improve this answer | |
  • $\begingroup$ Thanks a lot! This is how Gaussian calculates: gaussian.com/thermo $\endgroup$ – LonelyStudent Jan 2 '19 at 6:34
  • $\begingroup$ I now totally see this and understand. Thank you very much it makes sense, totally. One thing I am not sure (because English is not my first language) you said integrate P over each peak, meaning? $\endgroup$ – LonelyStudent Jan 2 '19 at 7:37
  • $\begingroup$ I've edited my answer to clarify. Let me know if things are still not clear. If you wish to accept the answer, please click on the check mark to the left (it is not obligatory; you may wish to wait for more suitable answers). $\endgroup$ – user197851 Jan 2 '19 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.