The conformationally averaged free energy

The free energy averaged over $N_\mathrm{conf}$ conformations of molecule X is

$G^\circ=-RT \ln \left( \sum_I^{N_\mathrm{conf}} e^{-G^\circ_I/RT} \right)$  (1)

However, the way to compute an average property of the conformations $\left\langle x\right\rangle$ is

$\left\langle x\right\rangle=\sum_I^{N_\mathrm{conf}} p_Ix_I$  (2)

where

$p_I=\frac{e^{-G_I^\circ/RT}}{\sum_J^{N_\mathrm{conf}} e^{-G_J^\circ/RT}}$ (3)

But this equation cannot be used to compute the conformationally averaged free energy because

$G^\circ \ne \left\langle G^\circ \right\rangle$

Why not? First, let's see where the first equation comes from.

Obtaining equation (1)

The first equation comes from the statistical mechanical definition of the Helmholtz free energy found any any p-chem textbook

$A^\circ = -RT \ln  \left( \sum_i^{\mathrm{states}} e^{-\varepsilon_i/kT} \right)=-RT \ln  \left( q \right)$

The sum over microstates can be split into sums over microstates of each conformation I

$A^\circ = -RT \ln  \left( \sum_I^{N\mathrm{conf}} \sum_{i\in I}^{\mathrm{states}} e^{-\varepsilon_i/kT} \right)$

$= -RT \ln  \left( \sum_I^{N\mathrm{conf}} q_I \right) = -RT \ln \left( \sum_I^{N_\mathrm{conf}} e^{-A_I^\circ/RT} \right)$

To get the corresponding expression for the Gibbs free energy [Eq (1)]:

$G^\circ = A^\circ +p^\circ V $

$=-RT \ln \left( \sum_I^{N_\mathrm{conf}} e^{-A^\circ_I/RT} \right)-RT \ln \left(  e^{-p^\circ V/RT} \right)$

$=-RT \ln \left( \sum_I^{N_\mathrm{conf}} e^{-G_I^\circ/RT} \right)$

The difference between equation (1) and (3) is the conformational entropy

$G^\circ - \left\langle G^\circ \right\rangle=-RT \ln \left( \sum_J^{N_\mathrm{conf}} e^{-G_J^\circ/RT} \right)\sum_I^{N_\mathrm{conf}} p_I$

$+RT\sum_I^{N_\mathrm{conf}} p_I\frac{-G_I}{RT}$

$=RT \sum_I^{N_\mathrm{conf}} p_I \left( \frac{-G_I}{RT}- \ln \left( \sum_J^{N_\mathrm{conf}} e^{-G_J^\circ/RT} \right) \right)$

$=RT \sum_I^{N_\mathrm{conf}} p_I \left( \ln \left(  e^{-G_I^\circ/RT} \right) - \ln \left( \sum_J^{N_\mathrm{conf}} e^{-G_J^\circ/RT} \right) \right)  $

$=RT \sum_I^{N_\mathrm{conf}} p_I \ln(p_I) =-TS_{\mathrm{conf}} $

Predicting binding free energies with electronic structure theory: thermodynamic considerations - Part 3

2015.01.25:  I have summarized discussion of this and related issues in this paper.

Predicting absolute binding free energies of biologically relevant molecules (in aqueous solution at pH 7 using electronic structure theory without empirical adjustments to within 1 kcal/mol) is one of the holy grails of computational chemistry. Recent work by Grimme and co-workers (here and here) have come close with mean absolute deviations of about 2 kcal/mol for host-guest complexes. This and other work has shed some light on why it is so difficult to predict binding free energies in aqueous solution, which I will discuss here and in two previous posts (Part 1 and Part 2). In addition I'll talk about further improvements that could potentially increase the accuracy further.

The approach is based on the following equation

$\Delta G^\circ = \Delta E_\mathrm{gas} + \Delta G^\circ_{\mathrm{gas,RRHO}}+\Delta \delta G^\circ_{\mathrm{solv}}$

where $\Delta E_\mathrm{gas}$ is the change in electronic energy and $\Delta G^\circ_{\mathrm{gas,RRHO}}$ the change in the translational, rotational, and vibrational free energy computed using the rigid-rotor and harmonic oscillator approximation, respectively - both evaluated in the gas phase. The vibrational frequencies are computed at a lower level of theory compared to that used to compute $\Delta E_{\mathrm{gas}}$. $\Delta \delta G^\circ_{\mathrm{solv}}$ is the change in solvation free energy computed using a continuum model of the solvent, such as COSMO.

Solvation Thermodynamics

Background. Most continuum models (CMs) of solvation compute the solvation free energy as the difference between the free energy in solution and the gas phase electronic energy

$\delta G^\circ\mathrm{_{solv}(X)}=G\mathrm{^{\circ,CM}_{soln}(X)}-E\mathrm{_{gas}(X)}$

$G\mathrm{^{CM}_{soln}(X)}$ typically contain energy terms describing the electrostatic interaction of the molecule and the continuum as well as the van der Waal inteactions with the solvent and free energy required to create the molecular cavity in the solvent (cavitation). The electrostatic interaction with the solvent alters the molecular wavefunction and is computed self-consistently.

Some software packages automatically compute $G\mathrm{^{\circ,CM}_{soln}(X)}$ and $E\mathrm{_{gas}(X)}$ in one run, while others only compute $G\mathrm{^{\circ,CM}_{soln}(X)}$. Also, some programs just compute the electrostatic component of $G\mathrm{^{\circ,CM}_{soln}(X)}$ by default.  However, the van der Waals and, especially, the cavitation component can make sizable contributions to the binding free energy and must be included for accurate results.  It is worth noting that any hydrophobic contribution to binding will derive primarily from the change in cavitation energy.

$G\mathrm{^{CM}_{soln}(X)}$ contain parameters (e.g. atomic radii) that are adjusted to reproduce experimentally measured solvation free energies

$\delta G\mathrm{^\circ_{solv}(X)}=G\mathrm{^{\circ,exp}_{soln}(X)}-G\mathrm{^{\circ,exp}_{gas}(X)}$

The standard state for both $G\mathrm{^{\circ,exp}_{soln}(X)}$ and $G\mathrm{^{\circ,exp}_{gas}(X)}$ is generally 1 M.  The latter is the reason why a 1 M reference state also must be used when computing $\Delta G^\circ_{\mathrm{RRHO}}$

Atomic radii. The solvation energy is computed using a set of atomic radii that define the solute-solvent boundary surface. These radii are usually obtained by fitting to experimentally measured solvation energies. Accurate solvation energies should not be expected from methods that use iso-electron density surfaces or van der Waals radii without additional empirical fitting. When using fitted radii one should use the same level of theory for the solute as was used in the parameterization.

Ions. For neutral molecules solvation free energies can be measured with an accuracy of roughly 0.2 kcal/mol and reproduced theoretically to within roughly 0.5 kcal/mol.  However, the solvation energies of ions cannot be directly measured and must be indirectly inferred relative to a standard (typically the solvation energy of the proton). The experimentally obtained solvation energies are typically accurate to within 3 kcal/mol and can be reproduced computationally with roughly the same accuracy.  The solvation energy of ions are therefore an especially likely source of error in binding free energies - especially if the ionic regions of the molecules become significantly desolvated due to binding.

Gas phase vs solution optimization.  The fitting of the radii described above is usually done using gas phase optimized structures only, i.e. any change in structure and corresponding rotational and vibrational effects are "included" in the radii via the parameterization.  However, for ionic species gas phase optimization can lead to significantly distorted structures or even proton transfer and in these cases solution phase optimizations and, hence, vibrational frequency calculations, tend to be used. However, numerical instability in the continuum models can make it necessary to increase (i.e. make less stringent) the geometry convergence criteria and can lead to more imaginary frequencies than in the gas phase. One option is to compute the vibrational contribution to $\Delta G^\circ_{\mathrm{RRHO}}$ using gas phase optimized structures.

