Incorrect surface energies using Environ & GC-DFT

[Brad Ayers]

Dear All,

Thank you in advance for any assistance.

I'm encountering an issue when calculating the surface energy of a lithium slab in implicit solvent (SSCS model) with size-modified PB for electrolyte and applied bias. I'm observing small and sometimes negative surface energies.

Using the standard method (G = 1/2A * (Eslab - NBulk)), I get expected results for small developed charges. However, following Alsunni et al. [1], who suggest matching ensembles (canonical-canonical or grand canonical-grand canonical), I applied the correction:

F = E - TS - μN

This changes my surface energy from 0.4 J/m^2 to 1.2 J/m^2. While this resolves negative energies, the value seems unrealistically large. Alsunni et al. found fluctuations of only ± 0.4 J/m^2, whilst I'm seeing values up to 2 J/m^2 for some surfaces.

I've successfully replicated Alsunni's results using their supplementary information, so the issue seems to lie with my Quantum ESPRESSO calculations.

Attached are the surface and bulk files, environ.in, and a Python script for surface energy calculations. I suspect the problem might be in the environ.in or my analysis. Also, frequent warnings about parabolic correction appear - could this be relevant?

Thank you for your help.

Best regards, Brad

[1] https://doi.org/10.1016/j.apsusc.2022.155147

[Brad Ayers]

I wasn't able to attach above so I will do so here :

[Mouyi Weng]

Hi Brad, 

I also encountered a similar energy problem while perfoming Poisson-Boltzmann calculation. In my case, there is a Pt 100 surface. I relaxed the structure with SCCS and linearized PB with -0.1, -0.05, 0, 0.05, 0.1 charge. Then I can get several energies as a function of surface charge. If I plot the total energy (from the output file) verses the surface charge, I can only observe a very linear line, but not in a parabolic shape (as shown in Figure 5 from J. chem. Phys. 150, 041730, 2019)

As the second order term in energy is related to the surface capacity, I am wondering if the interaction energy between the PB electrons and the surface charge is included in the total energy. But I can add a term of Q^2/2C into the total energy by hand. 

Because my project is not very related to the total energy term, so I temporarily ignored this problem. I am not very sure if the above finding is related to your case. But at least they both concern the total energy with PB solver. Perhaps we can wait for Oliviero to help us understand the actual physical meaning of the total energy in the output file. 

Best regards,

Mouyi

[Eric Read]

Could you run the calculations with verbose=1 in Environ.in and share the environ.debug files that are produced. They may shed some light as to what the issue is.