<div dir="ltr"><div dir="ltr">On Sat, Jan 16, 2021 at 9:09 PM Joshua Mann <<a href="mailto:jomann@physics.ucla.edu">jomann@physics.ucla.edu</a>> wrote:</div><div dir="ltr"><br></div><div class="gmail_quote"><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">From this I get the same exceedingly large energy results as before
, after removing the wg(ibnd,ik) term.
</div></blockquote><div><br></div><div>what makes you think that those values are not correct? Without the weights (band occupation times a symmetry weight) what you compute is just <\psi|T|\psi>, for the plane-wave part only (the augmentation terms are missing)<br></div><div><br></div><div>Paolo<br> </div><div><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>With wg the kinetic energy is too low, and my understanding is that this is used to weight each wavefunction when summing over the k-points and bands.</div><div><br></div><div>Is there some other consideration I am missing when I only want the kinetic energy of each orbital, and not the total contribution of a band or k-point?</div><div><br></div><div><br></div><div>Thanks,</div><div><br></div><div>Joshua Mann</div><div>Department of Physics and Astronomy, University of California, Los Angeles</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 15, 2021 at 5:04 AM Paolo Giannozzi <<a href="mailto:p.giannozzi@gmail.com" target="_blank">p.giannozzi@gmail.com</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>The attached piece of code computes (correctly the last time I tried) the kinetic energy. You may easily extend it to compute pieces of the kinetic energy for each band: use your routine "g2_kin_d" (pass variable g2kind as argument), replace ekin with ekin(nbnd,6), remove the sum over k-point (also the last mp_sum; move the previous one inside the loop over k-points). Don't modify anything else in QE: just add a call to this routine (preferably at the end, subroutine "punch")<br></div><div><br></div><div>Paolo<br> </div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 15, 2021 at 12:14 AM Joshua Mann <<a href="mailto:jomann@physics.ucla.edu" target="_blank">jomann@physics.ucla.edu</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">Hello all,<div><br></div><div>I wish to calculate the kinetic energy along a particular direction for each orbital (i.e., along n I want 1/2m <p_n^2>). To do this I need to calculate 1/2m <p_i p_j> for each orbital so I can later pick any arbitrary direction.</div><div><br></div><div>To do this I wrote a code to calculate an equivalent quantity to g2kin, g2kind, which has an extra dimension length 6 (to capture all combinations i,j). This code is at the end of this email.</div><div><br></div><div>And I wrote a code to calculate the expectation of this energy, also at the end of this email. I've modified a few other files to get it to compile.</div><div><br></div><div><br></div><div>I've ignored npol for now since it is just 1 in my case.</div><div>Performing this, and modifying wfck2r.f90 to output the results, I find that the kinetic energy (Txx+Tyy+Tzz) in bulk Tungsten can be way too large, 50-100 eV. Adding this kinetic energy to the potential energy <V> as calculated by unkr from wfck2r.x's output and the total potential from pp.x does not retrieve the total orbital energies eigs (from wfck2r).</div><div><br></div><div>The input wavefunction psi is just the output evc from read_collected_wfc, but I am not sure if this is correct. </div><div><br></div><div>Some oddities I've found are that tn, the running norm total for each wfc, ranges from 1 to 1.8 (I am using a PAW PP). I divide by this norm in case the wfc's are not necessarily normalized, but again I am not sure if this is correct.</div><div><br></div><div>I've found a very similar post ( <a href="https://lists.quantum-espresso.org/pipermail/users/2013-November/028007.html" target="_blank">https://lists.quantum-espresso.org/pipermail/users/2013-November/028007.html</a> ), but it quickly moved to a private conversation.</div><div><br></div><div>Any suggestions? I am quite new to Fortran and I've only gotten into QE's source code starting this week, so there may be several issues here...</div><div><br></div><div><br></div><div>Thanks,</div><div><br></div><div>Joshua Mann</div><div>Department of Physics and Astronomy, University of California, Los Angeles</div><div><br></div><div><br></div><div>Source code:</div><div><br></div><div><br></div><div>g2_kin_d.f90</div><div><br></div><div> SUBROUTINE g2_kin_d( ik )<br> !----------------------------------------------------------------------------<br> !! Calculation of kinetic energy along each dimension. Used just like g2_kin<br> !! Results placed in wvfct:g2kind<br> !<br> USE kinds, ONLY : DP<br> USE cell_base, ONLY : tpiba2 <br> USE klist, ONLY : xk, ngk, igk_k<br> USE gvect, ONLY : g<br> USE gvecw, ONLY : ecfixed, qcutz, q2sigma<br> USE wvfct, ONLY : g2kind<br> !<br> IMPLICIT NONE<br> !<br> INTEGER, INTENT(IN) :: ik<br> !<br> INTEGER :: npw<br> !<br> npw = ngk(ik)<br> g2kind(1:npw,1) = ( xk(1,ik) + g(1,igk_k(1:npw,ik)) )**2 * tpiba2<br> g2kind(1:npw,2) = ( xk(1,ik) + g(1,igk_k(1:npw,ik)) ) * ( xk(2,ik) + g(2,igk_k(1:npw,ik)) ) * tpiba2<br> g2kind(1:npw,3) = ( xk(1,ik) + g(1,igk_k(1:npw,ik)) ) * ( xk(3,ik) + g(3,igk_k(1:npw,ik)) ) * tpiba2<br> g2kind(1:npw,4) = ( xk(2,ik) + g(2,igk_k(1:npw,ik)) )**2 * tpiba2<br> g2kind(1:npw,5) = ( xk(2,ik) + g(2,igk_k(1:npw,ik)) ) * ( xk(3,ik) + g(3,igk_k(1:npw,ik)) ) * tpiba2<br> g2kind(1:npw,6) = ( xk(3,ik) + g(3,igk_k(1:npw,ik)) )**2 * tpiba2<br> !<br> RETURN<br> !<br>END SUBROUTINE g2_kin_d </div><div><br></div><div><br></div><div>psi_t_psi.f90 </div><div><br></div><div>SUBROUTINE psi_t_psi( ik, lda, n, m, psi, psitpsi )</div><div>!----------------------------------------------------------------------------<br> !! This routine computes the kinetic energies (<Tij>) of input wfcs<br> !<br> <br> USE kinds, ONLY : DP<br> USE wvfct, ONLY : g2kind<br> <br> !<br> IMPLICIT NONE<br> !<br> INTEGER, INTENT(IN) :: ik<br> !! Index of k-point of wfc<br> INTEGER, INTENT(IN) :: lda<br> !! leading dimension of arrays psi, spsi, hpsi<br> !! (Josh: Maximum possible number of PWs)<br> INTEGER, INTENT(IN) :: n<br> !! true dimension of psi, spsi, hpsi<br> !! (Josh: Number of PWs here)<br> INTEGER, INTENT(IN) :: m<br> !! number of states psi<br> !! (Josh: num states :) )<br> COMPLEX(DP), INTENT(IN) :: psi(lda,m) <br> !! the wavefunction<br> REAL(DP), INTENT(OUT) :: psitpsi(m, 6)<br> !! The kinetic energy (<Txx>, <Txy>, <Txz>, <Tyy>, <Tyz>, <Tzz>) of each wfc<br> <br> REAL(DP) :: tn<br> !! Running total norm<br> INTEGER :: i, j<br> <br> CALL g2_kin_d( ik )<br><br> psitpsi = 0.0<br><br> DO i = 1, m<br> tn = 0<br> DO j = 1, n<br> psitpsi(i, :) = psitpsi(i, :) + (g2kind(j, :) * (abs(psi(j, i))**2))<br> tn = tn + abs(psi(j, i))**2<br> ENDDO<br> psitpsi(i, :) = psitpsi(i, :) / tn<br> ENDDO<br> <br> RETURN<br><br>END SUBROUTINE psi_t_psi <br></div><div><br></div><div><br></div><div>Code placed in wfck2r.f90:</div><div><br></div><div> ALLOCATE( psitpsi(last_band-first_band+1,6) )<br><br> IF (loctave .and. ionode) THEN<br> write(6,*) 'writing directional kinetic energies...'<br> write(iuwfcr+1,'(/)')<br> write(iuwfcr+1,'("# name: ",A,/,"# type: matrix")') 'ti'<br> write(iuwfcr+1,'("# ndims: 3")')<br> write(iuwfcr+1,'(5I10)') 6, last_band-first_band+1, last_k-first_k+1<br> DO ik = first_k,last_k<br> npw = ngk(ik)<br> CALL read_collected_wfc ( restart_dir(), ik, evc )<br> CALL psi_t_psi( ik, npwx, npw, last_band-first_band+1, evc, psitpsi )<br> !write(6,*) ((psitpsi(j,i)*rytoev, i=1,6), j=1,last_band-first_band+1)<br> write(iuwfcr+1,'(E20.12)') ((psitpsi(j,i)*rytoev, i=1,6), j=1,last_band-first_band+1)<br> ENDDO<br></div><div>ENDIF</div></div>
_______________________________________________<br>
Quantum ESPRESSO is supported by MaX (<a href="http://www.max-centre.eu" rel="noreferrer" target="_blank">www.max-centre.eu</a>)<br>
users mailing list <a href="mailto:users@lists.quantum-espresso.org" target="_blank">users@lists.quantum-espresso.org</a><br>
<a href="https://lists.quantum-espresso.org/mailman/listinfo/users" rel="noreferrer" target="_blank">https://lists.quantum-espresso.org/mailman/listinfo/users</a></blockquote></div><br clear="all"><br>-- <br><div dir="ltr"><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>
_______________________________________________<br>
Quantum ESPRESSO is supported by MaX (<a href="http://www.max-centre.eu" rel="noreferrer" target="_blank">www.max-centre.eu</a>)<br>
users mailing list <a href="mailto:users@lists.quantum-espresso.org" target="_blank">users@lists.quantum-espresso.org</a><br>
<a href="https://lists.quantum-espresso.org/mailman/listinfo/users" rel="noreferrer" target="_blank">https://lists.quantum-espresso.org/mailman/listinfo/users</a></blockquote></div>
_______________________________________________<br>
Quantum ESPRESSO is supported by MaX (<a href="http://www.max-centre.eu" rel="noreferrer" target="_blank">www.max-centre.eu</a>)<br>
users mailing list <a href="mailto:users@lists.quantum-espresso.org" target="_blank">users@lists.quantum-espresso.org</a><br>
<a href="https://lists.quantum-espresso.org/mailman/listinfo/users" rel="noreferrer" target="_blank">https://lists.quantum-espresso.org/mailman/listinfo/users</a></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></div>