When using solution phase geometries the gas phase single point energies needed to evaluate $\delta G\mathrm{^\circ_{solv}(X)}$ represent added computational expense and it can be tempting to use solution phase free energies to evaluate the binding free energies

$\Delta G^\circ = \Delta G\mathrm{^{\circ,CM}_{soln}}- \Delta G\mathrm{^{\circ,CM,RRHO}_{soln}}$

One problem with this approach is that $\Delta G\mathrm{^{\circ,CM}_{soln}}$, unlike $\Delta E_{\mathrm{gas}}$, is not systematically improveable due to the empirical parameterization.

Cavities. The atomic radii and corresponding cavity generation algorithms are parameterized for small molecules. For more complex molecules such as receptors this can lead to continuum solvation of regions of molecules, e.g. deep in the binding pocket, that are not accessible to the molecular solvent. Furthermore, any solvent molecule inside such pocket is likely to be quite "un-bulk-like" and not well-represented by the bulk solvent or fixed by the underlying parametrization.  However, how big an error this may introduce to the binding free energy is not really known, but certain models for the cavitation energy have been shown to give unrealistically large contributions to the binding free energy.

Explicit water molecules. Adding explicit solvent molecules to the receptor and/or ligand can potentially lead to more accurate results. For example, including explicit water molecules around ionic sites reduces the strong dependence of the solvation energy on the corresponding atomic radii. Also, "un-bulk-like" water molecules now are treated more naturally and the risk of solvating non-solvent-acesssible regions is reduced somewhat.  However, adding explicit solvent molecules increases the computational cost by increasing the CPU time needed to compute energies, perform conformational searches, and compute vibrational frequencies.

There are several approaches to include the effect of explicit solvent molecules in the binding free energy.  Bryantsev and co-workers suggest computing the solvation energy by

$\delta G^\circ_{\mathrm{solv},n}\mathrm{(X)}=\Delta G^\circ_{\mathrm{gas}}\mathrm{(X(H_2O)_n})+\delta G^\circ_{\mathrm{solv}}\mathrm{(X(H_2O)_n)}-$

$\delta G^{\circ,liq}_{\mathrm{solv}}\mathrm{((H_2O)_n)}$

where

$\Delta G^\circ_{\mathrm{gas}}\mathrm{(X(H_2O)_n})=G^\circ_{\mathrm{gas}}\mathrm{(X(H_2O)_n})-G^\circ_{\mathrm{gas}}\mathrm{(X})-G^\circ_{\mathrm{gas}}\mathrm{((H_2O)_n})$

and

$\delta G^{\circ,liq}_{\mathrm{solv}}\mathrm{((H_2O)_n)}=\delta G^\circ_{\mathrm{solv}}\mathrm{((H_2O)_n)}+RT\ln(\mathrm{[H_2O]/n})$

with "$\circ$" and "$liq$" referring to a standard state of 1 M and 55.34, respectively.  The term $RT\ln(\mathrm{[H_2O]/n})$ is the free energy required to change the standard state of (H$_2$O)$_n$ from 1 M to 55.34/n M.

Bryantsev et al. have shown that using this water cluster approach leads to a smooth convergence of the solvation free energy with respect to the cluster size $n$.  The optimum choice of $n$ is this one where an additional water does changes the solvation energy by less than a certain user defined amount.  One can thereby compute the optimum number of water molecules for the receptor ($n$), ligand ($m$) and receptor-ligand complex ($l$) and then compute the binding free energy as

$\Delta \delta G^\circ_{\mathrm{solv},x}=\delta G^\circ_{\mathrm{solv},l}\mathrm{(RL)}-\delta G^\circ_{\mathrm{solv},n}\mathrm{(L)}-\delta G^\circ_{\mathrm{solv},m}\mathrm{(R)}$

and computing $\Delta E_{\mathrm{gas}}$ and $\Delta G^\circ_{\mathrm{RRHO}}$ as before.  One can show that this corresponds to the free energy change for this reaction

$\mathrm{L(H_2O)_n(aq)+R(H_2O)_m(aq)+(H_2O)_l(liq) \rightleftharpoons RL(H_2O)_l(aq)}+$
$\mathrm{(H_2O)_n(liq)+L(H_2O)_m(liq)}$ (1)

In principle the free energy change is zero for

$\mathrm{(H_2O)_l(liq) \rightleftharpoons+(H_2O)_n(liq)+L(H_2O)_m(liq)+sgn(d)(H_2O)_{|d|}(liq)}$

where $d=l-m-n$ and sgn(d) returns the sign of d. So the free energy change for Reaction (1) can also be computed as the free energy change for

$\mathrm{L(H_2O)_n(aq)+R(H_2O)_m(aq) \rightleftharpoons RL(H_2O)_l(aq)+sgn(d)(H_2O)_{|d|}}$ (2)

However, this is only approximately true in practice due to errors in the computed gas phase and solvation free energies. Furthermore, Reaction 2 does not really lead to any significant reduction in CPU time because the water cluster free energies only have to be computed once. However, if Reaction (2) is used then one must add an additional term correcting for the indistinguishability of water molecules

$G^\circ_{\mathrm{RRHO}}\mathrm{(X(H_2O)_n)} \rightarrow G^\circ_{\mathrm{RRHO}}\mathrm{(X(H_2O)_n)}-RT\ln (n!)$

and similarly for the water clusters.  Using Reaction (1) leads to a cancellation of this term and also maximizes error cancellation in the other energy terms. Similar considerations apply to when using individual water molecules to the balance the reaction instead of water clusters

$\mathrm{(H_2O)_l(liq) \rightleftharpoons (H_2O)_n(liq)+L(H_2O)_m(liq)+sgn(d)|d|H_2O}$

When using many explicit water molecules the error in the continuum solvation energies can be reduced by ensuring that the continuum solvation energy of a single water molecule matches the experimental value of -6.32 kcal/mol at 298.15K as close as possible.

This work is licensed under a Creative Commons Attribution 4.0

Predicting binding free energies with electronic structure theory: thermodynamic considerations - Part 2

2015.01.25:  I have summarized discussion of this and related issues in this paper.  I have also made some changes in the post to reflect what I have learned since I wrote it.

Predicting absolute binding free energies of biologically relevant molecules (in aqueous solution at pH 7 using electronic structure theory without empirical adjustments to within 1 kcal/mol) is one of the holy grails of computational chemistry. Recent work by Grimme and co-workers (here and here) have come close with mean absolute deviations of about 2 kcal/mol for host-guest complexes. This and other work has shed some light on why it is so difficult to predict binding free energies in aqueous solution, which I will discuss here and in a previous post. In addition I'll talk about further improvements that could potentially increase the accuracy further.

The general approach is based on the following equation

$\Delta G^\circ = \Delta E + \Delta G^\circ_{\mathrm{RRHO}}+\Delta \delta G_{\mathrm{solv}}$

where $\Delta E$ is the change in electronic energy and $\Delta G^\circ_{\mathrm{RRHO}}$ the change in the translational, rotational, and vibrational free energy computed using the rigid-rotor and harmonic oscillator approximation, respectively. The vibrational frequencies are computed at a lower level of theory compared to that used to compute $\Delta E$. $\Delta \delta G_{\mathrm{solv}}$ is the change in solvation free energy computed using a polarizable continuum model such as COSMO.

pH and molecular charge. Virtually all binding measurements in aqueous solution are performed in buffer with a constant pH and many ligands and or receptors contain one or more ionizable groups. The charge of an ionizable (acid/base) group in aqueous solution depends on its pK$_a$ and the pH:

$q= \frac{1}{1+10^{\mathrm{pH}-\mathrm{p}K_a}}-\delta$

where $\delta$ is 1 for an acid and 0 for a base.  In most cases, selecting the wrong charge for the ligand and/or host will result in a significant error in the computed binding free energy.  The pK$_a$ can be computed using electronic structure theory or empirically using software such Marvin. However, if the pK$_a$ value is perturbed by the binding the situation may be complicated further. Here I illustrate this point for a simple example where the ligand (L) has a basic group that is neutral when deprotonated and the receptor (R) is non-ionizable.

