calc_exx2.f90 8.52 KB
Newer Older
Marco Govoni's avatar
Marco Govoni committed
1
!
2
! Copyright (C) 2015-2017 M. Govoni 
Marco Govoni's avatar
Marco Govoni committed
3
4
5
6
7
8
9
10
11
12
13
! 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 .
!
! This file is part of WEST.
!
! Contributors to this file: 
! Marco Govoni
!
!-----------------------------------------------------------------------
14
SUBROUTINE calc_exx2( sigma_exx, nb1, nb2 )
15
16
17
18
19
20
21
22
23
24
  !-----------------------------------------------------------------------
  !
  ! store in sigma_exx(n,iks) = < n,iks | V_exx | n,iks >     n = nb1, nb2
  !
  USE kinds,                ONLY : DP 
  USE mp_global,            ONLY : inter_image_comm,intra_bgrp_comm,my_image_id
  USE mp,                   ONLY : mp_sum,mp_bcast
  USE io_global,            ONLY : stdout, ionode
  USE io_files,             ONLY : nwordwfc, iunwfc
  USE scf,                  ONLY : rho, rho_core, rhog_core
25
  USE gvect,                ONLY : g,nl,gstart,ngm_g,ngm
26
27
28
29
30
31
32
33
34
  USE gvecw,                ONLY : gcutw
  USE cell_base,            ONLY : tpiba2,omega,tpiba,at,alat
  USE fft_base,             ONLY : dfftp,dffts
  USE fft_interfaces,       ONLY : fwfft, invfft
  USE constants,            ONLY : tpi,fpi,rytoev,e2
  USE pwcom,                ONLY : npw,npwx,et,nks,current_spin,isk,xk,nbnd,lsda,igk_k,g2kin,nspin,current_k,ngk
  USE fft_at_gamma,         ONLY : single_invfft_gamma,single_fwfft_gamma
  USE fft_at_k,             ONLY : single_invfft_k,single_fwfft_k
  USE wavefunctions_module, ONLY : evc,psic,psic_nc
35
  USE westcom,              ONLY : iuwfc,lrwfc,npwq,nbnd_occ
36
37
38
39
40
41
42
43
44
45
46
47
  USE control_flags,        ONLY : gamma_only
  USE noncollin_module,     ONLY : noncolin,npol 
  USE buffers,              ONLY : get_buffer
  USE funct,                ONLY : dft_is_hybrid,init_dft_exxrpa,stop_exx
  USE exx,                  ONLY : x_gamma_extrapolation,exxdiv_treatment,exx_grid_init,exx_div_check,&
                                   &deallocate_exx,exxinit,vexx,exx_grid_initialized
  USE uspp,                 ONLY : vkb,nkb
  USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
  USE io_push,              ONLY : io_push_bar
  USE class_idistribute,    ONLY : idistribute
  USE coulomb_vcut_module,  ONLY : vcut_init, vcut_type, vcut_info, &
                                   vcut_get,  vcut_spheric_get, vcut_destroy
48
  USE types_bz_grid,        ONLY : k_grid, q_grid, compute_phase
49
  USE types_coulomb,        ONLY : pot3D
50
51
52
53
54
55
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  INTEGER, INTENT(IN) :: nb1, nb2
56
  REAL(DP),INTENT(OUT) :: sigma_exx( nb1:nb2, k_grid%nps) 
57
58
59
60
61
62
63
  !
  ! Workspace
  !
  COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:),pertr_nc(:,:)
  COMPLEX(DP), ALLOCATABLE :: evckmq(:,:), phase(:)
  REAL(DP), EXTERNAL :: DDOT
  COMPLEX(DP), EXTERNAL :: ZDOTC
64
  INTEGER :: ib,iv,i1,ir,iks,ik,is,ig,iv_glob,iq,ikqs
