[QE-users] Calculating kinetic energy along particular direction for each electron orbital
Joshua Mann
jomann at physics.ucla.edu
Tue Jan 19 22:22:58 CET 2021
Interpreting the output as the full kinetic energy (Txx+Tyy+Tzz) and
estimating the potential energy <V> (using unkr from wfck2r.x and full
potential from pp.x) results in energies much larger than and not
coinciding with the eigenvalues also from wfck2r.x.
However if there is more to the kinetic energy than that calculated from
the PW portion, then these results are likely correct and my interpretation
needs improvement.
I am ultimately trying to calculate the density of (occupied) states as a
function of only the kinetic energy in one direction + <V> (as opposed to
the total energy), which is why I can't sum over bands or k-points
beforehand.
How might I include the PAW augmentation terms? I'm very new to the
concept. Or is there some other way I can extract this DOS from QE?
Thanks,
Joshua Mann
Department of Physics and Astronomy, University of California, Los Angeles
On Mon, Jan 18, 2021 at 11:59 PM Paolo Giannozzi <p.giannozzi at gmail.com>
wrote:
> 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
>
> _______________________________________________
> 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/20210119/ea79e728/attachment.html>
More information about the users
mailing list