[QE-users] Calculating kinetic energy along particular direction for each electron orbital
Paolo Giannozzi
p.giannozzi at gmail.com
Fri Jan 15 14:03:41 CET 2021
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")
Paolo
On Fri, Jan 15, 2021 at 12:14 AM Joshua Mann <jomann at physics.ucla.edu>
wrote:
> 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
> _______________________________________________
> Quantum ESPRESSO is supported by MaX (www.max-centre.eu)
> users mailing list users at lists.quantum-espresso.org
> https://lists.quantum-espresso.org/mailman/listinfo/users
--
Paolo Giannozzi, Dip. Scienze Matematiche Informatiche e Fisiche,
Univ. Udine, via delle Scienze 206, 33100 Udine, Italy
Phone +39-0432-558216, fax +39-0432-558222
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20210115/bde00ab4/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: e_kin.f90
Type: text/x-fortran
Size: 2194 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20210115/bde00ab4/attachment.bin>
More information about the users
mailing list