[Pw_forum] Electrons kinetic energy

[Pietro Bonfa']

Dear Prof. Giannozzi,

thanks for your answer.

I believe that this is exactly what you suggested

   DO ik = 1, nk
       DO ibnd = 1, nbnd
           DO j = 1, npw
               IF(gamma_only)THEN !.and.j>1)then
                   ek = ek +  2*conjg(evc(j,ibnd)) * evc(j,ibnd) * &
                               g2kin(j) * wg(ibnd,ikk)
               ELSE
                   ek = ek +  conjg(evc(j,ibnd)) * evc(j,ibnd) * &
                               g2kin(j) * wg(ibnd,ikk)
               ENDIF
           ENDDO
        ENDDO
   ENDDO

I also started from the local kinetic energy density just to
double-check the code since I'm not familiar with the QE structure and
variables.

I encountered problems obtaining the total kinetic energy by summing
over the real grid the kinetic energy density.
I found that the problem is solved if I change the from 'Smooth' to
'Dense' the fft grid in the elf.f90 code (line 121) and I was wondering
if I did something wrong or not.

Anyway, it was almost only a curiosity. I'm happy with the code reported
above.

Thanks again,
Pietro Bonfa'



On 10/17/13 18:24, Paolo Giannozzi wrote:
> What do you want exactly to calculate? The "average electron
> kinetic energy" for norm-conserving pseudopotentials is a 
> rather simple object. In atomic (Ry) units:
>   \sum_k \sum_i f_{i,k} <\psi_{i,k}|\nabla^2|\psi_{i,k}> = 
>   \sum_k \sum_i f_{i,k} \sum_G |\psi_i(k+G)|^2 (k+G)^2
> where f_{i,k} is the occupancy of Kohn-Sham state \psi_{i,k}
> whose plane-wave components are \psi_i(k+G). YOu do not need
> the "local kinetic energy density": 
>   \sum_k \sum_i f_{i,k} |\nabla psi_{k,i}(r)|**2
> P.
> 
> On Wed, 2013-10-16 at 17:04 +0200, Pietro Bonfa' wrote:
>> Dear All,
>>
>> following the recipes given in this thread:
>>
>> http://qe-forge.org/pipermail/pw_forum/2009-May/087347.html
>>
>> I added two routines to the pp.x code to calculate the average
>> electronic kinetic energy. I started from both the elf.f90 subroutine in
>> PP code and from the pw2casino implementation.
>>
>> It's almost only a matter of copy/paste but there is a thing I cannot
>> understand.
>>
>> Let "kkin" be the variable with the kinetic energy density, following
>> the elf.f90 code I have:
>>
>>      0    ! symmetrize the local kinetic energy if needed
>>      1	  IF ( .not. gamma_only) THEN
>>      2	     !
>>      3	     CALL sym_rho_init ( gamma_only )
>>      4	     !
>>      5	     aux(:) =  cmplx ( kkin (:), 0.0_dp, kind=dp)
>>      6	     CALL fwfft ('Dense', aux, dfftp)
>>      7	     ALLOCATE (aux2(ngm))
>>      8	     aux2(:) = aux(nl(:))
>>      9	     !
>>     10	     ! aux2 contains the local kinetic energy in G-space to be
>> symmetrized
>>     11	     !
>>     12	     CALL sym_rho ( 1, aux2 )
>>     13	     !
>>     14	     aux(:) = (0.0_dp, 0.0_dp)
>>     15	     aux(nl(:)) = aux2(:)
>>     16	     DEALLOCATE (aux2)
>>     17	
>>     18	     CALL invfft ('Dense', aux, dfftp)
>>     19	     kkin (:) = dble(aux(:))
>>     20	     !
>>     21	  ENDIF
>>     22	  T = sum( kkin(:) )*omega / ( dfftp%nr1*dfftp%nr2*dfftp%nr3 )
>>
>> To get the correct kinetic energy (T), with respect to elf.f90, I had to
>> change
>> CALL fwfft ('Smooth', aux, dfftp)
>> to
>> CALL fwfft ('Dense', aux, dfftp)
>> on line 6 because otherwise I couldn't perform the sum correctly (even
>> if using dffts%nr1 nr2 nr3).
>>
>> With this trick the elf-derived and the pw2casino-derived routines give
>> similar results in all the cases I used for testing.
>> Nonetheless, I'm completely unfamiliar with the QE implementation, and
>> probably there is no need to change Smooth to Dense.
>> Moreover, I don't really know if performing sym_rho on the Dense grid is
>> fine or not.
>>
>> Can someone please help me figure out my mistake(s)?
>>
>> Thanks,
>> Pietro
>>
>> ps: the subroutines are attached.
>>
>>
>> _______________________________________________
>> Pw_forum mailing list
>> Pw_forum at pwscf.org
>> http://pwscf.org/mailman/listinfo/pw_forum
> 


-- 
Pietro Bonfa' - PhD student
Dipartimento di Fisica e Scienze della Terra "Macedonio Melloni"
Viale delle Scienze 7A
43124 Parma - Italy