I am using MOPAC2016 to calculate $\Delta H^\circ_\mathrm{r}$ (and then $\Delta G^\circ_\mathrm{r}$) for a reaction in aqueous phase, at $\pu{298 K}$. I will explain how, and then the question. Finally I give an example. Suppose this simple reaction:

\begin{align} \ce{A(aq)<=>B(aq)} \end{align} Then, $\Delta H^\circ_\mathrm{r}$ can be obtained as:

$$\Delta H^\circ_\mathrm{r} = \Delta H_\mathrm{f}^\circ(\ce{B}) - \Delta H_\mathrm{f}^\circ(\ce{A})$$ where each of those values is what we get on a thermodynamic calculation (also in an optimization) with the name of H.O.F. (heat of formation, HOF).

Formally the question should be: Is that correct to calculate $\Delta H^\circ_\mathrm{r}$ (and then $\Delta G^\circ_\mathrm{r}$) by using the value of HOF for each reactant?


Here is an example of MOPAC output for a thermodynamic calculation, suppose it is the one for molecule $\ce{B}$:

                                   KCAL/MOL   CAL/MOLE    CAL/K/MOL  CAL/K/MOL

 298.00  VIB.     0.2155D+06                  8205.6617    53.3353    51.9400
         ROT.     0.3572D+07                   888.2813     2.9808    32.9651
         INT.     0.7697D+12                  9093.9430    56.3161    84.9051
         TRA.     0.4009D+28                  1480.4688     4.9680    42.5424
         TOT.                     -670.232   10574.4119    61.2841   127.4474

              (=  Standard Enthalpy of Formation).

            Hvib:  Zero-point energy is not included.
                   frequencies of less than zero cm-1 are not included.
            Hrot = (3/2)RT
            Htra = (3/2)RT + pV = (5/2)RT

            Heat capacity is Cp, not Cv

          For more information, see: HTTP://OpenMOPAC.net/Manual/thermochemistry.html 

With that data table $\Delta G^\circ_\mathrm{f}$ of the molecule $\ce{B}$ should be: $$\Delta G^\circ_\mathrm{f}= \mathtt{H.O.F.} - TS^\circ= -670.232\times1000 - 298\times(127.4474- S_{atoms})$$ $S_{atoms}$ will be canceled when calculating $\Delta G^\circ_r$. On this mopac site, $\Delta H_\mathrm{f}^\circ$ appears to be obtained in a similar way, for the $\ce{NH4+}$ molecule.

Am I correct?

  • 2
    $\begingroup$ I have no experience with MOPAC, and I don't completely understand the output, but your method looks correct to me. $\endgroup$ Sep 19, 2017 at 9:30
  • $\begingroup$ great, @Martin-マーチン, and thanks for the comment. For some reasons I don't understand most of people working on it don't use this approach.. $\endgroup$
    – user43021
    Sep 19, 2017 at 9:39
  • 1
    $\begingroup$ Semi-empiric calculations are much too inaccurate for many systems and many people use other methods instead. DFT and other ab initio methods allow to calculate an absolute value for the electronic energy and thermodynamic corrections and using those to calculate the relative energies. (The linked manual section is very worth reading, there is more information on how thermochemistry is treated in MOPAC.) $\endgroup$ Sep 19, 2017 at 9:54
  • $\begingroup$ @Martin-マーチン I mean that they use mopac but calculate gibbs free energy in a different way. I didn't specify how because it is long..; this methods aren't too accurate but are parametrized to work good, so, they are quite accurate for most purposes (and also are fast).. $\endgroup$
    – user43021
    Sep 19, 2017 at 10:27
  • 3
    $\begingroup$ $\Delta G_f^\circ$, like, $\Delta H_f^\circ$, is relative to elements in their standard state, so the entropy is missing a term, but that will also cancel when you compute $\Delta G_r^\circ$. So, for 298K, the equation for $\Delta G_f^\circ$ will give the right $\Delta G_r^\circ$. For other temps you also need $(H_T - 10574.4119)$ where $H_T$ is the corresponding total $H$ for temperature $T$. $\endgroup$
    – Jan Jensen
    Sep 19, 2017 at 12:03

