Technical questions re: using Environ to calculate delta SCF PBE electron affinities

[Edward Linscott]

Hi all,

I'm using Environ to calculate delta SCF PBE electron affinities for small molecules as per Nattino et al 2019, and I have a few specific questions for those of you who are authors of that paper. (I would have sent this directly to the authors, but as some of you are here I thought it might be helpful for this conversation to be in the semi-public domain!)

1. It is important to be able to detect when a calculation for smaller values of epsilon_0 has "broken down" (and therefore should not be included in the polynomial fit). Is there a clever/correct way of doing this? Initially I was simply discounting calculations that failed to converge, but some cases* have convinced me to monitor the energies of the HOMO/LUMO, and discount calculations whenever the orbitals reorder relative to the large-epsilon_0-limit. I figure that the moment this happens, the extrapolation starts to stop making sense. I'm interested in hearing what you did and if you have any comments/suggestions!
2. The quartic fit to Delta E' gives some strange behaviour in many cases; often fitting the large-epsilon_0-limit comes at the expense of the much more relevant small-epsilon_0 cases. Did you explore alternative fitting schemes? (e.g. preferential weighting of small epsilon_0 data points?)
3. I'm encountering a handful of "too many bands not converged" errors for calculations within a cavity, when the same calculation without the cavity happily converges. Should these errors be resolved in the same way as with vanilla quantum espresso (that is, by increasing ecutwfc and ecutrho) or when Environ is involved is there something else I should try?
4. Finally, would you be happy to check the environ.in input file I've been using (content below)? It's simple enough and I'm pretty confident it's set up correctly for this method, but it would be nice to have this confirmed before I burn more CPU hours on these calculations!

* that is, cases where at small values of epsilon_0 the calculation converged but delta E' jumps relative to the calculations at large epsilon_0

Thanks so much!

Edward 

environ.in contents:

! For calculations determining EAs (Nattino et al 2019)

&ENVIRON

   verbose = 0

   environ_type = 'input'

   env_surface_tension = 0

   env_pressure = 0

   env_static_permittivity = <chosen value for epsilon_0>

/

&BOUNDARY

/

&ELECTROSTATIC

/

N.B. I am setting assume_isolated  = 'm-t' in the pw input file to remove periodicity.

-- 

Edward Linscott
Postdoctoral researcher
Theory and Simulation of Materials (THEOS)
EPFL

[Oliviero Andreussi]

Hi Edward,

Thanks for posting on the group, I do agree that it would be better to have these questions and answers shared openly, so that everybody can browse them, as the collection grows. 

Coming to your questions, some are a bit technical and, although I am one of the authors of that paper, I am not sure I can answer all of them with the required details. I will for sure try to forward them to Francesco Nattino, who did most of the calculations reported in the final paper. In the meantime, there are my comments.

1. I agree, sometimes the calculation converges, but the orbital switched order and the plain extrapolation would not work as well. I am not sure if there is any automatic way of doing this, I am pretty sure in the current era of machine learning algorithms there should be already something for this task, but we haven't gone that far yet, my understanding is that all the fitting and filtering we did was done manually.

2. It may be that a fit weighted more on the low epsilon values gives more accurate estimates, provided no band switching occurs. I struggled to find a reasonable way to estimate the accuracy of the extrapolation, but also in this case I am pretty sure that Baysian statistics could give some answers. We did tests comparing extrapolations with results for which QE converged and the estimates we reported in the paper are based on these. 

3. This is somehow not expected, in the sense that convergence should be a bit easier with the embedding than without. Thus, in general, I would expect that calculations that converge in vacuum also converge with the continuum environment. There are exceptions, a notable one being TiO2, that are particularly sensitive to the presence of the continuum. Most of these exceptions are due to something 'unphysical' happening between the electrons and the cavity, especially when using the SCCS cavity (i.e. the one defined on the electronic density). For TiO2 the problem seems to be that the electrons on the Ti atoms close to the interface are very sensitive to the local potential and change their density a lot when the continuum is nearby, making the SCF hard to converge. This was solved by pushing the dielectric a bit further away from the surface atoms, using a rigid interface (soft-spheres) or using the solvent-aware correction. The continuum getting too close to the ions can also cause problems with some acidic hydrogen atoms or with lithium atoms, which are particularly empty of electrons. I would think that the solvent-aware correction may improve in a few of these cases. The other solution could be to include solvent molecules in the system, near problematic atoms. Another reason a continuum calculation may be ill-defined is if the continuum interface is somehow too sharp in some regions. This could be solved similarly to what you say, by increasing ecutrho. Note that Environ is only defined on the density real-space grid, so changing ecutwfc does nothing to Environ. One last reason for which Environ may mess up the convergence is if the potential of Environ is added too late in the SCF step, when the electrons are already in a good state and the potential of Environ kicks them out too strongly. This can be solved by adding Environ earlier on during the simulation, thus increasing the value of environ_thr. This value is linked to the SCF accuracy reported by PW, i.e. Environ only computes its potentials when the SCF is below the SCF accuracy. 

