Thermodynamics in Solution (solvation) from quantum chemistry

Question

The gist of the question is "how do I get a heat of formation in solution" - but let me give it some context to be more specific.

I am interested in calculating the thermodynamic values of molecules in solution. Due to the number of jobs, time etc. I have to use DFT methods. These are not as accurate as for example CCSD related methods but still provide good answers in reasonable timeframes (plus I'm checking up on estimated values, so I don't immediately need the highest acutacy).

As atomisations are hard to calculate with DFT methods (see a book by Cramer & comment on the ORCA forum), isodesmic reactions offer a better option. So far so good. This allows me to use literature values for reference compounds to calculate the heat of formation of a species.

Now I need the heat of formation of the species in the liquid phase. I can repeat my quantum chemistry calculations with a solvation model and again obtain values - however I now lack the reference compounds. Based on my understanding of chemistry (see textbook by Cramer) I can relate the Gibbs free energy of the gas phase molecule with the solvated molecule. So using the following:

• $G_{liq} = G_{gas} + G_{sol}$

• $\delta G = \delta H_f - \delta S ~T_{LO}$

  I can work out the liquid phase heat of formation of the molecule. (Please correct me if I am wrong.)

Now here is a little problem: For organic molecules, there can be multiple local minima - but only one absolute total minimum.

The total minimum can correspond to a different structure in the gas and liquid phase, so do I compare structures and consider them or do I use the total minimum in both cases, I.e the lowest energy structure to calculate my thermodynamic data?

Reference: Essentials of Computational Chemistry: Theories and Models, 2nd Edition Christopher J. Cramer http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0470091827.html

Answer 1

The best way is to take all the structures into account and use Maxwell-Boltzmann distribution to weigh them accordingly. Here is what I would do
For each chemical structures/isomers (here I used F to denote Gibbs free energy), $ F_{liq} = F_{gas} + F_{solvation}$
Now, $ F_{gas} = E_{scf} - k_BTlnQ$
Where, $Q=q_{translational} \times q_{vibrational} \times q_{rotational} $
Using appropriate method (implicit solvation or QM/MM) you can find $F_{solvation}$. Now you have $F_{liq}$ for each isomers. Then use the following method to weigh them according to Maxwell-Boltzmann distribution.

In the above image, I use $\epsilon_i$ to denote $F_i$.