[QE-developers] cgsolve_all and orthogonality constraints

[Paolo Giannozzi]

Sorry for replying so late: I opened an issue on GitLab (
https://gitlab.com/QEF/q-e/-/issues/272) in the vain hope that somebody was
going to have a look at it. I cannot believe that it is a bug: it would be
too serious to go unnoticed for so much time. I suspect that the results
are orthogonalized somewhere to the occupied Kohn-Sham manifold before
being used to compute anything.

Paolo


On Tue, Jan 12, 2021 at 1:58 PM Stoppels Harmen <harmen.stoppels at cscs.ch>
wrote:

> Hi all,
>
>
> I was trying to understand the linear response / perturbation theory code,
> in particular how the Sternheimer equation is solved. I believe there is
> some issue with the application of the preconditioner, and loss of
> orthogonality of the perturbation dpsi to the wave functions psi. Let me
> share my understanding of all this first:
>
>
> As I see it the Sternheimer equation is the linear system you have to
> solve when you do a single Newton iteration on a perturbed eigenproblem,
> with the unperturbed eigenproblem as the initial guess.
>
>
> When deriving the Sternheimer equation that way, you end up with the
> requirement that the correction to the unperturbed wave functions is
> orthogonal to the unperturbed wave functions, so <dpsi|S|psi> = O.
>
>
> And the the Sternheimer equation can then be formulated like this:
>
>
> (I - SQQ')(H - eS)(I - QQ'B)dpsi = (I - SQQ') * -residual, where Q is all
> psi's as a tall and skinny matrix, e is an eigenvalue corresponding to psi,
> and residual = dV * psi when the potential is perturbed, subject to dpsi
> perpendicular to psi in the S-inner product.
>
>
> I think in QE you work out this product of projectors on both sides, and
> then approximate it even (?) to obtain this H - eS + alpha * SQQ'S
> operator, is that correct?
>
>
> Now if you solve this problem using conjugate gradients, and B = I, and
> the initial guess for dpsi satisfies the orthogonality constraint, you
> could in principle drop those projectors, since the Krylov subspace is
> S-orthogonal to Q by design. That's a neat property to exploit. However,
> due to rounding errors, orthogonality can be lost, so adding the projectors
> is always a good idea. In particular CG might behave badly if not, since
> the operator is only positive definite when acting on vectors orthogonal to
> Q.
>
>
> However, the cgsolve_all code allows for a preconditioner P, which is
> typically chosen to be the diagonal of H - eS. Adding this preconditioner
> destroys the orthogonality, since there's no guarantee that P * x does not
> introduce components in the span of Q.
>
>
> I've added an example here to monitor the orthogonality of dpsi w.r.t. psi
> when running the PHonon/examples/example_01 (if I recall correctly):
> https://gist.github.com/haampie/92aa648b22c8a42d7b53e93fa59cf066. It
> outputs <dpsi|S|psi> as some kind of overlap matrix (if that is the right
> terminology). Without preconditioning you get values of order 1e-11, with
> preconditioner you get values of order 1e-4.
>
>
> Now the question is: is this intentional? Or is this a problem? The way to
> solve it would be to surround P with projectors on both sides too (I -
> SQQ')P(I - QQ'B), which can be simplified to (I - SQQ')P because the
> result of H - eS with projector can be assumed to be S-orthogonal to Q
> already. So in fact it would come down to applying P, and then calling
> `orthogonalize` on the output. Is that something that has been used before?
>
>
> Best regards,
>
>
> Harmen Stoppels
>
> CSCS
>
>
>
> _______________________________________________
> developers mailing list
> developers at lists.quantum-espresso.org
> https://lists.quantum-espresso.org/mailman/listinfo/developers
>


-- 
Paolo Giannozzi, Dip. Scienze Matematiche Informatiche e Fisiche,
Univ. Udine, via delle Scienze 206, 33100 Udine, Italy
Phone +39-0432-558216, fax +39-0432-558222
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20210716/c54e7582/attachment.html>