$\mathrm{R + L(H^+) \rightleftharpoons RL(H^+)}$

The (apparent) experimental binding constant is then

$K^\prime=\mathrm{\frac{[RL]+[RLH^+]}{[R]([L]+[LH^+])}}$

and the corresponding binding free energy is

$\Delta G^{\prime\circ}=\Delta G^{\circ}(+)-RT\ln \left( \frac{1+10^{\mathrm{pH}-\mathrm{p}K_a^c}}{1+10^{\mathrm{pH}-\mathrm{p}K_a^f}} \right)$   (1)

where $\Delta G^{\circ}(+)$ is the binding free energy computed using the charged (protonated) form of the ligand and pK$_a^c$ and pK$_a^f$ are the pK$_a$ values the ligand bound to the receptor and the free ligand, respectively.

If the pK$_a$ is unaffected by the binding, the binding free energy is unaffected by pH and the chosen protonation state of the ligand. For the remaining scenarios it is instructive to plug in some numbers. For example, for pH = 7, pK$_a^f$ = 9 and pK$_a^c$ = 11, the ligand is protonated before and after binding and $\Delta G^{\prime\circ}=\Delta G^{\circ}(+)$ is a good approximation. However, if pH = 7, pK$_a^f$ = 5 and pK$_a^c$ = 3, the ligand is neutral before and after binding and assuming it is charged leads to a 2.7 kcal/mol error in the binding free energy (unless corrected by this equation).  Finally, if pH = 7, pK$_a^f$ = 8 and pK$_a^c$ = 6, the ligand is (91%) charged before and (91%) neutral after binding and assuming it remains charged leads to a 1.4 kcal/mol error in the binding free energy.

For many ligands of interest the pK$_a^f$ can be estimated fairly accurately in a matter of second using programs such as Marvin. The effect of binding on pK$_a^f$ can often be estimated by chemical intuition since hydrogen bonds to charged acid and basic groups tend to, respectively, lower or raise the pK$_a$ even further.  For example, if an amine with pK$_a^f$ = 9 binds to the receptor via hydrogen bonding, then pK$_a^c$ is likely higher than 9 and $\Delta G^{\prime\circ}=\Delta G^{\circ}(+)$ is a good approximation.  However, if pK$_a^f$ is close to 7 then pK$_a^c$ should be computed.  Also, it is possible for charged ligands to change to their neutral state if they bind to hydrophobic or similarly charged receptors.

If pK$_a^f$ is known with some degree of confidence then pK$_a^c$ can be estimated by

$\mathrm{p}K_a^c=\mathrm{p}K_a^f-\Delta G^\circ /RT\ln(10)$

where $\Delta A^\circ$ is the free energy change for

$\mathrm{RLH^+ + L \rightleftharpoons  RL + LH^+}$

If there are several ionizable groups then Eq (1) generalizes to

$\Delta G^{\prime\circ}=\Delta G^{\circ}(-/+)-RT\ln \left(\sum_i \frac{1+10^{n_i(\mathrm{pH}-\mathrm{p}K_{a,i}^c)}}{1+10^{n_i(\mathrm{pH}-\mathrm{p}K_{a,i}^f)}} \right)$

where $\Delta G^{\circ}(-/+)$ is the binding free energy when all acids and bases are deprotonated and protonated, respectively, the sum runs over all ionizable groups and $n_i$ is 1 and -1 if $i$ is a base or acid, respectively.

However, this assumes that the ionizable groups titrate independently of one another, i.e. that the pK$_a$ value of one group is independent of the protonation states of all other ionizable groups.  If that is not the case then it is difficult to give a general expression for the pH-dependent free energy correction in terms of pK$_a$ values (though it can easily be derived for a specific case).

Instead a general expression can be written in terms of Legendre transformed free energies as suggested by Alberty (modified here to electronic structure calculations):

$G^{\prime\circ}(\overline{X})=-RT\ln \left( \sum_i \exp{(-G^{\prime\circ}(X_i)/RT)} \right)$

Here the sum runs over all possible protonation states and

$G^{\prime\circ}(X_i)=G^{\circ}(X_i)-n(\mathrm{H^+})[\delta G^\circ (\mathrm{H^+}) -RT\ln(10)\mathrm{pH}]$

where $G^{\circ}(X_i)$ is the usual standard free energy of protonation state $i$, $n(\mathrm{H^+})$ is the number of ionizable proton in pronation state $i$, and $\delta G^\circ (\mathrm{H^+})$ is the solvation free energy of the proton.  So in the case of ligand L considered above, $n(\mathrm{H^+})$ is 0 and 1 for L and LH$^+$, respectively. $\delta G^\circ (\mathrm{H^+})$ is usually taken from the literature where estimates vary between -265.8 and -268.6 kcal/mol.

Thus, Eq (1) can be rewritten as

$\Delta G^{\prime\circ}=G^{\prime\circ}(\overline{RL})-G^{\prime\circ}(\overline{L})-G^{\circ}(R)$

Since the electronic energy contribution to the standard free energy can be very large in magnitude this form is more easily evaluated

$G^{\prime\circ}(\overline{X})=G^{\prime\circ}({X_0})-RT\ln \left( 1+\sum_{i\ne0} \exp{(-G^{\prime\circ}(X_i)/RT)} \right)$

where $X_0$ is some arbitrarily chosen reference protonation state, for example that for which $n(\mathrm{H^+})$ = 0.  The sum can be combined with that over different conformations, discussed in a previous post.

Other ions
The buffers that are commonly used to regulate the pH also contain other ions, such as Na$^+$, Mg$^{2+}$, Cl$^-$.  At high ion concentrations, it is possible that these ions bind at certain sites in the ligand, receptor, or ligand-receptor complex with sufficient probability that they must be included in the thermodynamics. If so the exact same equations and considerations outlined above for H$^+$ also apply to, e.g. Cl$^-$ and pCl$^-$ (computed from the specified buffer concentration) is used instead of pH.

This work is licensed under a Creative Commons Attribution 4.0

Predicting binding free energies with electronic structure theory: thermodynamic considerations - Part 1

2015.01.25:  I have summarized discussion of this and related issues in this paper.  I have also made some changes in the post to reflect what I have learned since I wrote it.

Predicting absolute binding free energies of biologically relevant molecules (in aqueous solution at pH 7 using electronic structure theory without empirical adjustments to within 1 kcal/mol) is one of the holy grails of computational chemistry. Recent work by Grimme and co-workers (here and here) have come close with mean absolute deviations of about 2 kcal/mol for host-guest complexes. This and other work has shed some light on why it is so difficult to predict binding free energies in aqueous solution, which I will discuss here and in subsequent posts. In addition I'll talk about further improvements that could potentially increase the accuracy further.

The approach is based on the following equation

$\Delta G^\circ = \Delta E + \Delta G^\circ_{\mathrm{RRHO}}+\Delta \delta G_{\mathrm{solv}}$

where $\Delta E$ is the change in electronic energy and $\Delta G^\circ_{\mathrm{RRHO}}$ the change in the translational, rotational, and vibrational free energy computed using the rigid-rotor and harmonic oscillator approximation, respectively. The vibrational frequencies are computed at a lower level of theory compared to that used to compute $\Delta E$. $\Delta \delta G_{\mathrm{solv}}$ is the change in solvation free energy computed using a polarizable continuum model such as COSMO.

Electronic energy

Grimme has shown that dispersion typically makes a very big (>10 kcal/mol) contribution to binding free energies of host-guest complexes. Dispersion corrections are therefore a must if DFT is used to compute the electronic binding energy. Furthermore, Grimme has shown that three-body dispersion makes a non-negligible (2-3 kcal/mol) contribution to the electronic binding energy.  For convergent methods this effect is only included in rather expensive methods that involve triple-excitations such as MP4 and CCSD(T).

Molecular Thermodynamics