4. Your input looks more than reasonable. MT should work fine with Environ. Note that you can also use the confining potential (env_confine) to achieve the same task, i.e. localizing the extra electron of anions on the system. 

Best,

Oliviero

[Mauro Sgroi]

Dear Oliviero,

regarding point (3) I'm observing the same problem of convergence. 

I'm concentrating on two systems: a doped graphene monolayer with an aluminum atom sustituting one graphitic carbon and the same system with a Li-polysulfide bonded to the aluminum atom. The Al atom is out of plane with respect to the graphene layer.

The first system does not converge so I suspect that the polysulfide has the same effect of placing explicit solvent molecules near the Al atom.

Best regards,

Mauro Sgroi.

[Oliviero Andreussi]

I have not seen this problem with bare aluminum, but it may very well be the same as with lithium or acid hydrogens. One check could be to use the solvent-aware correction, to see if the problem disappears. If you need help with that, just let me know.

Best

Oliviero

[Mauro Sgroi]

Dear Oliviero,

my convergence problem was solved using SSCS.

I tried to use SCCS setting solvent_radius to 7 bohr but the calculation was not cenverging.

Is there something else that I can try if I want to use the SCCS?

Best regards,

Mauro.

[andr...@gmail.com]

Sorry for missing to reply on your last question, I missed it. The solvent-aware model is activated by the solvent radius (7 Bohr seems a bit large, but that should not be too relevant). However, the intensity of the solvent-awareness is controlled by filling_threshold. The idea of the solvent-aware correction is that regions of continuum that are surrounded by the QM system (i.e. vacuum) should be treated as vacuum. To decide if a continuum region has too much vacuum around the filled fraction is computed within a sphere of fixed radius, proportional to the solvent radius. By lowering the filling_threshold you are smoothing out the regions of continuum that penetrate inside the QM system. In theory, reasonable values should be larger than 0.7, the default is 0.825 (from geometric considerations assuming spherical regions). For TiO2 anatase, which is a typical example that does not converge with plain SCCS, we saw convergence with filling_threshold of 0.65-0.7. Note that these values start to be in the unphysical region, so this is a good test to check if the problem is really related to the closeness of the interface. The 'physical' solution in this case is to add a solvent molecule where the problem occurs. 

I am not sure how clear the above explanation is, please let me know if you need more details.

Best,

Oliviero

[Mauro Sgroi]

Dear Oliviero,

thanks a lot for the explanation.

Best regards,

Mauro.

[Edward Linscott]

Dear Oliviero,

Thanks for your comprehensive response. As you suggested I switched over to using SSCS and I no longer get any "too many bands not converged" errors, so thanks for that suggestion!

However, there are still several calculations which fail to converge if (and only if) epsilon_0 > 1. Decreasing "tol" resolved some cases, but not all. I have noticed that in these failed runs, there is a negative rho that persists beyond the first few iterations. Doubling ecutrho did not diminish the negative rho component, and increasing environ_thr did not avoid negative rho, either. Is there anything else you suggest I try? (I can provide example failed calculations if that would help)

N.B. These failures are presenting in charged calculations where I have tot_charge = -1 and tot_magnetization = 1. Obviously this demands a PBC correction; I have tried Martyna-Tuckerman and parabolic corrections but neither converged (though parabolic got closer to reaching the desired scf accuracy)

Best wishes,
Edward

[andr...@gmail.com]

Hi Edward,

Thanks for the feedbacks, they are very useful. I am not sure what you mean for a negative rho, do you refer to the electronic density or the polarization density of Environ? If you could send the output files it would be great. It seems strange that the calculations converge in vacuum and not in solution for these systems, as the dielectric is meant to stabilize the extra electron on the anions, which instead should be unbound in vacuum. 

How big are your systems? Are these the ions of the G2-1 set? If the molecules are small, I would expect no problems with strange dielectric interfaces. For isolated systems, Martyna Tuckermann should be the safest way to go, provided the cell size is twice as big as the system size. 

Best,

Oliviero

[Edward Linscott]

Hi Oliviero,

The systems are small molecules, some similar to the G2-1 set but others are marginally larger -- e.g. I've attached an example run for ethylbenzene which does not converge.

I should have been more specific sorry: by negative rho I mean negative electronic density as reported by QE e.g.

     negative rho (up, down):  1.135E-03 5.259E-01

Having looked into it more I think I've managed to mostly* eliminate the problem of the persistent negative electronic density (for future reference I think removing poorly-chosen starting_magnetization keywords were what fixed it).

However, the calculations still stagnate at SCF accuracies of around ~ 1e-6, which is far higher than the levels to which I usually try to converge things (though I'm aware the QE default for conv_thr is 1e-6...). If you have the time to look at the attached output file and have any suggestions that would be amazing!

* I said "mostly" because in this ethylbenzene calculation you'll notice that at several steps a negative electronic density is detected and the SCF accuracy jumps upwards. However, this gets "fixed" immediately i.e. the subsequent SCF step does not exhibit a negative density, compared to what I was seeing before which was a persistent negative electronic density at most/every SCF step.

Best wishes,

Edward

P.S. my messages have not been sending, so if you get this multiple times sorry for the spam!