<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">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>