[Pw_forum] Total energy of solvated molecules

[Oliviero Andreussi]

Dear Robert,

Sorry for the late reply, I have been out of office for a while and have not checked the PWforum mailing list so often in the meantime.

You are right, the output is confusing in that the total energy is not the sum of the reported terms and it is not immediately clear what solvation energy means, but before going into the details of why, I want first to make clear a point which appears to be confused in your question, sorry if I restate something obvious. dGsol, and in particular dGel, are not quantities that you get from a single calculation in solution, but rather are obtained as the difference in energy from a calculation fully optimised (nuclei and electrons) in solution minus a calculation fully optimised (nuclei and electrons) in vacuum. So you need to perform two geometry optimisation calculations to get the dGel (or dGsol) reported in the article, and this is what the examples are set up to perform. So, the "solvation energy" output from a single calculation is not dGsol.

For the solvent contributions, what are reported are the correct expressions of the energies that need to be added to the total energy of a system in vacuum. In particular the "solvation energy" reported is 1/2 integral of ( rho^el + rho^ion) * phi^pol [that is E^pol as reported in Eq. 20-21 of the reference you cite, Andreussi Dabo Marzari JCP 2012], thus the polarisation contribution to the electrostatic energy of the molecule. The cavitation energy is env_surface_tension * Surface, while the pressure energy is env_pressure * Volume. The name "solvation energy" is consistent with what is used in the quantum-chemistry community (in particular in PCM), I agree that is somehow misleading and in most cases the quantity itself is rather useless, it is just an arbitrary component of the total electrostatic energy of the system.

Now, coming to the reason why the different terms don't add up to the correct total energy. This is due to an extra term which is spurious and is subtracted from the reported one-electron contribution, but is not reported in the output. The reason is the following: total energy is computed in pw making use of the secular expression (see for example W.E. Pickett, Computer Physics Reports, 9(3), 115-198, 1989, equations 5.21-5.22 etc), thus from the sum of occupied states eigenvalues (so called band structure energy, "eband" in pw.x), in which some contributions are included correctly (the kinetic and the electron-nuclei interaction), some are double counted (hartree), some need to be removed and added back with the correct expression (xc term), and some are missing (nuclei-nuclei interaction via ewald sum). What PW reports in the output as the one-electron contribution is in fact ebands minus the double counted hartree and minus the wrong xc term (these spurious terms are summed up in pw.x into a term named "deband").

Similarly, when performing the solvated clalculation, there are some terms which are double counted (the electronic part of the solvation energy), some that are missing (the ionic part of solvation energy), some that have the wrong expression (pressure and cavitation) and some that need to be removed (for example a term coming from the rho-dependence of the dielectric constant). All these spurious terms are collected in a term named "deenviron", which has the same meaning of "deband" seen above. To go to the point, what is wrong in the output is the one-electron contribution, which is reported including the spurious deenviron term. As the Environ module was designed not to affect the standard printout of QE, I did not modify the one-electron output, with the result that the terms do not sum to the correct total energy. Reporting deenviron, similarly to reporting deband, does not make too much sense, as there are no clear information into these terms, which collect several different contributions. Maybe I can add an extra line in the Environ output reporting the "true-one-electron contribution =", but this also depends on the future development of the Environ inside QE. Still, I apologise for the confusion this may have caused.

I am happy to provide more details on the output, implementation or usage of the Environ module, through the PW forum mailing list or privately, fee free to contact me directly in case something is not clear or running as expected.

If I misreported something on how PW does things or if anyone else in the community finds something unclear/wrong in what I reported, please let me know.

Thanks for your feedbacks,

Oliviero Andreussi
--
Senior Postdoctoral Researcher
École Polytechnique Fédérale de Lausanne (EPFL) and
Università della Svizzera Italiana (USI) of Lugano
USI Campus, Via G. Buffi 17, 6904 Lugano, Switzerland
Emails: oliviero.andreussi @ epfl.ch -or- usi.ch
Tel: +41-(0)58-666-4810 / Skype: olivieroandreussi
Web: https://sites.google.com/site/olivieroandreussi

________________________________
From: pw_forum-bounces at pwscf.org [pw_forum-bounces at pwscf.org] on behalf of Robert Wexler [rwexler at sas.upenn.edu]
Sent: Wednesday, August 31, 2016 4:41 PM
To: PWSCF Forum
Subject: [Pw_forum] Total energy of solvated molecules

Dear QE users,

I just worked example 1 in espresso-5.3.0/Environ/examples/ and obtained the following result for a water molecule in water solvent:

!    total energy              =     -34.06234148 Ry
     Harris-Foulkes estimate   =     -33.95165940 Ry
     estimated scf accuracy    <          6.8E-12 Ry

     The total energy is the sum of the following terms:

     one-electron contribution =     -81.79999129 Ry
     hartree contribution      =      42.43346037 Ry
     xc contribution           =      -8.39845654 Ry
     ewald contribution        =      13.81332806 Ry
     solvation energy          =      -0.03439636 Ry
     cavitation energy         =       0.01360379 Ry
     PV energy                 =      -0.00691218 Ry

When I sum the energy components, I get -33.97936415 Ry, which is not equal to the total energy above. I read

O. Andreussi, I. Dabo and N. Marzari, J. Chem. Phys. 136 064102 (2012)

and the authors state that dGsol = dGel + Gcav + a*S + b*V + PdV where a and b are fitting parameters to reproduce the experimental solvation energies for a set of reference molecules. Is the "solvation energy" above the same as dGsol? Is the total energy the sum of these terms or does the "solvation energy" contain the cavitation and PV energies?

Thanks in advance for your assistance!

Best,
Rob
---------------------------------------------
Robert B. Wexler
PhD Candidate in Chemistry
University of Pennsylvania, 2018
rwexler at sas.upenn.edu<mailto:rwexler at sas.upenn.edu>
(215) 801-8741
---------------------------------------------