The translational, rotational and vibrational thermodynamic contribution to the binding free energy is very large (>10 kcal/mol) and must be included for accurate results.  Some years ago there was a bit of confusion in the literature about whether the RRHO approach was appropriate for condenses phase systems, but Zhou and Gilson have clarified this beautifully.

Standard state.  Most electronic structure codes compute the RRHO energy corrections for an ideal gas, where the standard state is a pressure of 1 bar.  The standard state for solution is 1 M, so the free energy of a molecule X must be corrected accordingly

$G^\circ_{\mathrm{RRHO}}(\mathrm{X}) = G^{\circ,\mathrm{1 bar}}_{\mathrm{RRHO}}(\mathrm{X})-RT\ln (V^{-1})$

where $V$ is the volume of an ideal gas at 1 bar and temperature $T$. At 298K $-RT\ln (V^{-1})$ = 1.90 kcal/mol.  This correction is already included in the solvation energy $\delta G_{\mathrm{solv}}$

Symmetry. Many host molecules and some guest molecules are symmetric, which affects the rotational entropy through the symmetry number ($\sigma$) which is a function of the point group.

$S_{rot}=R\ln\left(\frac{8\pi^2}{\sigma}\left(\frac{2\pi ekT}{h^2}\right)^{3/2}\sqrt{I_1I_2I_3}\right)$

It can be very difficult to build large molecules with the correct point group and most studies use $C_1$ symmetry.  In this case the effect of symmetry must be added manually to the free energy

$G^\circ_{\mathrm{RRHO}}(\mathrm{X}) \rightarrow G^\circ_{\mathrm{RRHO}}(\mathrm{X})+RT \ln(\sigma_X)$

As an example, the popular host molecule corcubit[7] has $D_{7h}$ symmetry and a corresponding $\sigma$ value of 14, in which case the correction contributes 1.56 kcal/mol to the free energy at 298K.

Anharmonicity and low frequency modes. Host-guest complexes can exhibit very low frequency vibrations on the order of 50 cm$^{-1}$ or less, which tend to dominate the vibrational entropy contribution.  Grimme (and many others) have questioned whether the harmonic approximation is valid for such low frequency modes and this is an open research question.  The main problem is that it is very difficult to compute the vibrational entropy exactly.  Most methods for computing anharmonic vibrational are developed to obtain the 2 or 3 lowest energy states, but for very low frequency modes 10-20 states likely significantly populated at room temperature and therefore contribute to the entropy.

In the absence of theoretical benchmarks, comparison to experiment can prove constructive. Kjærgaard and co-workers have recently measured standard binding free energies for small gas phase compounds and compared them to CCSD(T)/aug-cc-pV(T+d) calculations.  For example, in the case of acetronitrile-HCl the measured binding free energy at 295K is between 1.2 and 1.9 kcal/mol, while the predicted value is 2.3 kcal/mol using the harmonic approximation.  Since the errors in $\Delta E$ and the rigid-rotor approximation presumably are quite low, this suggest and error in the vibrational free energy of at most 1.1 kcal/mol, despite the fact that the lowest vibrational frequency is only about 30 cm$^{-1}$.  Furthermore, the error can be reduced by 0.4 kcal/mol by scaling the harmonic frequencies by anharmonic scaling factors suggested by Shields and co-workers.  Similar results were found for dimethylsulfide-HCl.  So there are some indications that the harmonic approximation yields free energy corrections that are reasonable and that can be improved upon by relatively minor corrections.

Grimme has taken a different approach by arguing that low-frequency modes resemble free rotations and using the corresponding entropy term for low frequency modes.  This changes the RRHO free energy correction by 0.5 - 4 kcal/mol, depending on the system.

Imaginary frequencies. Low frequencies are especially susceptible to numerical error and it is not unusual to see 1 or 2 imaginary frequencies of low magnitude in a vibrational analysis of a host-guest complex.  Since imaginary frequencies are excluded from the vibrational free energy this effectively removes 1 or 2 low frequency contributions to the vibrational free energy. For example, a 30 cm$^{-1}$ frequency contributes about 1.7 kcal/mol to the free energy at 298K.

Imaginary frequencies resulting from a flat PES and numerical errors can often be removed by making the convergence criteria for the geometry optimization and electronic energy minimization more stringent and making the grid size finer in the case of DFT calculations. If the Hessian is computed using finite difference it is important to use double-differencing.  If all else fails, it is probably better to pretend that the imaginary frequency is real and add the corresponding vibrational free energy contribution.

Conformations. One of the main problems in computing accurate binding free energies is to identify the structures of the host, guest and (especially) the host-guest complex with the lowest free energy. Because both the RRHO and solvation energy contributions contribute greatly to the binding free energy change, simply finding the structure with the lowest electronic energy and computing the free energy only for that conformation is probably not enough.

The free energy for a molecule (X) with N conformations the standard free energy is

$G^\circ(\mathrm{X})= G^\circ _0(\mathrm{X})-RT\ln\left(\sum^N_{i=1} \left( 1+\exp{\left( -(G^\circ _i- G^\circ _0)/RT\right)} \right)\right)$

where $G^\circ _0(\mathrm{X})$ is the conformation with lowest free energy.  Conformations with free energies higher than 1.36 kcal/mol contribute less than 0.1 to the sum at 298K.  So a significant number of very low free energy structures is needed to make even a 0.5 kcal/mol contribution to the free energy.  Conformations related by symmetry should not be included here as their effects are accounted for in the rotational entropy (see above).  Also, conformationally flexible regions of the ligand or host that are unaffected by binding need not be explored since the effect on the binding free energy will cancel.

This work is licensed under a Creative Commons Attribution 4.0

QM computed standard free energy changes and pH

2015.01.25:  I have summarized discussion of this and related issues in this paper.  I have also made some changes in the post to reflect what I have learned since I wrote it.

The problem
Thermodynamics involving ionizable functional group at constant pH require special considerations. Robert Alberty has written extensively on this topic (example) but I confess I had some trouble relating the work to quantum chemical calculations.  This post aims to do just that.

Let's say you want to compute the standard free energy change at pH 7 for this equilibrium

$\mathrm{B \rightleftharpoons HA}$

where molecule $\mathrm{B}$ can be converted to an acid $\mathrm{HA}$, which is in equilibrium with it conjugate base and thus pH-dependent.

$\mathrm{HA \rightleftharpoons A^- + H^+}$

The apparent equilibrium constant ($K^\prime$) is

$K^\prime=\mathrm{\frac{[HA]+[A^-]}{[B]}}$

How to compute the corresponding standard free energy change?:

$\Delta G^{\prime\circ}=-RT\ln(K^\prime)$

The Solution
$K^\prime$ can be rewritten as

$K^\prime=K+\frac{KK_a}{[\mathrm{H^+}]}$

where

$K=\mathrm{\frac{[HA]}{[B]}}=\exp{(-\Delta G^\circ/RT)}$

and

$K_a=\mathrm{\frac{[A^-][H^+]}{[HA]}}=\exp{(-\Delta G^\circ_a/RT)}$

Thus,

$K^\prime=\exp{(-\Delta G^\circ/RT)} + \exp{(-\Delta G^\circ/RT)}\exp{(-\Delta G^\circ_a/RT)}10^{\mathrm{pH}}$

$=\exp{(-\Delta G^\circ/RT)} + \exp{(-(\Delta G^\circ+\Delta G^\circ_a-RT\ln(10)\mathrm{pH})/RT)}$

$=\frac{\exp{(-G^\circ(\mathrm{HA})/RT)}+\exp{(-G^{\prime\circ}(\mathrm{A^-})/RT)} }{\exp{(-G^\circ(\mathrm{B})/RT)}}$

where

$G^{\prime\circ}(\mathrm{A^-})=G^{\circ}(\mathrm{A^-})+[G^{\circ}(\mathrm{H^+})-RT\ln(10)\mathrm{pH}]$  (1)

Thus

$\Delta G^{\prime\circ}= G^{\prime\circ}(\mathrm{\overline{HA}})-G^\circ(\mathrm{B})$ (2)

where

