[QE-users] Calculating kinetic energy along particular direction for each electron orbital
Joshua Mann
jomann at physics.ucla.edu
Sat Jan 16 21:08:28 CET 2021
Hello,
Thanks for the code! So what I am trying to calculate is the kinetic energy
of each orbital (the kinetic portion of the variable et, eigs in wfck2r's
output), for each k-point and band index. So ekin has the shape ekin(nks,
nbnd, 6) and I do not perform any summation (except over pw's). The
modified version of e_kin is attached.
>From this I get the same exceedingly large energy results as before, after
removing the wg(ibnd,ik) term. 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.
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?
Thanks,
Joshua Mann
Department of Physics and Astronomy, University of California, Los Angeles
On Fri, Jan 15, 2021 at 5:04 AM Paolo Giannozzi <p.giannozzi at gmail.com>
wrote:
> 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
>
> _______________________________________________
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20210116/280d001d/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: e_kin.f90
Type: application/octet-stream
Size: 2431 bytes
Desc: not available
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20210116/280d001d/attachment.obj>
More information about the users
mailing list