1!
2! Copyright (C) 2001-2013 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9
10!
11!----------------------------------------------------------------------------
12SUBROUTINE energies_xc( lda, n, m, psi, e_xc, e_h,ispin )
13  !----------------------------------------------------------------------------
14  !
15  ! computes the expectation values of the exchange and correlation potential
16  ! and of hartree potential
17  ! ... input:
18  ! ...    lda   leading dimension of arrays psi, spsi, hpsi
19  ! ...    n     true dimension of psi, spsi, hpsi
20  ! ...    m     number of states psi
21  ! ...    psi
22  !
23  ! ... output:
24  !       e_xc
25  !       e_h
26  USE kinds,    ONLY : DP
27  USE uspp,     ONLY : vkb, nkb
28  USE gvecs,  ONLY : doublegrid
29  USE gvect,                ONLY : ngm, gstart, g, gg, gcutm
30  USE cell_base,            ONLY :  alat, omega
31  USE lsda_mod,             ONLY : nspin
32  USE ldaU,     ONLY : lda_plus_u
33  USE lsda_mod, ONLY : current_spin
34  USE gvect,    ONLY : gstart
35  USE io_global, ONLY :stdout
36  USE scf,       ONLY : rho, vltot, vrs, rho_core,rhog_core, scf_type
37  USE constants,  ONLY :rytoev
38  USE io_files, ONLY : diropn
39  USE mp, ONLY : mp_sum, mp_barrier
40  USE mp_world, ONLY : world_comm
41  USE control_flags,        ONLY : gamma_only
42  USE funct,            ONLY : dft_is_meta
43  USE fft_base,             ONLY : dfftp, dffts
44  USE fft_interfaces,       ONLY : fwfft, invfft, fft_interpolate
45
46  USE exx,      ONLY : vexx !Suriano
47  USE funct,    ONLY : exx_is_active,dft_is_hybrid
48  USE klist, ONLY : igk_k
49
50  !
51  IMPLICIT NONE
52  INTEGER, EXTERNAL :: find_free_unit
53  !
54  ! ... input/output arguments
55  !
56  INTEGER          :: lda, n, m
57  COMPLEX(DP) :: psi(lda,m)
58  REAL(kind=DP) :: e_xc(m), e_h(m)
59  INTEGER, INTENT(in) :: ispin !spin 1,2
60
61  REAL(kind=DP), ALLOCATABLE :: vr(:,:)
62  !
63  !
64  CALL start_clock( 'h_psi' )
65  allocate(vr(dfftp%nnr,nspin))
66  !
67  IF ( gamma_only ) THEN
68     !
69     CALL energies_xc_gamma()
70     !
71  ELSE
72     !
73     CALL energies_xc_k( )
74     !
75  END IF
76
77  !
78  CALL stop_clock( 'h_psi' )
79  deallocate(vr)
80  !
81  RETURN
82  !
83  CONTAINS
84     !
85     !-----------------------------------------------------------------------
86     SUBROUTINE energies_xc_k( )
87       !-----------------------------------------------------------------------
88       !
89       ! ... k-points version
90       !
91       USE wavefunctions, ONLY : psic
92       USE becmod,  ONLY : becp
93       !
94       IMPLICIT NONE
95       !
96
97       INTEGER :: ibnd, j,is, ig
98       REAL(dp) :: etxc,vtxc
99       REAL(kind=DP) :: ehart, charge
100       ! counters
101       !
102       !
103       !
104       ! ... Here we apply the kinetic energy (k+G)^2 psi
105       !
106       !
107       !
108       ! ... Here we add the Hubbard potential times psi
109       !
110       !
111       ! ... the local potential V_Loc psi. First the psi in real space
112!set exchange and correlation potential
113          if(.not.allocated(psic)) write(stdout,*) 'psic not allocated'
114      !
115       if (dft_is_meta()) then
116!         call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r )
117      else
118         CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vr )
119      endif
120      !
121       do is=1,nspin
122          vrs(:,is)=vr(:,is)
123          if(doublegrid) call fft_interpolate(dfftp, vrs(:,is),dffts,vrs(:,is)) ! interpolate from dense to smooth
124       enddo
125       !
126       DO ibnd = 1, m
127          !
128          CALL start_clock( 'firstfft' )
129          !
130          psic(1:dffts%nnr) = ( 0.D0, 0.D0 )
131          !
132          psic(dffts%nl(igk_k(1:n,1))) = psi(1:n,ibnd)
133          !
134          CALL invfft ('Wave', psic, dffts)
135
136          !
137          CALL stop_clock( 'firstfft' )
138          !
139          ! ... product with the potential vrs = (vltot+vr) on the smooth grid
140          !
141          psic(1:dffts%nnr) = psic(1:dffts%nnr) * vrs(1:dffts%nnr,1)
142          !
143          ! ... back to reciprocal space
144          !
145          CALL start_clock( 'secondfft' )
146          !
147          CALL fwfft ('Wave', psic, dffts)
148          !
149          ! ... addition to the total product
150          !
151          e_xc(ibnd)=0.d0
152          do ig=1,n
153             e_xc(ibnd)=e_xc(ibnd)+real(conjg(psi(ig,ibnd))*psic(dffts%nl(igk_k(ig,1))))
154          enddo
155          call mp_sum(e_xc(ibnd),world_comm)
156          write(stdout,*) 'energies_xc :', ibnd, e_xc(ibnd)*rytoev
157!
158          CALL stop_clock( 'secondfft' )
159          !
160       END DO
161       vr(:,:)=0.d0
162       call  v_h(rho%of_g , ehart, charge, vr )
163       do is=1,nspin
164          vrs(:,is)=vr(:,is)
165          if(doublegrid) call fft_interpolate(dfftp, vrs(:,is), dffts, vrs(:,is)) ! interpolate from dense to smooth
166       enddo
167
168
169       DO ibnd = 1, m
170
171          CALL start_clock( 'firstfft' )
172          psic(1:dffts%nnr) = ( 0.D0, 0.D0 )
173          psic(dffts%nl(igk_k(1:n,1))) = psi(1:n,ibnd)
174
175          CALL invfft ('Wave', psic, dffts)
176
177
178          CALL stop_clock( 'firstfft' )
179
180
181          psic(1:dffts%nnr) = psic(1:dffts%nnr) * vrs(1:dffts%nnr,1)
182
183          CALL start_clock( 'secondfft' )
184          CALL fwfft ('Wave', psic, dffts)
185          e_h(ibnd)=0.d0
186          do ig=1,n
187             e_h(ibnd)=e_h(ibnd)+real(conjg(psi(ig,ibnd))*psic(dffts%nl(igk_k(ig,1))))
188          enddo
189          call mp_sum(e_h(ibnd),world_comm)
190          write(stdout,*) 'energies_h :', ibnd, e_h(ibnd)*rytoev
191
192          CALL stop_clock( 'secondfft' )
193
194       enddo!
195
196       !
197       !
198       !
199       RETURN
200       !
201     END SUBROUTINE energies_xc_k
202     SUBROUTINE energies_xc_gamma
203
204       USE uspp, ONLY : okvan
205       USE wannier_gw, ONLY : becp_gw, restart_gww,l_whole_s,l_verbose,&
206                               &l_scissor,scissor,num_nbndv,num_nbnds
207      ! USE realus,  ONLY : adduspos_gamma_r
208       USE wvfct,    ONLY : npwx,npw,nbnd, et,g2kin
209       USE wavefunctions, ONLY : evc
210       USE klist,                ONLY : xk
211       USE mp, ONLY : mp_sum
212       USE mp_world, ONLY : world_comm
213       USE gvect,  ONLY : gstart,g
214       USE constants, ONLY : rytoev
215       USE becmod,           ONLY : becp, calbec,allocate_bec_type,deallocate_bec_type
216       USE cell_base,            ONLY : tpiba2
217       USE io_global, ONLY : ionode
218       USE io_files, ONLY :prefix,tmp_dir
219     USE exx, ONLY : exxalfa
220
221       implicit none
222
223       INTEGER, EXTERNAL :: find_free_unit
224
225       INTEGER :: ibnd,jbnd,ir,ig
226       INTEGER :: iunwfcreal,iunu
227       REAL(kind=DP) :: etxc,vtxc,ehart, charge
228       REAL(kind=DP), ALLOCATABLE :: psi_r(:),psi_rs(:)
229       LOGICAL :: exst
230       REAL(kind=DP), ALLOCATABLE :: rho_fake_core(:)
231       COMPLEX(kind=DP), ALLOCATABLE :: hpsi(:,:)
232       REAL(kind=DP), ALLOCATABLE :: exact_x(:)
233       REAL(kind=DP), ALLOCATABLE :: e_hub(:)!Hubbard contribution to KS energies
234       REAL(kind=DP), ALLOCATABLE :: et_off(:,:)!off-diagonal energies
235
236
237       allocate(exact_x(nbnd))
238       allocate(e_hub(nbnd))
239       if(l_whole_s) then
240          allocate (et_off(nbnd,nbnd))
241       endif
242!if required calculates also the KS energies
243 !     if(restart_gww==-1) then
244         if(l_verbose) write(stdout,*) 'ATTENZIONE1'
245         FLUSH(stdout)
246         !allocate( becp%r( nkb, nbnd ) )
247         call allocate_bec_type ( nkb, nbnd, becp)
248         if(l_verbose) write(stdout,*) 'ATTENZIONE2'
249         FLUSH(stdout)
250
251         IF ( nkb > 0 )  CALL init_us_2( npw, igk_k(1,1), xk(1,1), vkb )
252         !call ccalbec( nkb, npwx, npw, nbnd, becp%r, vkb, evc )
253        !if(nkb> 0)CALL calbec ( npw, vkb, psi, becp, nbnd)
254
255        if(l_verbose)write(stdout,*) 'ATTENZIONE3'
256         FLUSH(stdout)
257
258         allocate(hpsi(npwx,nbnd))
259         if(l_verbose)write(stdout,*) 'ATTENZIONE4'
260         FLUSH(stdout)
261
262         g2kin(1:npw) = ( ( xk(1,1) + g(1,igk_k(1:npw,1)) )**2 + &
263              ( xk(2,1) + g(2,igk_k(1:npw,1)) )**2 + &
264              ( xk(3,1) + g(3,igk_k(1:npw,1)) )**2 ) * tpiba2
265
266
267         if(l_verbose)write(stdout,*) 'ATTENZIONE5'
268          FLUSH(stdout)
269
270
271 !         exxalfa=0.d0!ATTENZIONE
272         call h_psi( npwx, npw, nbnd, psi, hpsi )
273         et(:,ispin)=0.d0
274         if(l_verbose)write(stdout,*) 'ATTENZIONE6'
275         if(l_verbose)write(stdout,*) 'EXXALFA', exxalfa
276         FLUSH(stdout)
277
278          do ibnd=1,nbnd
279
280           !   call dgemm('T','N',1,1,2*npw,2.d0,evc(:,ibnd),2*npwx,hpsi(:,ibnd),2*npwx,&
281           !        &0.d0,et(ibnd,1),1)
282              do ig=1,npw
283                 et(ibnd,ispin)=et(ibnd,ispin)+2.d0*dble(conjg(evc(ig,ibnd))*hpsi(ig,ibnd))
284              enddo
285              if(gstart==2) then
286                et(ibnd,ispin)=et(ibnd,ispin)-dble(conjg(evc(1,ibnd))*hpsi(1,ibnd))
287             endif
288          enddo
289          call mp_sum(et(:,ispin),world_comm)
290          if(l_scissor) then
291             et(1:num_nbndv(ispin),ispin)=et(1:num_nbndv(ispin),ispin)+scissor(1)/rytoev
292             et(num_nbndv(ispin)+1:num_nbnds,ispin)=et(num_nbndv(ispin)+1:num_nbnds,ispin)+scissor(2)/rytoev
293          endif
294
295          if(l_verbose)write(stdout,*) 'ATTENZIONE7'
296          FLUSH(stdout)
297!if required calculate Hubbard U contribution to eigen-energies
298         e_hub(:)=0.d0
299         if ( lda_plus_u ) then
300            hpsi(:,:)=(0.d0,0.d0)
301            CALL vhpsi( npwx, npw, nbnd, psi, hpsi )
302
303            do ibnd=1,nbnd
304
305               do ig=1,npw
306                  e_hub(ibnd)=e_hub(ibnd)+2.d0*dble(conjg(psi(ig,ibnd))*hpsi(ig,ibnd))
307               enddo
308               if(gstart==2) then
309                  e_hub(ibnd)=e_hub(ibnd)-dble(conjg(psi(1,ibnd))*hpsi(1,ibnd))
310               endif
311            enddo
312            call mp_sum(e_hub(:),world_comm)
313            do ibnd=1,nbnd
314               write(stdout,*) 'Hubbard U energy:',ibnd,e_hub(ibnd)*rytoev
315            enddo
316            FLUSH(stdout)
317
318         endif
319          do ibnd=1,nbnd
320             write(stdout,*) 'KS energy:',ibnd,et(ibnd,ispin)*rytoev
321          enddo
322          FLUSH(stdout)
323
324!in case of hybrid functionals and HF we have to calculated also the exact exchange part
325
326          if(dft_is_hybrid()) then
327!NOT_TO_BE_INCLUDED_START
328             hpsi(:,:)=(0.d0,0.d0)
329             call vexx( npwx, npw, nbnd, psi, hpsi )
330             do ibnd=1,nbnd
331                call dgemm('T','N',1,1,2*npw,2.d0,evc(:,ibnd),2*npwx,hpsi(:,ibnd),2*npwx,&
332                     &0.d0,exact_x(ibnd),1)
333                if(gstart==2) then
334                   exact_x(ibnd)=exact_x(ibnd)-dble(conjg(evc(1,ibnd))*hpsi(1,ibnd))
335                endif
336                call mp_sum(exact_x(ibnd),world_comm)
337                write(stdout,*) 'Exact exchange :',ibnd, exact_x(ibnd)
338             enddo
339!NOT_TO_BE_INCLUDED_END
340          end if
341
342 !         deallocate(hpsi,becp%r)
343          call deallocate_bec_type ( becp)
344!       endif
345
346
347
348       allocate(psi_r(dfftp%nnr),psi_rs(dfftp%nnr))
349
350       iunwfcreal=find_free_unit()
351       CALL diropn( iunwfcreal, 'real_whole', dffts%nnr, exst )
352
353
354!calculate xc potential on fine grid
355
356
357       if(.not.allocated(vr)) write(stdout,*) 'vr not allocated'
358       allocate(rho_fake_core(dfftp%nnr))
359       rho_fake_core(:)=0.d0
360       !
361       if (dft_is_meta()) then
362      !    call v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v%of_r, v%kin_r )
363       else
364          CALL v_xc( rho, rho_core, rhog_core, etxc, vtxc, vr )
365       endif
366       !
367     deallocate(rho_fake_core)
368
369
370       if(l_whole_s) then
371!NOT_TO_BE_INCLUDED_START
372          allocate(hpsi(npwx,nbnd))
373          hpsi(:,:)=(0.d0,0.d0)
374          CALL vloc_psi_gamma ( npwx, npw, nbnd, evc, vr(1,ispin), hpsi )
375          call dgemm('T','N',nbnd,nbnd,2*npw,2.d0,evc,2*npwx,hpsi,2*npwx,&
376               &0.d0,et_off,nbnd)
377          if(gstart==2) then
378             do ibnd=1,nbnd
379                do jbnd=1,nbnd
380                   et_off(ibnd,jbnd)=et_off(ibnd,jbnd)-dble(conjg(evc(1,ibnd))*hpsi(1,jbnd))
381                enddo
382             enddo
383          endif
384          deallocate(hpsi)
385          call mp_sum(et_off,world_comm)
386!write on file
387          if(ionode) then
388             iunu = find_free_unit()
389             if(ispin==1) then
390                open(unit=iunu,file=trim(tmp_dir)//trim(prefix)//'.exc_off',status='unknown',form='unformatted')
391             else
392                open(unit=iunu,file=trim(tmp_dir)//trim(prefix)//'.exc_off2',status='unknown',form='unformatted')
393             endif
394             write(iunu) nbnd
395             do ibnd=1,nbnd
396                write(iunu) et_off(1:nbnd,ibnd)
397             enddo
398             close(iunu)
399          endif
400!NOT_TO_BE_INCLUDED_END
401       endif
402
403
404
405
406       do ibnd=1,m!loop on states
407!read from disk wfc on coarse grid
408         CALL davcio( psi_rs,dffts%nnr,iunwfcreal,ibnd+(ispin-1)*nbnd,-1)
409         if(doublegrid) then
410           call fft_interpolate(dffts, psi_rs, dfftp, psi_r) ! interpolate from smooth to dense
411         else
412           psi_r(:)=psi_rs(:)
413         endif
414
415         do ir=1,dfftp%nnr
416            psi_r(ir)=psi_r(ir)**2.d0
417         enddo
418
419         !if(okvan) call adduspos_gamma_r(ibnd,ibnd,psi_r,1,becp_gw(:,ibnd),becp_gw(:,ibnd))
420
421         e_xc(ibnd)=0.d0
422         do ir=1,dfftp%nnr
423           e_xc(ibnd)=e_xc(ibnd)+psi_r(ir)*vr(ir,ispin)!the 1 is for the spin NOT IMPLEMENTED YET
424         enddo
425         e_xc(ibnd)=e_xc(ibnd)/dble(dfftp%nr1*dfftp%nr2*dfftp%nr3)
426
427         call mp_sum(e_xc(ibnd),world_comm)
428
429!ifrequired add the contribution from exact exchange for hybrids and HF
430         if(dft_is_hybrid()) then
431!NOT_TO_BE_INCLUDED_START
432            e_xc(ibnd)=e_xc(ibnd)+exact_x(ibnd)
433!NOT_TO_BE_INCLUDED_END
434         endif
435
436         write(stdout,*) 'Routine energies_xc :', ibnd, e_xc(ibnd)*rytoev
437
438!now hartree term
439
440
441      enddo
442
443!if required add to e_xc Hubbard U terms
444
445      if(lda_plus_u) then
446
447         e_xc(1:nbnd)=e_xc(1:nbnd)+e_hub(1:nbnd)
448      endif
449
450
451       vr(:,:)=0.d0
452
453       call  v_h(rho%of_g , ehart, charge, vr )
454
455
456        do ibnd=1,m!loop on states
457!read from disk wfc on coarse grid
458           CALL davcio( psi_rs,dffts%nnr,iunwfcreal,ibnd+(ispin-1)*nbnd,-1)
459           if(doublegrid) then
460              call fft_interpolate(dffts, psi_rs, dfftp, psi_r) ! interpolate from smooth to dense
461           else
462              psi_r(:)=psi_rs(:)
463           endif
464
465           do ir=1,dfftp%nnr
466              psi_r(ir)=psi_r(ir)**2.d0
467           enddo
468
469           !if(okvan) call adduspos_gamma_r(ibnd,ibnd,psi_r,1,becp_gw(:,ibnd),becp_gw(:,ibnd))
470
471           e_h(ibnd)=0.d0
472           do ir=1,dfftp%nnr
473              e_h(ibnd)=e_h(ibnd)+psi_r(ir)*vr(ir,ispin)
474           enddo
475           e_h(ibnd)=e_h(ibnd)/dble(dfftp%nr1*dfftp%nr2*dfftp%nr3)
476
477           call mp_sum(e_h(ibnd),world_comm)
478           write(stdout,*) 'Routine energies_h :', ibnd, e_h(ibnd)*rytoev
479
480!now hartree term
481
482
483        enddo
484
485
486
487
488
489       deallocate(psi_r,psi_rs)
490       deallocate(exact_x)
491      close(iunwfcreal)
492      deallocate(e_hub)
493      if(l_whole_s) then
494!NOT_TO_BE_INCLUDED_START
495         deallocate(et_off)
496!NOT_TO_BE_INCLUDED_END
497      endif
498
499      return
500
501     END SUBROUTINE energies_xc_gamma
502     !
503END SUBROUTINE energies_xc
504
505
506SUBROUTINE write_energies_xc(e_xc)
507
508  USE kinds, ONLY : DP
509  USE wannier_gw, ONLY : num_nbnds, l_verbose
510  USE io_files, ONLY : prefix,tmp_dir
511  USE io_global, ONLY : ionode
512  USE wvfct,    ONLY : nbnd
513  USE lsda_mod, ONLY : nspin
514
515  implicit none
516
517  INTEGER, EXTERNAL :: find_free_unit
518
519  REAL(kind=DP) :: e_xc(nbnd,nspin)!exchange and correlation energies
520
521  INTEGER :: iunu, iw
522
523
524  if(ionode) then
525     iunu = find_free_unit()
526
527     open(unit=iunu,file=trim(tmp_dir)//trim(prefix)//'.dft_xc',status='unknown',form='unformatted')
528
529     write(iunu) nbnd
530
531
532     do iw=1,nbnd
533        write(iunu) e_xc(iw,1)
534        if(l_verbose) WRITE(*,*) 'SCRITTO e_XC 1', e_xc(iw,1)
535     enddo
536     if(nspin==2) then
537        do iw=1,nbnd
538           write(iunu) e_xc(iw,2)
539           if(l_verbose) WRITE(*,*) 'SCRITTO e_XC 2', e_xc(iw,2)
540        enddo
541     endif
542     close(iunu)
543  endif
544
545  return
546
547END SUBROUTINE write_energies_xc
548