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?
Example
Here is an example of MOPAC output for a thermodynamic calculation, suppose it is the one for molecule $\ce{B}$:
CALCULATED THERMODYNAMIC PROPERTIES
*
TEMP. (K) PARTITION FUNCTION H.O.F. ENTHALPY HEAT CAPACITY ENTROPY
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
*: NOTE: HEATS OF FORMATION ARE RELATIVE TO THE
ELEMENTS IN THEIR STANDARD STATE AT 298K
(= 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?