1 Answer 1


The free energy change of a reaction occurring in solution can be computed using a thermodynamic cycle: $$\Delta G^\circ_{r,soln} = \Delta G^\circ_{r,gas} + \Delta\Delta G^\circ_{r,solv} $$ $\Delta G^\circ_{r,soln}$ can also be written as $$\Delta G^\circ_{r,soln} = G^\circ_{soln}(B)- G^\circ_{soln}(A)$$ where $$G^\circ_{soln}(X) = G^\circ_{gas}(X) + \Delta G^\circ_{solv} $$ and $$G^\circ_{gas}(X) = E_{gas}(X) + G^\circ_{gas,RRHO}(X)$$ The solvation energy is usually computed using a continuum as the difference in electronic energy with and without the continuum. $\Delta G^\circ_{solv}(X) = E^\prime_{soln}(X)-E^\prime_{gas}(X)$ The parameters in the continuum model are usually optimized for a particular level of theory, which may not be appropriate for the computation of $\Delta G^\circ_{r,gas}$, so $E_{gas}(X)$ is not necessarily the same as $E^\prime_{gas}(X)$

However, if $E_{gas}(X) = E^\prime_{gas}(X)$ then $$G^\circ_{soln}(X) = E_{soln}(X) + G^\circ_{gas,RRHO}(X) $$

The parameters in the continuum model are usually determined using the gas phase geometry $(X)$. But if $E_{gas}(X) = E^\prime_{gas}(X)$ then many people often choose to use the geometry optimized in solution $G^\circ_{soln}(X_{soln}) = E_{soln}(X_{soln}) + G^\circ_{soln,RRHO}(X_{soln})$.

For semiempirical calculations, the electronic energy corresponds to the enthalpy of formation at 298K, so if one simply makes the substitution $$E(X) \rightarrow \Delta H^\circ_{f,298}(X)$$ then (at 298K) $$G^\circ_{soln}(X) = \Delta H^\circ_{f,298}(X) -TS^\circ_{gas,RRHO}(X)$$ $G^\circ_{soln}(X)$ doesn't correspond to a measurable free energy because $S^\circ_{gas,RRHO} \ne S^\circ_{f}$ but it can still be used to compute $\Delta G^\circ_{r,soln}$ since the missing terms cancel.

For temperatures other than 298K $$G^\circ_{soln}(X) = \Delta H^\circ_{f,298}(X) + (H^\circ_{gas,RRHO}(X,T)-H^\circ_{gas,RRHO}(X,T=298K)) -TS^\circ_{gas,RRHO}(X)$$

  • 1
    $\begingroup$ I don't understand why the gas-solution cycle is needed, although I understand the procedure you wrote. But if mopac give us "$G_f$ "for each reactant in solution, why don't we direcly sum those values?(with corresponding index for reactant and products) – $\endgroup$
    – user43021
    Sep 20, 2017 at 17:04
  • 2
    $\begingroup$ the cycle is needed because the continuum model has parameters that need to be fitted to measured values. You can't measure $G^\circ_{soln}(X)$, but you can measure $\Delta G^\circ_{solv}(X)$. So you need the cycle to relate $\Delta G^\circ_{solv}(X)$. to $G^\circ_{soln}(X)$ $\endgroup$
    – Jan Jensen
    Sep 20, 2017 at 17:06
  • $\begingroup$ damn..thats complicated. Thanks for the help.. $\endgroup$
    – user43021
    Sep 20, 2017 at 17:15
  • 3
    $\begingroup$ Yes. Yes, it is. You're welcome $\endgroup$
    – Jan Jensen
    Sep 20, 2017 at 17:21

Your Answer

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