A: There may be "physical" reason or "numerical" reasons that explain why your calculation explodes. The latter are easier to solve, by tuning some of the parameters in environ.in, the former are a bit more tricky and I would leave them to a second analysis. Roughly speaking, in a calculation with Environ you need to include some extra terms to the Kohn-Sham potential of the system in vacuum. These terms depend on the electronic density itself, in particular due to the fact that the interface between the system and the environment is defined in terms of the electronic density. Moreover, to compute one of these terms (probably the most important, the electrostatic interaction with the dielectric continuum) Environ uses an iterative approach, so every time you need to add the environ potential in an SCF step, Environ performs an additional iterative procedure to compute it: let's call this the polarisation iteration, in contrast to the SCF iteration.
Having said this, there are two environ parameters which are crucial to help convergence of the total SCF calculation:
tolrhopol: this parameter controls the convergence of the polarisation iteration. If you decrease it, the extra term coming from the interaction with the dielectric will be computed more accurately, but it will require more time to be computed. On the other hand, if you don't compute it accurately enough, the SCF accuracy will probably not go down as fast as it should and you may end up with an SCF accuracy that oscillates above the selected threshold (so that pw.x spends lots of time without converging). Since a single polarisation step only costs as much as a few FFTs, while a single SCF step costs more, it is generally better to INCREASE THE ACCURACY OF THE POLARISATION CALCULATION (I.E. DECREASE TOLRHOPOL) until you have a steep convergence of the SCF.
environ_thr: since the extra terms depend in a delicate way on the electronic density (the shape of the environment depends on it), it is advisable not to add these terms at the very beginning of the SCF, which in pw.x usually initialise the electronic density in a random fashion. Thus, to enforce that the environment potentials are only added when the electronic density has a sensible shape, Environ only acts if the SCF accuracy is below a certain threshold, defined by environ_thr. On the other hand, you don't want to wait too long to add these extra terms, because it makes no sense to converge the density in vacuum and then add a very strong unexpected potential: the SCF accuracy will go banana and you need to redo all the SCF again. The net result is the one that you see in your output, you spend quite a lot SCF iterations to achieve the required SCF accuracy, then environ adds its additional potential and the SCF accuracy goes up again, with an oscillating down and up trend which may eventually explode. Thus, if you see poor convergence of the SCF the way to go is to INCREASE ENVIRON_THR, NOT TO DECREASE IT. As a side problem, SCF accuracy is an extensive quantity (it becomes larger for bigger systems), so that it is not possible to fix a default hardcoded value for this parameter. Nonetheless, as a rule of thumb, environ_thr should be set so as to skip only a few (three, four, five) SCF steps, so you may want to set environ_thr to the value of the SCF_accuracy at the third/fourth/fifth SCF step in vacuum. NOTE: if your SCF starts from a previous run in vacuum and the electronic density has an already well defined shape, you may want to start adding Environ potentials from the very beginning: you can do this by setting environ_restart = .true. inside the environ.in input file.
Apart from these two numerical parameters, a further source of errors may arise due to the use of pseudo potentials to handle core electrons. This problem is more common in transition metals or alogens, and it comes from the neglect of the core electrons in the charge density used to build the interface with the environment. In some cases there appears to be a hole in the charge density on the nuclei, which is then filled with the environment: this clearly creates big troubles when computing the environ contributions. To avoid this, you can try to set eps_mode = 'full' in the environ.in input file: this tells environ to use the electronic density plus a fictitious density of gaussians centred on each atom (to describe the core electrons).
Apart from these numerical tricks, you may want to simplify a bit the physics of your problem. In particular, since in most applications the most important effects of the environment are the electrostatic ones, you may want to switch off all the other non-electrostatic terms. You may do this by using environ_type = 'input' instead of the default ones (environ_type = 'water') and then by only specifying the dielectric constant of the environment env_static_permittivity = 78.3 for water at room temperature. Since by default all environ contributions are turned off, by doing this you are only including the dielectric continuum in your environ calculation. This is analogous to what is done in example02 and example03 in the Environ/examples directory.
Moreover, if you are simulating a 2D system, convergence of the polarisation potential may be made more difficult by the artificial finite electric field coming from periodic boundary conditions. To avoid this, it would be better to use a PBC correction, in particular for slabs in Environ it is implemented a simple parabolic correction, which corrects most of the artefacts and helps converging the polarisation iterations. You can set it up by setting assume_isolated = 'slabz' in the &SYSTEM namelist in your pw input file (NOT the environ.in input file). Note that you can use the same correction also for a calculation in vacuum, you still need to run pw.x --environ and set environ_type = 'vacuum' inside environ.in (similarly to what is done in example03).
Eventually, one reason of poor convergence of an Environ calculation may arise if your system's interface with the environment has a strange/irregular shape. This problem can be solved, but is a bit more cumbersome and I would only analyse it if the other things don't work.
Q. Why is the Harris-Foulkes estimate completely off with respect to the total energy when using Environ?
A. The Harris-Foulkes estimate does not take into account the extra terms coming from Environ, as it would cost more than necessary for a quantity that is only used to have an estimate. As a matter of fact, also the scf accuracy does not take into full account of the environment, but again it would need a secondary correction that would cost more than what its utility.
Q. Why is the total energy in an Environ calculation not the sum of the reported terms? What is the meaning of the environ contributions to the total energy?
A. The energy output in an Environ calculation may be 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 often confused. 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 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.