$G^{\prime\circ}(\mathrm{\overline{HA}})=-RT\ln \left( \exp{(-G^\circ(\mathrm{HA})/RT)}+\exp{(-G^{\prime\circ}(\mathrm{A^-})/RT)}  \right)$ (3)

Here $\mathrm{\overline{HA}}$ refers to "$\text{HA and A}$"

An example
How do you compute the standard free energy change in aqueous solution at pH 7 using QM for this reaction?:

$\mathrm{\text{HC(=O)-}NH_2+ H_2O \rightleftharpoons HCOOH + NH_3}$

First you optimize each molecule and perform the vibrational analysis to get the free energies.  You can account for solvation effects using a method like PCM or COSMO.  $\mathrm{HCOOH}$ and $\mathrm{NH_3}$ are acids and basis, respectively so you will also need the energies for $\mathrm{HCOO^-}$ and $\mathrm{NH_4^+}$

Most QM programs assume the molecules are in the gas phase when computing the Gibbs free energies so the standard state is 1 bar.  The free energies must be corrected for the solution standard state of 1 M:

$G^\circ(\mathrm{X}) = G^{\circ,\mathrm{1 bar}}(\mathrm{X})-RT\ln (V^{-1})$

where $V$ is the molar volume of an ideal at gas at your chosen temperature.  The exception is water since that is also the solvent.  Here the standard state is 55.34 M at 25C.

$G^\circ(\mathrm{H_2O}) = G^{\circ,\mathrm{1 bar}}(\mathrm{H_2O})-RT\ln ((55.34V)^{-1})$

$G^{\prime\circ}(\mathrm{HCOO^-})$ and is then computed using Eq (1). $G^{\circ}(\mathrm{H^+})$ is usually taken from the literature though estimates vary between -265.8 and -268.6 kcal/mol.

$G^{\prime\circ}(\mathrm{\overline{HCOOH}})$ can then be computed using Eq (3).  The electronic energy component of $G^{\circ}(\mathrm{X})$ can be quite large in magnitude and give some numerical problems when computing the exponential function so Eq (3) and be rewritten as

$G^{\prime\circ}(\mathrm{\overline{HA}})= G^\circ(\mathrm{HA})-RT\ln \left(1+\exp{(-(G^{\prime\circ}(\mathrm{A^-})-G^\circ(\mathrm{HA}))/RT)}  \right)$

Following a derivation similar to the one at the beginning of this post,

$G^{\prime\circ}(\mathrm{NH_4^+})=G^{\circ}(\mathrm{NH_4^+})-[G^{\circ}(\mathrm{H^+})-RT\ln(10)\mathrm{pH}]$

and $G^{\prime\circ}(\mathrm{\overline{NH_3}})$ is computed by Eq (3)

Finally, we put everything together

$\Delta G^{\prime\circ}=G^{\prime\circ}(\mathrm{\overline{HCOOH}})+G^{\prime\circ}(\mathrm{\overline{NH_3}})- G^\circ(\mathrm{\text{HC(=O)-}NH_2})-G^\circ(\mathrm{H_2O})  $

Another way
If you happen to know the pKa values (e.g. from experiment) of the ionizable species you can use this expression

$\Delta G^{\prime\circ}=\Delta A^{\circ}-RT\ln \left( 1+10^{\mathrm{pH}- pK_{a,\mathrm{HCOOH}}} \right)-RT\ln \left( 1+10^{ pK_{a,\mathrm{NH_4^+}}-\mathrm{pH}} \right)$

where

$\Delta G^{\circ}=G^{\circ}(\mathrm{HCOOH})+G^{\circ}(\mathrm{NH_3})- G^\circ(\mathrm{\text{HC(=O)-}NH_2})-G^\circ(\mathrm{H_2O})$



This work is licensed under a Creative Commons Attribution 4.0

Calculation of Solvation Free Energies of Charged Solutes Using Mixed Cluster/Continuum Models

2015.01.25:  I have summarized discussion of this and related issues in this paper.  Point 1 is wrong and Gibbs free energies should be used throughout.

This post takes it title from this paper by Bryantsev et al., which provides an excellent discussion of the issue.  The point of this post is to suggest four modifications to the one of the approaches discussed in the paper:

1. The use of standard Helmholz free energies instead of standard Gibbs free energies

2. Correction for indistinguishability of water molecules

3. Parameterization of a water oxygen solvent radius

4. The general use of the cluster model

Background
Continuum solvation models do not account for strong and specific interactions of solvent molecules with the solute (notably ions) so one must include some explicit solvent molecules in the continuum calculations.  The question is how many explicit solvent molecules to use.

Pliego and co-workers suggested that the optimum number of water molecules is the one for which $\Delta G^*_{\text{solv}}(\text{A}^{m\pm})$, computed using Scheme 1, is most negative (refer to the paper for a definition of terms).

Taken from 10.1021/jp802665d. (c) American Chemical Society

However, Bryantsev et al. argue that $\Delta G^*_{\text{solv}}(\text{A}^{m\pm})$ should be computed according to Scheme 2 and that computed in this way, $\Delta G^*_{\text{solv}}(\text{A}^{m\pm})$ will converge smoothly to the best value.  Thus, the optimum number of explicit solvent molecules is that for which $\Delta G^*_{\text{solv}}(\text{A}^{m\pm})$ changes very little when you add another one.

In principle, Scheme 1 and 2 should yield the same $\Delta G^*_{\text{solv}}(\text{A}^{m\pm})$ but Bryantsev et al. argue that Scheme 2 provides better cancellation of error compared to Scheme 1. The authors demonstrate this by point showing that the computed free energy of water cluster formation in bulk water, computed using Scheme 3, deviates significantly from 0 with increasing cluster size.

Taken from 10.1021/jp802665d. (c) American Chemical Society

In Figure 1 I show a plot of the residual error ($RE$, large dots), 

$RE=-n\Delta G^*_{\text{solv}}(\text{H}_2\text{O}) + \Delta G^\circ_{\text{g,bind}}-(n-1)\Delta G^{\circ-*}$

$ + \Delta G^*_{\text{solv}}((\text{H}_2\text{O})_n)-RT\ln (n[\text{H}_2\text{O}]^{n-1})$

computed using two different continuum solvation models (SM6 red and COSMO blue) for increasing number of $n$.

Figure 1. Plot of RE without (large dots) and with (small dots) the Helmholz free energy and correction for indistinguishability of water molecules

Next I'll suggest some changes that significantly reduce $RE$. 

1. The use of standard Helmholz free energies instead of standard Gibbs free energies

The change in Helmholz free energy ($\Delta A^\circ$) is related to the change in Gibbs free energy ($\Delta G^\circ$) by

$\Delta G^\circ=\Delta A^\circ+p^\circ \Delta V$

In the gas phase $p^\circ \Delta V = \Delta n RT$ where $\Delta n$ is the difference between the number of product and reactant molecules.  In solution $\Delta V$ is hard to estimate, but the best approximation is $\Delta V \approx 0$ and $\Delta G^\circ \approx \Delta A^\circ$.  So, I suggest using 

$\Delta A^\circ_{\text{g,bind}}=\Delta G^\circ_{\text{g,bind}}-(n-1)RT$

2. Correction for indistinguishability of water molecules

Let's consider formation of the water dimer in solution, i.e. Scheme 3 for $n=2$.  As I have written about earlier, there are two ways of making the water dimer: one in which molecule $A$ is the H-donor and one where molecule $B$ is the H-donor (I assume the gas phase Hessian calculation for the water molecule is done in $C_{2v}$ symmetry).  In general, a (H$_2$O)$_n$ cluster can be made $n!$ ways so I suggest using 

$\Delta A^\circ_{\text{g,bind}}=\Delta G^\circ_{\text{g,bind}}-(n-1)RT-RT\ln(n!)$ 

