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 <[email protected]> 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 [email protected] > 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
! ! Copyright (C) 2019 Quantum ESPRESSO group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! !---------------------------------------------------------------------------- SUBROUTINE e_kin ( ) !---------------------------------------------------------------------------- ! ! ... Calculation of kinetic energy contribution to total energy ! ... Written by Paolo Giannozzi between Trieste and Palmanova ! ... To be called in electrons_scf, for instance after "sum_band" ! USE kinds, ONLY : DP USE mp, ONLY : mp_sum USE mp_bands, ONLY : intra_bgrp_comm USE mp_pools, ONLY : inter_pool_comm USE klist, ONLY : nks, ngk, wk USE buffers, ONLY : get_buffer USE control_flags, ONLY : gamma_only USE io_files, ONLY : iunwfc, nwordwfc USE noncollin_module, ONLY : noncolin USE wvfct, ONLY : g2kin, wg, nbnd, npwx USE wavefunctions, ONLY : evc ! IMPLICIT NONE ! REAL(dp) :: ekin INTEGER :: npw, ik, ig, ibnd ! ekin = 0.0_dp ! k_loop: DO ik = 1, nks ! npw = ngk(ik) IF ( nks > 1 ) & CALL get_buffer ( evc, nwordwfc, iunwfc, ik ) CALL g2_kin (ik) ! !$omp parallel do default(shared), private(ibnd,ig) reduction(+:ekin) band_loop: DO ibnd = 1, nbnd ! pw_loop: DO ig = 1, npw ekin = ekin + wg(ibnd,ik) * g2kin(ig) * & ( REAL(evc(ig,ibnd),KIND=dp)**2 + AIMAG(evc(ig,ibnd))**2 ) END DO pw_loop IF ( noncolin ) THEN pw_loop2: DO ig = 1, npw ekin = ekin + wg(ibnd,ik) * g2kin(ig) * & ( REAL(evc(ig+npwx,ibnd),KIND=dp)**2 + & AIMAG(evc(ig+npwx,ibnd))**2 ) END DO pw_loop2 END IF ! END DO band_loop !$omp end parallel do ! END DO k_loop IF ( gamma_only ) ekin = ekin*2.0_dp ! CALL mp_sum (ekin, intra_bgrp_comm) CALL mp_sum (ekin, inter_pool_comm) ! END SUBROUTINE e_kin
_______________________________________________ Quantum ESPRESSO is supported by MaX (www.max-centre.eu) users mailing list [email protected] https://lists.quantum-espresso.org/mailman/listinfo/users
