This post is about Quantum ESPRESSO (QE), an open source Density Functional Theory (DFT) code implemented in the framework of pseudopotentials and plane waves.
If you don't have a clue of what I am talking about, you may want to stop reading. Even though I will try to keep it as simple as I can, it might still require some background.
There is much to grok about QE but I will just tell you about one error message this time:
from mix_rho : error # 1
negative dr2
The error is obviously from mix_rho.f90 routine in PW. You would most likely encounter it while doing a self-consistent-field (scf) step. Here, dr2 is a measure of the error one makes in total energy in the scf cycle.
As you all know Kohn-Sham equations are solved self-consistently here: Start with a guess wavefunction, calculate the density (rho_in), write potential as functional of density, insert this potential into KS hamiltonian, solve the hamiltonian to get the new density(rho_out).
In an ideal world, you would keep doing this until self-consistency is reached, where rho_in = rho_out, exactly so. In the numerical world we do this a bit differently due to obvious reasons. So how do we do?
The convergence criterion
1. A naive way of checking if things are converged enough would be to compare two densities, rho_in and rho_out. Density is defined on a grid. So one may choose the maximum change in density as the convergence criterion (like max_change <1d-14) Or a vector that measures the total change between these two arrays can be used to check for the convergence.Why it would not be my choice is because, density changing a little does not necessarily mean that it is close to convergence. (Think why) Density changing a bit also does not mean that energy is converging. (The change might be small but its effect combined with potential might result in considerably big errors in energy)
2.What pw.x does: As far as I can understand from the code, pw.x calculates an upperbound for the change in energy due to change in density, it is proportional to the hartree potential energy density due to density rho_out ; assuming that the change due to exchange and correlation part is negligible (When I have a good latex widget I can also work out the details of this, it might get more clear.)
The scf cycle is assumed converged when this quantity (dr2) is less than a user provided value (conv_thr=1d-6 as default)
Now lets go back to the error message: negative dr2. As you probably guessed from the above description this term should be positive except for pathological cases. However I have realized that one can bump into this problem rather often in the case of PAW pseudopotentials. And the reason for that is as follows:
The accuracy of the projector augmented wave (PAW) method relies on good projectors in the "core" region. It is only with accurate projectors one can hope to reconstruct the all electron density. In the case of PAW pseudopotentials, dr2 term is corrected for the core region using the all electron density. How? Well, one subtracts the contribution of pseudo part in dr2 and adds the all electron one. In principle it should really work like superposition principle: You have the hartree energy due to rho_out in all space, then you calculate the contribution of this pseudo density in the core region and subtract it , and then add the contribution due to all electron density in the core region. When dr2 is negative that means you are subtracting too much and adding too few: The projectors are not good enough and/or the pseudo density is more localized than the all electron one, which is generally not the case except norm conserving pseudopotentials.
"Good enough" is a relative thing though. Therefore I am not suggesting you to throw away the pseudopotential just yet. Whenever things are going wrong in the radial grid in the core region it is a good idea to check plane wave expansion first: the ecuts. Also, it is possible that dr2 is negative but too small. Or if this happens in the very beginning of an scf run, it might as well be because projectors are not good enough for that arbitrary guess wavefunction, especially if one is very far from the generation configuration of the pseudopotential. What I would suggest you instead would be to check the pseudopotential thoroughly, do all the tests you can. If you believe it is "ghost-free" and actually not that bad of a pseudo, then go back to your failing calculation. This time run it with the "uncorrected" dr2; simply by commenting the following line in scf_mod.f90
IF (okpaw) rho_ddot = rho_ddot + paw_ddot(rho1%bec, rho2%bec)
inside function rho_ddot. Now remake the pw.x executable and run it again.
If your calculation goes alright without any further error and converges to a reasonable place, well you might have just saved yourself the hassle of generating a new pseudopotential. Beware though, I have not tested this very well but I've seen on pw_forum that some people are actually trying this approach (it has been discussed in today's posts).
Well, I hope that was helpful for some of you out there. I would be glad if you can let me know of the results. And please do not hesitate to correct me if you think I have said some BS.
Good luck!
No comments:
Post a Comment