65
66
67
68
69
70
  INTEGER :: nbndval
  INTEGER :: npwkq
  TYPE(idistribute) :: vband
  TYPE(bar_type) :: barra
  INTEGER :: barra_load
  LOGICAL :: l_gammaq
71
  REAL(DP) :: g0(3), peso
72
73
74
75
76
77
  !
  WRITE(stdout,'(5x,a)') ' '
  CALL io_push_bar()
  WRITE(stdout,'(5x,a)') '(X)-Sigma'
  CALL io_push_bar()
  !
78
  ALLOCATE( pertg( ngm ) )
79
80
81
82
83
84
85
  IF (gamma_only) THEN
     peso = 2._DP
  ELSE
     peso = 1._DP
     ALLOCATE( phase(dffts%nnr) )
     ALLOCATE( evckmq(npwx*npol,nbnd) )
  ENDIF
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
  IF(noncolin) THEN 
     ALLOCATE( pertr_nc( dffts%nnr, npol ) )
  ELSE
     ALLOCATE( pertr( dffts%nnr ) )
  ENDIF
  !
  !
  ! Set to zero
  !
  sigma_exx = 0._DP
  !
  barra_load = k_grid%nps * ( nb2-nb1 + 1 )
  CALL start_bar_type( barra, 'sigmax', barra_load )
  !
  ! LOOP 
  !
  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
103
104
     !
     ik = k_grid%ip(iks) 
105
     is = k_grid%is(iks)
106
107
108
109
110
111
112
113
114
     !
     ! ... Set k-point, spin, kinetic energy, needed by Hpsi
     !
     current_k = iks
     IF ( lsda ) current_spin = isk(iks)
     call g2_kin( iks )
     !
     ! ... More stuff needed by the hamiltonian: nonlocal projectors
     !
115
     IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), k_grid%p_cart(1,ik), vkb )
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
     npw = ngk(iks)
     !
     ! ... read in wavefunctions from the previous iteration
     !
     IF(k_grid%nps>1) THEN
        !iuwfc = 20
        !lrwfc = nbnd * npwx * npol 
        !!CALL get_buffer( evc, nwordwfc, iunwfc, iks )
        IF(my_image_id==0) CALL get_buffer( evc, lrwfc, iuwfc, iks )
        !CALL mp_bcast(evc,0,inter_image_comm)
        !CALL davcio(evc,lrwfc,iuwfc,iks,-1)
        CALL mp_bcast(evc,0,inter_image_comm)
     ENDIF
!     !
!     ! ... Needed for LDA+U
!     !
!     IF ( nks > 1 .AND. lda_plus_u .AND. (U_projection .NE. 'pseudo') ) &
!          CALL get_buffer ( wfcU, nwordwfcU, iunhub, iks )
!     !
!     current_k = iks
!     current_spin = isk(iks)
!     !
!     CALL gk_sort(xk(1,iks),ngm,g,gcutw,npw,igk,g2kin)
!     g2kin=g2kin*tpiba2
!     !
!     ! reads unperturbed wavefuctions psi_k in G_space, for all bands
!     !
!     !
!     CALL init_us_2 (npw, igk, xk (1, iks), vkb)
     !
     !
     DO ib = nb1, nb2
        !
        sigma_exx(ib,iks) = 0._DP
        !
151
152
153
        IF (gamma_only) THEN
           CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ib),psic,'Wave')
        ELSEIF (noncolin) THEN
154
           CALL single_invfft_k(dffts,npw,npwx,evc(1     ,ib),psic_nc(1,1),'Wave',igk_k(1,current_k))
155
           CALL single_invfft_k(dffts,npw,npwx,evc(npwx+1,ib),psic_nc(1,2),'Wave',igk_k(1,current_k))
156
157
158
159
        ELSE
           CALL single_invfft_k(dffts,npw,npwx,evc(1,ib),psic,'Wave',igk_k(1,current_k))
        ENDIF
        !
160
        DO iq = 1, q_grid%np