With these corrections $RE$ vs $n$ is significantly reduced as shown in Figure 1 (small dots). For example, for $n=18$ and COSMO $RE$ has been reduced from 41 to 10 kcal/mol.  The remaining error comes from errors in the solvation energies, hydrogen bond strengths, anharmonicity, and conformational entropy.  The $n=18$ cluster has around 30 hydrogen bonds to the 10 kcal/mol error could easily be explained by a 0.3 kcal/mol error in the computed hydrogen bond strength. 

3. Parameterization of a water oxygen solvent radius

The solvation energies are clearly another source of error.  For example, the COSMO value for $\Delta G^*_{\text{solv}}(\text{H}_2\text{O})$ is 0.35 kcal/mol lower than the experimental value, corresponding to an error of 6.3 kcal/mol for 18 water molecules.  Part of that error is probably cancelled by similar errors in $\Delta G^*_{\text{solv}}((\text{H}_2\text{O})_n)$ but not all.  

In fact, the higher errors for SM6 must be due to its underestimating $\Delta G^*_{\text{solv}}(\text{H}_2\text{O})$ by 2.51 kcal/mol (corresponding to ca 45 kcal/mol in the gas phase for $n=18$!).  Thus if SM6 has to be used it is important to re-parameterize it so that $\Delta G^*_{\text{solv}}(\text{H}_2\text{O})$ is reproduced reasonably well.  The easiest way is probably to change the radius of the oxygen water molecule sphere.

Notice that results are not improved simply by using the experimental value of $\Delta G^*_{\text{solv}}(\text{H}_2\text{O})$ instead of the computed one.  For example, for $n=18$ the difference in $RE$ for SM6 and COSMO (where the error in the solvation energy is quite small) is "only" ca 20 kcal/mol so introducing a 45 kcal/mol correction will not improve things for the right reasons. 

Changes in Scheme 1 and 2
So, I suggest the following corrections to Scheme 1 and 2, respectively

$\Delta A^\circ_{\text{g,bind}}(I)=\Delta G^\circ_{\text{g,bind}}(I)-nRT-RT\ln(n!)$ 

$\Delta A^\circ_{\text{g,bind}}(II)=\Delta G^\circ_{\text{g,bind}}(II)-RT$ 

Notice that there is a $RT\ln(n!)$ term for both $(\text{H}_2\text{O})_n$ and $[\text{A}(\text{H}_2\text{O})_n]^{m\pm}$ leading to cancellation of this term for Scheme 2.

Clearly the corrections are much more important for Scheme 1 and Scheme 2 is still the more accurate approach due to better cancellation of errors.

4. The general use of the cluster model
While Bryantsev et al. advocate the use of the cluster cycle to find the optimum cluster size for the ions they still seem to advocate the use of a monomer cycle for pKa predictions (Scheme 4).

Taken from 10.1021/jp802665d. (c) American Chemical Society

I think better cancellation would be achieved by using two water clusters of size $n$ and $m$ instead of $n+m$ water monomers.  Thus

$\Delta G^\circ_{\text{g,deprot}}-(n+m-1)\Delta G^{\circ-*} \rightarrow \Delta A^\circ_{\text{g,deprot}} - \Delta G^{\circ-*}$

$\Delta A^\circ_{\text{g,deprot}}=\Delta G^\circ_{\text{g,deprot}}-RT$ 

$(n+m)\Delta G^*_{\text{solv}}(\text{H}_2\text{O}) \rightarrow \Delta G^*_{\text{solv}}((\text{H}_2\text{O})_n)+ \Delta G^*_{\text{solv}}((\text{H}_2\text{O})_m)$

$(n+m)RT\ln([\text{H}_2\text{O}]) \rightarrow 2RT\ln([\text{H}_2\text{O}]/(nm))$

This work is licensed under a Creative Commons Attribution 4.0

Protein unfolding at high temperature happens because of entropy

Here is a video in which I use Molecular Workbench to illustrate protein unfolding.  I then go on to explain how conformational entropy contributes to the spontaneous unfolding with increasing temperature. The videos are part of a series that I am working on.


This work is licensed under a Creative Commons Attribution 4.0 International License.

MolCalc and the ideal gas enthalpy and entropy contributions

Here are two videos in which I introduce the enthalpy and entropy contributions for an ideal gas with two concrete examples illustrated using MolCalc. The videos are part of a series that I am working on.


This work is licensed under a Creative Commons Attribution 4.0 International License.

Illustrating enthalpy, entropy, and free energy changes using MD

Here is a video I made that uses an MD simulation I found on Youtube to illustrate enthalpy, entropy, and free energy changes. I also use the simulation in these two videos (here and here) to illustrate enthalpy and entropy changes separately.  The videos are part of a series that I am working on.



This work is licensed under a Creative Commons Attribution 4.0 International License.

Illustrating energy states and estimating enthalpy changes

Here are two videos in which I use Jmol and Molecular Workbench to illustrate energy states and MolCalc to estimate enthalpy changes. The videos are part of a series that I am working on.


This work is licensed under a Creative Commons Attribution 4.0 International License.

Illustrating dynamic equilibrium

Here is a video in which I use Molecular Workbench to illustrate dynamic equilibrium.  The video is part of a series that I am working on.


This work is licensed under a Creative Commons Attribution 4.0 International License.

Chemistry assignments that use Molecule Calculator (MolCalc)

1. One of the reviewers of our J. Chem. Ed. paper on MolCalc included the following tutorial: Molecular Orbital Calculations of Molecules I.Diatomics, Triatomics and Reactions

2.  n-Butane can exist in two different conformations called gauche and anti (Google butane and conformation).  Use Molecular Calculator to estimate the fraction of molecules in the gauche conformation at 25 $^\circ$C. $\Delta H^\circ$  can be computed as the difference in heat of formation.

3. Estimate $\Delta H^\circ$  the for the following reaction at 25 $^\circ$C

NH$_2$CHO + H$_2$O $\rightleftharpoons$ NH$_3$ + HCOOH

a. Using bond energies
b. Using Molecule Calculator

4. How does the molecular structure determine the rotational entropy?  Find out by constructing a molecule with the largest possible rotational entropy using Molecule Calculator.  The largest value I could find was 133 J/molK.  Can you beat that?

5. How well do the simple solvation models work?
a. Estimate the solvation energy of NH$_4^+$ using MolCalc?
b. What is the polar solvation energy of NH$_4^+$ in water at 25 $^\circ$C assuming that it is spherical?

6. Why do ionic compounds dissolve in water?  Use MolCalc to estimate $\Delta G^\circ$ at 25 oC for the following equilibrium 

 

N(CH$_3$)$_4^+\cdot$Cl$^-$ $\rightleftharpoons$ N(CH$_3$)$_4^+$ + Cl$^-$

a. in the gas phase

b. in aqueous solution

This blogpost contains some useful hints for part b. 

7. Solvent screening: charge-charge interactions are weaker in aqueous solution than in the gas phase.  Compute the difference in G$^\circ$ at 25 $^\circ$C between these two molecules using MolCalc

  

a. in the gas phase

b. in aqueous solution

This blogpost contains some useful hints for part b. 

8. Build a molecule with a solvation energy that is as close to 0 as possible.  The closest I got is -1.3 kJ/  How close can you get? 

If you have other suggestions please leave a comment

This work is licensed under a Creative Commons Attribution 3.0 Unported License.

Conformational and rotational entropy and point group symmetry

The point group affects the free energy
The use of point group symmetry in quantum chemical calculations can speed up calculations significantly, but it is often difficult to input symmetric coordinates correctly so many people opt to run calculations without symmetry (i.e. in $C_1$ symmetry) anyway.  However, the lack of symmetry changes the rotational entropy and, hence, the free energy you compute. So if you have symmetric molecules but choose to run in $C_1$ symmetry you must correct the entropies and free energies.

