Does MOPAC allows calculation of Standard Enthalpy of reaction in this way? 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)&lt;=&gt;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?
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)$$