161
           !
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
           IF (gamma_only) THEN
              l_gammaq = .TRUE.
              CALL pot3D%init('Dense',.FALSE.,'gb')
              nbndval = nbnd_occ(iks)
           ELSE
              l_gammaq = q_grid%l_pIsGamma(iq)
              CALL pot3D%init('Dense',.FALSE.,'gb',iq)
              !
              CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
              CALL compute_phase( g0, 'cart', phase )
              !
              nbndval = nbnd_occ(ikqs)
              npwkq = ngk(ikqs)
              IF ( my_image_id == 0 ) CALL get_buffer( evckmq, lrwfc, iuwfc, ikqs )
              CALL mp_bcast( evckmq, 0, inter_image_comm )
           ENDIF
178
           !
179
180
           vband = idistribute()
           CALL vband%init(nbndval,'i','nbndval',.FALSE.)
181
182
183
184
185
186
           !
           DO iv = 1, vband%nloc
              !
              iv_glob = vband%l2g(iv)
              !
              ! Bring it to R-space
187
188
189
190
191
192
193
              IF (gamma_only) THEN
                 CALL single_invfft_gamma(dffts,npw,npwx,evc(1,iv_glob),pertr,'Wave')
                 DO ir=1,dffts%nnr
                    pertr(ir)=psic(ir)*pertr(ir)
                 ENDDO
                 CALL single_fwfft_gamma(dffts,ngm,ngm,pertr,pertg,'Dense')
              ELSEIF(noncolin) THEN
194
195
196
                 CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1     ,iv_glob),pertr_nc(1,1),'Wave',igk_k(1,ikqs))
                 CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1+npwx,iv_glob),pertr_nc(1,2),'Wave',igk_k(1,ikqs))
                 DO ir=1,dffts%nnr 
197
                    pertr_nc(ir,1)=DCONJG(pertr_nc(ir,1)*phase(ir))*psic_nc(ir,1)+DCONJG(pertr_nc(ir,2)*phase(ir))*psic_nc(ir,2)
198
                 ENDDO
199
                 CALL single_fwfft_k(dffts,ngm,ngm,pertr_nc(1,1),pertg,'Dense') ! no igk
200
201
202
203
204
              ELSE
                 CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1,iv_glob),pertr,'Wave',igk_k(1,ikqs))
                 DO ir=1,dffts%nnr 
                    pertr(ir)=DCONJG(pertr(ir)*phase(ir)) * psic(ir)
                 ENDDO
205
                 CALL single_fwfft_k(dffts,ngm,ngm,pertr,pertg,'Dense') ! no igk
206
207
              ENDIF 
              !
208
              DO ig = 1,ngm
209
                 pertg(ig) = pertg(ig) * pot3D%sqvc(ig) 
210
              ENDDO
211
              sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - peso*DDOT( 2*ngm, pertg(1), 1, pertg(1), 1)/omega*q_grid%weight(iq)
212
              !IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
213
              IF( ib == iv_glob .AND. gstart == 2 .AND. l_gammaq ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
              !
           ENDDO ! iv
           !
        ENDDO ! iq
        !
        CALL update_bar_type( barra, 'sigmax', 1 )
        !
     ENDDO ! ib
     !
  ENDDO ! iks
  !
  CALL stop_bar_type( barra, 'sigmax' )
  !
  CALL mp_sum( sigma_exx, intra_bgrp_comm )
  CALL mp_sum( sigma_exx, inter_image_comm ) 
  !
  DEALLOCATE( pertg ) 
  IF( noncolin ) THEN 
    DEALLOCATE( pertr_nc ) 
  ELSE
    DEALLOCATE( pertr ) 
  ENDIF
236
237
238
239
  IF (.NOT.gamma_only) THEN
     DEALLOCATE( phase )
     DEALLOCATE( evckmq )
  ENDIF
240
241
  !
END SUBROUTINE