[QE-developers] cgsolve_all and orthogonality constraints

[Stoppels Harmen]

Excuse me, a few typo's in the math expressions:

(I - SQQ')P(I - QQ'B) should be (I - SQQ')P(I - QQ'S)


and similarly B = I should read S = I.

________________________________
From: developers <developers-bounces at lists.quantum-espresso.org> on behalf of Stoppels Harmen <harmen.stoppels at cscs.ch>
Sent: Tuesday, January 12, 2021 1:57:48 PM
To: developers at lists.quantum-espresso.org
Subject: [QE-developers] cgsolve_all and orthogonality constraints


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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20210112/9edee240/attachment.html>