The rotational entropy is given by$$S_{rot}=R\ln\left(\frac{8\pi^2}{\sigma}\left(\frac{2\pi ekT}{h^2}\right)^{3/2}\sqrt{I_1I_2I_3}\right)$$ $\sigma$ is called the symmetry number and depends on the point group: for example, $\sigma=1$ for $C_1$ and $C_s$, $\sigma=n$ for $C_{nv}$ and $C_{nh}$, $\sigma= 2n$ for $D_{nh}$ and $D_{nd}$, and $\sigma= 12$ or $T_d$.  So the rotational entropy calculated with and without symmetry will differ by $R\ln(\sigma)$:$$S_{rot}=S_{rot}^{C_1}-R\ln(\sigma)$$and similarly for the free energy$$G^\circ=G^{\circ,C_1}+RT\ln(\sigma)$$
An example: $H_2O+Cl^-\rightleftharpoons HOH\cdots Cl^-$
The free energy change for this reaction is$$\Delta G^\circ=\Delta G^{\circ,C_1}+RT\ln\left(\frac{\sigma_{HOH\cdots Cl^-(C_s)}}{\sigma_{H_2O(C_{2v})}\sigma_{Cl^-(C_1)}}\right)\\ \Delta G^\circ=\Delta G^{\circ,C_1}+RT\ln\left(\frac{1}{2\times 1}\right)=\Delta G^{\circ,C_1}-RT\ln (2)$$The corresponding equilibrium constants are$$K=e^{-\Delta G^\circ/RT}=2K^{C_1}$$So, $\sigma$'s accounts for the fact that, because of the symmetry of water, there are two ways of making $HOH\cdots Cl^-$and another way to view $R\ln(2)$ is that it is the conformational entropy of the complex.

A test: $H_2O+NH_3\rightleftharpoons HOH\cdots NH_3$
For the above equilibrium what is $X$ in$$\Delta G^\circ=\Delta G^{\circ,C_1}-RT\ln (X)$$

     

An exception: $2H_2O\rightleftharpoons HOH\cdots OH_2$
Based on the rules outlined so far one would expect the free energy change for this equilibrium to be$$\Delta G^\circ=\Delta G^{\circ,C_1}-RT\ln (4)$$The factor of four accounts for the fact that there are four ways of making $HO^AH\cdots O^BH_2$. However, since the water molecules are identical there are actually four additional water dimer possibilities for  $HO^BH\cdots O^AH_2$, so$$\Delta G^\circ=\Delta G^{\circ,C_1}-RT\ln (8)$$In general, for $A+A\rightleftharpoons Product$ reactions the symmetry number for $A+A$ is $2\sigma_A^2$rather than $\sigma_A^2$

All these considerations also apply to activation free energies and rate constants as outlined in this excellent paper by Fernández-Ramos et al., which inspired this post.  See also this excellent paper by Gilson and Irikura.


This work is licensed under a Creative Commons Attribution 3.0 Unported License.

Thermodynamics in solution: a brief guide for quantum chemists

2015.01.25: Please read this paper instead of this post. It turns out most of what I write here is wrong.

The thermodynamic properties such as enthalpy, entropy, and free energy you get from a vibrational analysis by programs such as GAMESS and Gaussian are those of an ideal gas at 1 bar pressure (and usually 298 K).  If you calculate free energy changes in solution (using methods such as PCM or COSMO) there are a few things you should do differently compared to gas phase calculations.  These things will only make a difference if you are computing free energy changes for processes where the number of particles change, such as binding free energies.  The corrections will cancel out for things such as conformational free energy differences.

Use Helmholtz free energies instead of Gibbs free energies
Experimental studies typically report Gibbs free energy changes ($\Delta G^\circ$), which are related to Helmholtz free energy changes ($\Delta A^\circ$) by$$\Delta G^\circ = \Delta A^\circ + p^\circ \Delta V$$ $\Delta V$ is the change in volume of the solution due to the reaction.  This volume change is negligible so $\Delta G^\circ = \Delta A^\circ$ is a good approximation.

The Gibbs free energy printed by quantum chemistry program correspond to an ideal gas where $pV=RT$ and the difference between $\Delta G^\circ$ and $\Delta A^\circ$ is much larger.

Computing free energies with semi-empirical methods
Most semiempirical methods such as AM1 and PM6 are parameterized such that the electronic energy matches experimental heat of formations ($\Delta H_f^\circ$) at 298 K.  So a gas phase free energy change should be computed as$$\Delta G^\circ=\Delta\Delta H_f^\circ-T\Delta S^\circ$$i.e. you don't need any of the enthalpy information printed out as part of the vibrational analysis.  A solution free energy change should be computed as $$\Delta A^\circ=\Delta\Delta U_f-T\Delta S^\circ$$where $\Delta U_f=\Delta H_f^\circ-RT$

Added 2013.08.16: However, dispersion and/or hydrogen bond corrected semi-empirical methods such as PM6-DH+ are parameterized against electronic binding energies.  So if you are computing binding free energies with such methods you need to add the translational, rotational, and vibrational enthalpy corrections to  $\Delta\Delta H_f^\circ$.

Change the standard state to 1 mol per liter
The translational entropy printed out depends on the volume of the system, which is computed as $V=RT/p$ = 24.79 liters.  For solution that volume should be 1 liter.  Most programs do not allow you to change volume so you must apply the correction to the entropy$$S_{soln}^\circ=S_{gas}^\circ+R\ln\left(\frac{1}{24.79}\right)$$ yourself

Is it OK to use the rigid rotor-harmonic oscillator approximation in solution?
Short answer: yes.  The derivations of the translational, rotational, and vibrational free energies are done for an isolated molecule in vacuum, while a molecule in solution interacts with solvent molecules.  However, a molecule in solution still has the exact same degrees of freedom as in the gas phase: it is free to explore the entire three dimensional volume and is free to adapt any rotational angle just like in the gas phase, so the resulting free energy expressions are the same.  Similarly, the individual molecules still have $3N-6$ internal vibrational frequencies.  The energy from the solute-solvent vibrations are included in the solvation free energy.

This blogpost was inspired by this excellent paper by Muddana and Gilson, and I thank Mike Gilson for helpful discussions.


This work is licensed under a Creative Commons Attribution 3.0 Unported License.

Simulations in teaching physical chemistry: thermodynamics and statistical mechanics

In this post I summarize the simulations and I have used in teaching thermo and stat mech, and talk a bit about how I use them.

I co-teach two quite similar courses on this topic: one for nano-students and another for chemistry and biochemistry students.  In the nano course we use the book Molecular Driving Forces by Dill and Bromberg, and in the other Quanta, Matter, and Change by Atkins, de Paula, and Friedman.  At the end of this post I have organized the simulations by chapter for each book.

Some of simulations I have made (or modified extensively) and most of these have been discussed in previous blog posts, so I simply give the link to the respective blog post where there is more information.

The other simulations are from the Molecular Workbench (MW) library of models, and here I provide links that will open in MW, so you need to install MW before clicking on the links.  For some of them I also provide a brief description of what concepts try to demonstrate using the simulations.

How do I use the simulations?
All simulations are used during lecture to visualize concepts, start discussions, and motivate equations. I'll take Illustrating energy states as an example: instead of saying "Molecules in a gas translate, rotate, vibrate, and ....", I say "Here is a zoomed-in view of butane gas where you can see the molecules.  You can see that individual molecules move differently.  How do they move differently?  Anyone?  Right, they have different speeds.  This kind of motion is called translation.  What else? ..."

Practical tips
On a very practical note, my own simulations are all on web sites and I make sure to open all of them before the lecture, while I have all the MW simulations for the course indexed on a single MW page (click here to open in MW). It is not possible to embed these simulations in Powerpoint slides, but you can switch between Powerpoint and other applications without quitting Powerpoint (on a Mac you use command-tab and on Windows i believe it is windowskey-tab).  Note that you need access to the internet in the lecture room.

While I have screencasts of most of simulations on the blog posts, I don't use these during lecture.  I think it is too passive, and puts the students to sleep.  But I believe the screencasts are a good way for the students to review the main points of simulations after the lecture.  I put links to the blog posts on the course web site and in the lecture notes.

