<div dir="ltr"><div>Sorry for replying so late: I opened an issue on GitLab (<a href="https://gitlab.com/QEF/q-e/-/issues/272">https://gitlab.com/QEF/q-e/-/issues/272</a>) 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. <br></div><div><br></div><div>Paolo<br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Jan 12, 2021 at 1:58 PM Stoppels  Harmen <<a href="mailto:harmen.stoppels@cscs.ch">harmen.stoppels@cscs.ch</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">




<div dir="ltr">
<div id="gmail-m_4886213212850651419divtagdefaultwrapper" style="font-size:12pt;color:rgb(0,0,0);font-family:Calibri,Helvetica,sans-serif" dir="ltr">
<p>Hi all,</p>
<p><br>
</p>
<p>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:<br>
</p>
<p><br>
</p>
<p>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.</p>
<p><br>
</p>
<p>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.</p>
<p><br>
</p>
<p>And the the Sternheimer equation can then be formulated like this: <br>
</p>
<p><br>
</p>
<p>(I - SQQ')(H - eS)(I - QQ'B)dpsi = <span>(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</span> when the potential is perturbed,
<span>subject to dpsi perpendicular to psi in the S-inner product.</span></p>
<p><span><br>
</span></p>
<p><span><span>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</span> operator, is that correct?<br>
</span></p>
<p><span><br>
</span></p>
<p><span>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, si</span>nce 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.</p>
<p><br>
</p>
<p>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.</p>
<p><br>
</p>
<p>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):
<a href="https://gist.github.com/haampie/92aa648b22c8a42d7b53e93fa59cf066" id="gmail-m_4886213212850651419LPlnk632815" target="_blank">
https://gist.github.com/haampie/92aa648b22c8a42d7b53e93fa59cf066</a>. It outputs <dpsi|S|psi> as some kind of overlap matrix (if that is the right terminology). Without preconditioning you get values of order 1<span>e-11</span>, with preconditioner you get
 values of order 1e-4.</p>
<p><br>
</p>
<p>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
<span>(I - SQQ')P(I - QQ'B)</span>, which can be simplified to <span>(I - SQQ')P</span> because the result of
<span>H - eS</span> 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?</p>
<p><br>
</p>
<p>Best regards,</p>
<p><br>
</p>
<p>Harmen Stoppels</p>
<p>CSCS<br>
</p>
<p><br>
</p>
<p><br>
</p>
</div>
</div>

_______________________________________________<br>
developers mailing list<br>
<a href="mailto:developers@lists.quantum-espresso.org" target="_blank">developers@lists.quantum-espresso.org</a><br>
<a href="https://lists.quantum-espresso.org/mailman/listinfo/developers" rel="noreferrer" target="_blank">https://lists.quantum-espresso.org/mailman/listinfo/developers</a><br>
</blockquote></div><br clear="all"><br>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>Paolo Giannozzi, Dip. Scienze Matematiche Informatiche e Fisiche,<br>Univ. Udine, via delle Scienze 206, 33100 Udine, Italy<br>Phone +39-0432-558216, fax +39-0432-558222<br><br></div></div></div></div></div>