[QE-users] Calculating kinetic energy along particular direction for each electron orbital
Paolo Giannozzi
p.giannozzi at gmail.com
Tue Jan 19 08:58:37 CET 2021
On Sat, Jan 16, 2021 at 9:09 PM Joshua Mann <jomann at physics.ucla.edu> wrote:
>From this I get the same exceedingly large energy results as before , after
> removing the wg(ibnd,ik) term.
>
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)
Paolo
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
>
> _______________________________________________
> 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/20210119/f96d1d15/attachment.html>
More information about the users
mailing list