[QE-users] Calculating kinetic energy along particular direction for each electron orbital

[Joshua Mann]

Hello all,

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.

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.

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.


I've ignored npol for now since it is just 1 in my case.
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).

The input wavefunction psi is just the output evc from read_collected_wfc,
but I am not sure if this is correct.

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.

I've found a very similar post (
https://lists.quantum-espresso.org/pipermail/users/2013-November/028007.html
), but it quickly moved to a private conversation.

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


Thanks,

Joshua Mann
Department of Physics and Astronomy, University of California, Los Angeles


Source code:


g2_kin_d.f90

  SUBROUTINE g2_kin_d( ik )

!----------------------------------------------------------------------------
  !! Calculation of kinetic energy along each dimension. Used just like
g2_kin
  !! Results placed in wvfct:g2kind
  !
  USE kinds,                ONLY : DP
  USE cell_base,            ONLY : tpiba2
  USE klist,                ONLY : xk, ngk, igk_k
  USE gvect,                ONLY : g
  USE gvecw,                ONLY : ecfixed, qcutz, q2sigma
  USE wvfct,                ONLY : g2kind
  !
  IMPLICIT NONE
  !
  INTEGER, INTENT(IN) :: ik
  !
  INTEGER :: npw
  !
  npw = ngk(ik)
  g2kind(1:npw,1) = ( xk(1,ik) + g(1,igk_k(1:npw,ik)) )**2 * tpiba2
  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
  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
  g2kind(1:npw,4) = ( xk(2,ik) + g(2,igk_k(1:npw,ik)) )**2 * tpiba2
  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
  g2kind(1:npw,6) = ( xk(3,ik) + g(3,igk_k(1:npw,ik)) )**2 * tpiba2
  !
  RETURN
  !
END SUBROUTINE g2_kin_d


psi_t_psi.f90

SUBROUTINE psi_t_psi( ik, lda, n, m, psi, psitpsi )
!----------------------------------------------------------------------------
  !! This routine computes the kinetic energies (<Tij>) of input wfcs
  !

USE kinds, ONLY : DP
USE wvfct, ONLY : g2kind

!
    IMPLICIT NONE
    !
INTEGER, INTENT(IN) :: ik
!! Index of k-point of wfc
    INTEGER, INTENT(IN) :: lda
    !! leading dimension of arrays psi, spsi, hpsi
!! (Josh: Maximum possible number of PWs)
    INTEGER, INTENT(IN) :: n
    !! true dimension of psi, spsi, hpsi
!! (Josh: Number of PWs here)
    INTEGER, INTENT(IN) :: m
    !! number of states psi
!! (Josh: num states :) )
    COMPLEX(DP), INTENT(IN) :: psi(lda,m)
    !! the wavefunction
    REAL(DP), INTENT(OUT) :: psitpsi(m, 6)
    !! The kinetic energy (<Txx>, <Txy>, <Txz>, <Tyy>, <Tyz>, <Tzz>) of
each wfc

REAL(DP) :: tn
!! Running total norm
INTEGER :: i, j

CALL g2_kin_d( ik )

psitpsi = 0.0

DO i = 1, m
tn = 0
DO j = 1, n
psitpsi(i, :) = psitpsi(i, :) + (g2kind(j, :) * (abs(psi(j, i))**2))
tn = tn + abs(psi(j, i))**2
ENDDO
psitpsi(i, :) = psitpsi(i, :) / tn
ENDDO

RETURN

END SUBROUTINE psi_t_psi


Code placed in wfck2r.f90:

  ALLOCATE( psitpsi(last_band-first_band+1,6) )

  IF (loctave .and. ionode) THEN
write(6,*) 'writing directional kinetic energies...'
    write(iuwfcr+1,'(/)')
    write(iuwfcr+1,'("# name: ",A,/,"# type: matrix")') 'ti'
    write(iuwfcr+1,'("# ndims: 3")')
    write(iuwfcr+1,'(5I10)') 6, last_band-first_band+1, last_k-first_k+1
DO ik = first_k,last_k
npw = ngk(ik)
CALL read_collected_wfc ( restart_dir(), ik, evc )
CALL psi_t_psi( ik, npwx, npw, last_band-first_band+1, evc, psitpsi )
!write(6,*) ((psitpsi(j,i)*rytoev, i=1,6), j=1,last_band-first_band+1)
write(iuwfcr+1,'(E20.12)') ((psitpsi(j,i)*rytoev, i=1,6),
j=1,last_band-first_band+1)
ENDDO
ENDIF
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20210114/3993e2c5/attachment.html>