Is using simulations a good idea?
If possible I try to use a simulation within the first five minutes of a lecture, and have a maximum of 20 minutes between simulations.  I now only have one (45 minute) lecture left where I don't use a single simulation and I can just feel how I loose the student's attention after about 30 minutes.  You can just see it.  That being said, no one has ever mentioned the simulations in their course evaluations (good or bad), so I have no hard evidence that it improves my teaching.  But I can tell you that I enjoy lecturing much more with the simulations, so unless I get complaints I'll keep doing it. 

Making room for simulations in the lecture
I have taught the topics for many years without any simulations, and was never at a loss for material to cover.  Lecture time is precious, and these simulations take time to present and discuss.  You really have to introduce the simulation carefully (don't rush this part!) before you start them, and very often you want the students to speculate about what will happen before you start them.  Furthermore, they tend to stimulate many more questions, that you can hopefully turn into a discussion instead of simply answering them, than derivations - that's the whole point.

So how do you "make room" for the simulations?  I have cut out most of the derivations from the lectures.  To pay for my sins, I provide relatively detailed (typed) lecture notes ahead of lecture (I generally don't use Powerpoint), which include step-by-step derivations. So I'll say things like "Starting with these assumptions we can write down this equation.  This can be rewritten as this equation, which is much simpler.  The details on how we got from here to there are in your notes, but note that in step 3 we assume that ... which is an approximation."  No complaints so far.  If only more progress had been made on simulating derivations ...

Here are the simulations organized by chapter

Molecular Driving Forces by Dill and Bromberg (1st edition)

Ch 6: Entropy and the Boltzmann distribution law
Illustrating entropy

Ch 10: Boltzmann distribution law
Polymer unfolding: The book uses two simple bead models of polymers in this chapter to illustrate micro and macrostates and model protein melting.  I use this example extensively both in lectures and homework problems.  So I made this simulation to illustrate how higher energy macrostates become more likely at higher temperatures.


Ch 11: Statistical mechanics of simple gasses and solids
Illustrating energy states
Energy states in the water molecule: a slightly more complicated molecule than HCl (used in Illustrating energy states) with more than one vibrational mode and 3 rotational degrees of freedom.
Internal energy and molecular motion
Entropy, volume, and temperature


Ch 12: Temperature, heat capacity
The molecular basis of differential scanning calorimetry: heat capacity and energy fluctuations

Ch 13: Chemical equilibria
Seeing chemical equilibrium (opens in MW)
Dalton's law of partial pressure (opens in MW)


Ch 14: Equilibria between solids, liquids, and gasses
Seeing specific and latent heat (opens in MW): I use this simulation to illustrate how the same substance can be solid, liquid, and gas depending on the temperature.
A gas under a piston (opens in MW): I use this simulation to show that, for example, decreasing the pressure can have the same effect as increasing the temperature.
The phase diagram explorer (opens in MW)
Raoult's law: ideal solutions (opens in MW): Here, I use the simulation of the pure liquid to illustrate vapor pressure.


Ch 15: Solution and Mixtures
Mixing gasses, and mixing of ideal and non-ideal liquids
Raoult's law: ideal solutions (opens in MW)
Raoult's law: negative deviation (opens in MW) 
Raoult's law: positive deviation (opens in MW)


Ch 16: Solvation and transfers of molecules between phases
Visualizing osmotic pressure in an osmotic equilibrium (opens in MW)
Desalination using reverse osmosis (opens in MW)




Quanta, Matter, and Change by Atkins, de Paula and Friedman (1st edition)

Ch 13: The Boltzmann distribution
Illustrating energy states
Energy states in the water molecule: a slightly more complicated molecule than HCl (used in Illustrating energy states) with more than one vibrational mode and 3 rotational degrees of freedom.
Internal energy and molecular motion

Ch 14: The first law of thermodynamics#
The molecular basis of differential scanning calorimetry: heat capacity and energy fluctuations
  
Ch 15: The second law of thermodynamics
Illustrating entropy
Entropy, volume, and temperature
  
Ch 16: Physical equilibria
Seeing specific and latent heat (opens in MW): I use this simulation to illustrate how the same substance can be solid, liquid, and gas depending on the temperature.

A gas under a piston (opens in MW): I use this simulation to show that, for example, decreasing the pressure can have the same effect as increasing the temperature.

The phase diagram explorer (opens in MW)
Raoult's law: ideal solutions (opens in MW): Here, I use the simulation of the pure liquid to illustrate vapor pressure.

Raoult's law: negative deviation (opens in MW) 
Raoult's law: positive deviation (opens in MW)

Mixing gasses, and mixing of ideal and non-ideal liquids 

Visualizing osmotic pressure in an osmotic equilibrium (opens in MW)
Desalination using reverse osmosis (opens in MW)

Ch 17: Chemical equilibria#
Seeing chemical equilibrium (opens in MW)
Dalton's law of partial pressure (opens in MW)

# I don't teach this part of the course, but if I did I would use these simulations

Related posts:
An Atkins Diet of Molecular Workbench 
One, Two, Three, MD 
Tunneling and STM (a first stab at using Molecular Workbench to teach quantum mechanics)

Illustrating mixing

This screencast shows Molecular Workbench simulations I have made to illustrate mixing.

The first simulation illustrates the mixing of 2 ideal gases, which mix readily.  Since the gas particles don't interact you can think of the mixing as each gas expanding to fill both containers independently of each other.  As I have shown in this simulation, the driving force for this expansion is an increase in entropy.  Therefore, the driving force for mixing two ideal gasses is also purely entropic.

The second set of simulations illustrates the mixing of 2 liquids.  Since they are liquids there must be attractive interactions between the atoms.  If there were no interactions they would be (ideal) gasses.  The strength of the interactions (and the temperature) determine whether they mix or not.

In the first liquid simulation, the attraction between two green atoms (ε_GG), between two blue atoms (ε_BB), and between a green and a blue atom (ε_GB) are the same.

This means that a green atom doesn't care whether it is sitting next to a blue atom or another green atom.  The net effect is that green and red atoms are equally likely to be on the right or left side of the container, and the liquids mix for the same reason as the ideal gasses mix: the driving force is purely entropic. That means the enthalpy of mixing is zero:

This is the definition of an ideal mixture (or ideal solution).  The two liquids will mix at any temperature.

In the second liquid simulation, the attraction between two blue atoms (ε_BB) is stronger than between two green atoms (ε_GG) and between a green and a blue atom (ε_GB).

Note that the ε's are negative: a smaller ε means a stronger attraction.  This means that the blue particles would rather be with other blue particles, i.e. the enthalpy increases if the particles are mixed.

(z is the number of contacts between particles in solution, and x_G is the mole fraction of green atoms).  This is an example of a non-ideal mixture, where the definition for a non-ideal mixture is

Because Δ_mixH > 0 this non-ideal solution mixes spontaneously (i.e. Δ_mixG < 0) only for

Oil and water is a common example of such an non-ideal mixture: the oil-oil interactions are stronger than the oil-water and water-water interactions.

Salt and water is another example of on idea mixture, but here Δ_mixH <  0 so salt and water almost always mixes spontaneously.  The interpretation is that the interactions between the salt ions and water is stronger than the average interaction between salt ions and between water molecules.

Implications and limitations
The definition of Δ_mixH in terms of the ε's suggest that liquids should also mix if

which would be a more general definition of an ideal mixture.

This is tested in the third liquid simulation.  As you can see the liquids mix more than in the second simulation, but not quite as much as in the first simulation.  This is mostly because of the simulation runs only for 100 picoseconds, which is to short to mix fully.  But another reason is that there is less space (on average) between the blue atoms compared to the green atoms, because the blue atoms attract each other more.  The next effect is that the blue particles tend to stay together to lower the enthalpy.  More mathematically,  z (the number of contacts between particles in solution) is not exactly the same for the blue and green particles so the interpretation of Δ_mixH in terms of the ε's breaks down.  The "safest" definition of an ideal mixture thus remains:

i.e. "like dissolves like".

Accessing the simulations
You can play around with the simulations here and here, or you can download the models here and here if you have Molecular Workbench installed on your computer.