I have generated at least 10 unique conformers for a peptide that I would like to calculate the "theoretical" proton affinity. Using DFT, each conformer is assigned a ΔH298K value.
By arranging the relative ΔH298K values from lowest to highest I can determine the weighted average (ΔwtH298K) and use this value in an isodesmic reaction scheme.
Now, through some textbooks and discussion with my advisor, the correct term to use to obtain the Boltzmann weights is the relative ΔG298K. I have tried using the weights obtained from relative ΔG298K to my ΔwtH298K value, but get widely different results.
http://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution#Derivation_and_related_distributions
http://en.wikipedia.org/wiki/Conformational_isomerism#Population_distribution_of_conformers
The weighted average is about the microstates (Ei), not macrostates (G). Please find the attached paper from my former diploma supervisor.
It depends on how the conformers were generated. Samples drawn from widely separated time-points from a well equilibrated MD simulation should already be Boltzmann-weighted and don't need another round of weighting, but in that case, 10 conformers is probably too few. On the other hand, if they come from 10 random points in configuration space each taken to the nearest local energy minimum, then in principle, they need re-weighting and ΔG298K is the way to do it because the idea is that each conformer represents an independent (i.e. non-overlapping) region of configuration space and these regions have different populations which are determined by their relative free energies.
But what is it that you should average? It's NOT PA! That's because the conformer populations can change between the protonated and deprotonated states. Instead, you need the enthalpies of each conformer. So you need to take whatever your QM tells you is the enthalpy (typically some sort of enthalpy of formation from elements) of each protonated conformer, do the weighted average of these to get the enthalpy of the protonated state, then do the same for the deprotonated state and subtract them to get PA. But this will need a correction for the fact that the deprotonated-side calculations contain one less proton. So this will depend on what your program thinks is the enthalpy (of formation?) of a proton at 298 K is. I guess you could find this out by asking your program to do a PA calculation for a single conformer of some small molecule and then compare that to the difference in the calculated enthalpies of the protonated and deprotonated species.
Oh, in all this I'm assuming that by PA you mean the enthalpy change, not the free energy change, associated with protonation.
@Donald Bashford and Milán Szőri - If I understand you correctly, the procedure I used to calculate PA is similar to what you described. I generate 10 lowest energy conformers for the neutral and protonated (each at a local minima) and obtain their ΔH298K values. I order them by relative energy( Δ(ΔH298K) ), determine the weights, and derive the weighted average.
However, instead of subtracting the Protonated and Neutral directly, I use an isodesmic proton-transfer reaction w/ Ethylamine (well established PA value):
MH+ + EtNH2 → M + EtNH3+
PAM = PAEtNH2 + ΔHrxn
Using this method I have found that the PAcalc values are nicely comparable to my PAexp values. However if I use ΔG298K and Δ(ΔG298K) orderings and weights, my PAcalc values do not even come close to making any sense.
I have not taken statistical mechanics in my undergraduate course but the response and article that was attached from Milán Szőri seem to make sense. Your comments seem to confirm my method of deriving accurate PAcalc.
- Badges
- Science method