Compute solvation energies with (ESM) RISM-SCF

Question

I want to use the 3D-RISM-SCF and ESM-RISM implementation in Quantum Espresso to compute solvated adsorption energies and absolute solvation energies.

Before starting I read the main paper by S. Nishihara and M. Otani and read (not completely) the code, at least to understand what was printed in the output file.

In the code you can find the most important variable "esol" that is being printed as "solvation energy" in the energy detail. This variable is the sum of an array usol which are the chemical potentials of solvation for each site within the RISM formalism.

It is explicitly stated in the code "only esol has physical meaning not usol".

Everything is perfect, I do my calculations and get my solvation energies. Until recently I have a discussion with one of the authors who tells me that "esol" has no physical meaning because it contains a divergent term that is cancelled by another DFT term.

I read the theory again and I find the supplementary material where it is written that solvation energy (the esol term) diverge only when the unit cell is charged.

I have a questions:

Can someone with RISM experience confirm what I just wrote?

Answer 1

There are certain quantities that become poorly defined when you have a charged system. This could be what the author was referring to if you are studying the adsorption of ions or studying adsorption processes under an applied electrical bias. If you are not studying a charged system, then I believe esol and usol carry physical meaning. However, if you are interested in computing solvation free energies of neutral species, you should compute the energy difference between the solvated and non-solvated system. I believe esol refers to the excess chemical potential of solvation. Absolute solvation free energies for charged species are less well-defined unfortunately since you must account for the air-water dipole potential contribution, which itself is ill-defined and not included in the model. I can double check some of this and get back to you. If this is still an open question, I would also encourage you to post on the QE user forum. You can get a direct (and usually prompt) response from the developers there.

Some more details

The divergence of usol is a consequence of how ESM-RISM solves the Poisson equation of the composite quantum-classical cell. Essentially, the quantum and solvent charge densities are determined self-consistently, where the electron density is iteratively optimized for a frozen set of solvent distributions (solvent charge density), the potential generated by the new electron density is then used to re-optimize the solvent distributions, and this cycle continues until you reach your convergence criterion. As a result, the two charge densities are necessarily determined separately, and the total electrostatic potential of the composite quantum-classical cell is given by the sum of the quantum and classical potentials.

If the quantum region carries a net charge, the electrostatic potentials of the solvent and quantum solute regions must diverge (in opposite directions) since the integrated charges of the solute and solvent carry opposite signs. When summed, the diverging potentials cancel and you obtain a flat potential far from the interface. This is a non-issue though if the quantum region is neutral, since no divergence in the potential occurs.