Does heat of formation include dielectric energy?

[aizvorski]

I'm wondering if the value reported in "FINAL HEAT OF FORMATION" includes the contribution from "DIELECTRIC ENERGY".

The following experiment with an ion in implicit solvent vs vacuum suggests it does, but I'd like to confirm since it doesn't seem the documentation states that either way. Thanks!

Input1:

XYZ PM6-D3H4X 1SCF CHARGE=-1 EPS=78.5


Cl 0 1 0 1 0 1

Output1:

          FINAL HEAT OF FORMATION =       -140.54436 KCAL/MOL =    -588.03762 KJ/MOL
          DIELECTRIC ENERGY       =         -3.48258 EV

Input 2:

XYZ PM6-D3H4X 1SCF CHARGE=-1


Cl 0 1 0 1 0 1

Output 2:

          FINAL HEAT OF FORMATION =        -60.23417 KCAL/MOL =    -252.01975 KJ/MOL

Dielectric energy = -3.48258 EV == -336.018 kJ/mol
Difference in heats of formation = -588.03762 - (-252.01975) kJ/mol == -336.01787 kJ/mol, this matches exactly, suggesting the dielectric energy is included

[godotalgorithm]

I've converted this into a Discussion because it doesn't pertain to a bug or feature request.

Yes, the heat of formation includes a solvation energy that is meant to approximate the molecule-solvent interaction and how the molecule changes the energy of the solvent itself (modeled crudely as charge rearrangements in a dielectric continuum). I would need to review the COSMO model more carefully to provide a more precise answer, but the dielectric energy contains the direct contributions to the heat of formation from solvent effects. In your simple example where symmetry prevents the solvent from causing any charge rearrangements, it is exactly the difference between the heat of formation of the molecule in and out of the implicit solvent. In general, the solvent model will cause charge rearrangements in the molecule, and so the difference between heats in and out of solvent will not match the dielectric energy exactly (although it is usually close).

[aizvorski]

Thanks Jonathan! Yeah, I picked a monoatomic example exactly so there would be no charge redistribution. This is good to know.

P.S. One related question I just thought of: do reported gradients also include the gradient of the dielectric energy wrt coordinates? ie, are they the true gradient of the "FINAL HEAT OF FORMATION" quantity wrt coordinates?

[godotalgorithm]

Yes, the forces should be gradients of the overall heat of formation, including all contributions like implicit solvent. MOPAC's gradients are a form of numerical derivative that use the Hellmann-Feynman theorem with numerical derivatives of the electronic Hamiltonian, which contains the solvent-molecule interactions. MOPAC's numerical gradients are not as accurate as the heat of formation itself, but any appreciable numerical discrepancy would constitute a bug.

While there are some adjustable numerical parameters in MOPAC, it definitely has a "numerical error floor" because of approximations it makes that cannot be systematically improved. In the case of implicit solvent for example, I don't think that numerical first derivatives account for the change in the solvent-accessible-surface mesh caused by perturbing nuclear coordinates. As far as I understand its history, MOPAC's numerical choices have been pragmatic and empirical - they were made for representative sets of practical example calculations at production-level accuracy. If a formal contribution to the numerical derivative (in propagating chains of partial derivatives) was very complicated or expensive but numerically negligible under these conditions, then it was likely ignored. In my own research, I have a lot of experience with implementing and validating numerical and semi-numerical derivatives of very complicated functions, and to be thorough you usually have to create highly sensitive artificial tests and chase terms down to machine precision. MOPAC's numerics are not that rigorous, and it would be difficult to scale that sort of rigor to production software with as many features as MOPAC has.