do_sxx.f90 12.4 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
14
15
16
17
18
19
! 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, Nicholas Brawand 
!
!----------------------------------------------------------------------------
SUBROUTINE do_sxx ( )
  !----------------------------------------------------------------------------
  !
  USE kinds,                 ONLY : DP
  USE uspp,                  ONLY : vkb,nkb
  USE io_global,             ONLY : stdout
Govoni's avatar
Govoni committed
20
  USE pwcom,                 ONLY : current_spin,wk,nks,nelup,neldw,isk,igk_k,xk,npw,npwx,lsda,nkstot,current_k,ngk,et
Marco Govoni's avatar
Marco Govoni committed
21
  USE io_push,               ONLY : io_push_title,io_push_bar
22
  USE westcom,               ONLY : iuwfc,lrwfc,westpp_range,westpp_save_dir,nbnd_occ,iks_l2g,westpp_epsinfty,dvg,ev,&
23
                                  & npwq,npwqx,fftdriver,logfile,ngq,westpp_n_pdep_eigen_to_use
Marco Govoni's avatar
Marco Govoni committed
24
25
26
27
28
29
30
31
32
33
34
  USE mp_global,             ONLY : inter_image_comm,my_image_id,intra_bgrp_comm
  USE mp,                    ONLY : mp_bcast,mp_sum
  USE fft_base,              ONLY : dfftp,dffts
  USE wvfct,                 ONLY : nbnd
  USE buffers,               ONLY : get_buffer
  USE wavefunctions_module,  ONLY : evc,psic,psic_nc
  USE bar,                   ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
  USE fft_at_gamma,          ONLY : single_invfft_gamma,single_fwfft_gamma
  USE fft_at_k,              ONLY : single_invfft_k,single_fwfft_k
  USE distribution_center,   ONLY : pert
  USE control_flags,         ONLY : gamma_only 
35
  USE gvect,                 ONLY : g,nl,gstart,ngm_g,ngm
Govoni's avatar
Govoni committed
36
  USE cell_base,             ONLY : omega,at,alat
Marco Govoni's avatar
Marco Govoni committed
37
38
39
  USE noncollin_module,      ONLY : noncolin,npol 
  USE mp_world,              ONLY : mpime,root
  USE constants,             ONLY : rytoev
Govoni's avatar
Govoni committed
40
  USE json_module,           ONLY : json_file
41
42
  USE pdep_db,               ONLY : pdep_db_read
  USE types_bz_grid,         ONLY : k_grid, q_grid, compute_phase
43
  USE types_coulomb,         ONLY : pot3D
Marco Govoni's avatar
Marco Govoni committed
44
45
46
47
48
  !
  IMPLICIT NONE
  !
  ! ... LOCAL variables
  !
Marco Govoni's avatar
Marco Govoni committed
49
  INTEGER :: ir, ip, ig, iks, ib, iv, ip_glob, ik, is, ikqs, ikq, iq, nbndval, npwkq
Marco Govoni's avatar
Marco Govoni committed
50
  COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:),pertr_nc(:,:)
51
52
53
  COMPLEX(DP), ALLOCATABLE :: evckmq(:,:), phase(:)
  LOGICAL :: l_gammaq
  REAL(DP) :: g0(3)
Marco Govoni's avatar
Marco Govoni committed
54
55
56
57
58
  TYPE(bar_type) :: barra
  REAL(DP),ALLOCATABLE :: sigma_exx( :, : ) 
  REAL(DP),ALLOCATABLE :: sigma_sxx( :, : ) 
  REAL(DP) :: peso
  REAL(DP), EXTERNAL :: DDOT
Govoni's avatar
Govoni committed
59
  CHARACTER(LEN=6) :: myglobk
Marco Govoni's avatar
Marco Govoni committed
60
61
62
  REAL(DP),ALLOCATABLE :: out_tab(:,:)
  COMPLEX(DP),ALLOCATABLE :: zproj(:,:)
  REAL(DP),ALLOCATABLE :: dproj(:,:)
Govoni's avatar
Govoni committed
63
64
  TYPE(json_file) :: json
  INTEGER :: iunit
65
  LOGICAL :: l_print_pdep_read
Marco Govoni's avatar
Marco Govoni committed
66
67
68
  !
  CALL io_push_title("(S)creened eXact eXchange")
  !
69
70
  ALLOCATE( sigma_exx( westpp_range(1):westpp_range(2), k_grid%nps) )
  ALLOCATE( sigma_sxx( westpp_range(1):westpp_range(2), k_grid%nps) )
Marco Govoni's avatar
Marco Govoni committed
71
72
73
74
  !
  sigma_exx = 0._DP
  sigma_sxx = 0._DP
  !
75
76
  !ALLOCATE( pertg(npwqx) )
  ALLOCATE( pertg(ngm) )
Marco Govoni's avatar
Marco Govoni committed
77
78
79
80
81
82
83
84
85
86
87
88
  IF(noncolin) THEN 
     ALLOCATE( pertr_nc( dffts%nnr, npol ) )
  ELSE
     ALLOCATE( pertr( dffts%nnr ) )
  ENDIF
  !
  IF( gamma_only ) THEN 
     peso = 2._DP  
     ALLOCATE( dproj( 1, pert%nloc ) )
  ELSE
     peso = 1._DP
     ALLOCATE( zproj( 1, pert%nloc ) )
89
90
     ALLOCATE( phase(dffts%nnr) )
     ALLOCATE( evckmq(npwx*npol,nbnd) )
Marco Govoni's avatar
Marco Govoni committed
91
92
  ENDIF
  !
93
  CALL start_bar_type( barra, 'westpp', k_grid%nps * (westpp_range(2)-westpp_range(1)+1)  ) 
Marco Govoni's avatar
Marco Govoni committed
94
  !
95
96
97
98
  DO iks = 1, k_grid%nps  ! KPOINT-SPIN LOOP
     !
     ik = k_grid%ip(iks)
     is = k_grid%is(iks)
Marco Govoni's avatar
Marco Govoni committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
     !
     ! ... 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
     !
     !IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), xk(1,iks), vkb )
     npw = ngk(iks)
     !
     ! ... read in wavefunctions from the previous iteration
     !
113
     IF(k_grid%nps>1) THEN
Marco Govoni's avatar
Marco Govoni committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
        !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
     !
     !nbndval = nbnd_occ(iks)
     !
     DO ib = 1, nbnd
        !
        IF( ib < westpp_range(1) .OR. ib > westpp_range(2) ) CYCLE 
        !
        IF(gamma_only) THEN
           CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ib),psic,'Wave') 
        ELSEIF(noncolin) THEN
           CALL single_invfft_k(dffts,npw,npwx,evc(1     ,ib),psic_nc(1,1),'Wave',igk_k(1,current_k))
           CALL single_invfft_k(dffts,npw,npwx,evc(1+npwx,ib),psic_nc(1,2),'Wave',igk_k(1,current_k))
        ELSE
           CALL single_invfft_k(dffts,npw,npwx,evc(1,ib),psic,'Wave',igk_k(1,current_k))
        ENDIF
        !
138
        DO iq = 1, q_grid%np
Marco Govoni's avatar
Marco Govoni committed
139
           !
140
141
142
143
144
145
           IF ( gamma_only ) THEN
              !
              l_gammaq = .TRUE.
              nbndval = nbnd_occ(iks)
              CALL pot3D%init('Dense',.FALSE.,'gb')
              !
Marco Govoni's avatar
Marco Govoni committed
146
           ELSE
147
148
149
              !
              l_gammaq = q_grid%l_pIsGamma(iq)
              !
150
151
152
              !CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )  !M
              CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), 'cart', ikq, g0 )        
              ikqs = k_grid%ipis2ips(ikq,is)                                                        
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
              CALL compute_phase( g0, 'cart', phase )
              !
              nbndval = nbnd_occ(ikqs)
              npwq = ngq(iq)
              npwkq = ngk(ikqs)
              !
              CALL pot3D%init('Dense',.FALSE.,'gb',iq)
              !
              IF ( my_image_id == 0 ) CALL get_buffer( evckmq, lrwfc, iuwfc, ikqs )
              CALL mp_bcast( evckmq, 0, inter_image_comm )
              !
              IF (iks==1 .AND. ib==1 .AND. iq==1) THEN 
                 l_print_pdep_read = .TRUE.
              ELSE
                 l_print_pdep_read = .FALSE.
              ENDIF
              !
           ENDIF
Marco Govoni's avatar
Marco Govoni committed
171
           !
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
           DO iv = 1, nbndval
              !
              ! Bring it to R-space
              IF(gamma_only) THEN
                 CALL single_invfft_gamma(dffts,npw,npwx,evc(1,iv),pertr,'Wave')
                 DO ir=1,dffts%nnr 
                    pertr(ir)=psic(ir)*pertr(ir)
                 ENDDO
                 CALL single_fwfft_gamma(dffts,npwq,npwqx,pertr,pertg,TRIM(fftdriver))
              ELSEIF(noncolin) THEN
                 CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1     ,iv),pertr_nc(1,1),'Wave',igk_k(1,ikqs))
                 CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1+npwx,iv),pertr_nc(1,2),'Wave',igk_k(1,ikqs))
                 DO ir=1,dffts%nnr 
                    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)
                 ENDDO
                 !CALL single_fwfft_k(dffts,npwq,npwqx,pertr_nc(1,1),pertg,'Dense') ! no igk
                 CALL single_fwfft_k(dffts,ngm,ngm,pertr_nc(1,1),pertg,'Dense') ! no igk
              ELSE
                 CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1,iv),pertr,'Wave',igk_k(1,ikqs))
                 DO ir=1,dffts%nnr 
                    pertr(ir)=DCONJG(pertr(ir)*phase(ir)) * psic(ir)
                 ENDDO
                 !CALL single_fwfft_k(dffts,npwq,npwqx,pertr,pertg,'Dense') ! no igk
                 CALL single_fwfft_k(dffts,ngm,ngm,pertr,pertg,'Dense') ! no igk
              ENDIF 
              !
              !DO ig = 1,npwq
              DO ig = 1,ngm
                 pertg(ig) = pertg(ig) * pot3D%sqvc(ig) 
Marco Govoni's avatar
Marco Govoni committed
201
              ENDDO
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
              !sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - peso * DDOT( 2*npwq, pertg(1), 1, pertg(1), 1) / omega
              sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - peso*DDOT(2*ngm, pertg(1), 1, pertg(1), 1)/omega * q_grid%weight(iq)
              !IF(gstart==2) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) + REAL( pertg(1), KIND = DP )**2 / omega
              IF( ib == iv .AND. gstart == 2 .AND. l_gammaq ) sigma_exx( ib, iks ) = sigma_exx( ib, iks ) - pot3D%div
              !
              ! -- < SXX >
              !
              IF (gamma_only) THEN
                 CALL pdep_db_read(westpp_n_pdep_eigen_to_use)
              ELSE
                 CALL pdep_db_read(westpp_n_pdep_eigen_to_use,iq,l_print_pdep_read)
              ENDIF
              !
              IF( gamma_only ) THEN  
                 CALL glbrak_gamma( pertg, dvg, dproj, npwq, npwqx, 1, pert%nloc, 1, npol)
                 CALL mp_sum( dproj, intra_bgrp_comm )
                 DO ip = 1, pert%nloc
                    ip_glob = pert%l2g(ip)
                    sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - dproj(1,ip)**2 * (ev(ip_glob)/(1._DP-ev(ip_glob))) / omega
                 ENDDO
                 IF( ib == iv ) sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - (1._DP/westpp_epsinfty-1._DP) * pot3D%div
              ELSE
                 CALL glbrak_k( pertg, dvg, zproj, npwq, npwqx, 1, pert%nloc, 1, npol)
                 CALL mp_sum( zproj, intra_bgrp_comm )
                 DO ip = 1, pert%nloc
                    ip_glob = pert%l2g(ip)
                    sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - REAL(zproj(1,ip)*CONJG(zproj(1,ip)),KIND=DP) &
                                         & * (ev(ip_glob)/(1._DP-ev(ip_glob))) / omega * q_grid%weight(iq)
                 ENDDO
                 !
                 IF( ib == iv ) sigma_sxx( ib, iks ) = sigma_sxx( ib, iks ) - (1._DP/westpp_epsinfty-1._DP) * pot3D%div
              ENDIF 
              ! -- </ SXX > 
              !
           ENDDO ! iv
Marco Govoni's avatar
Marco Govoni committed
237
           !
238
        ENDDO ! iq
Marco Govoni's avatar
Marco Govoni committed
239
240
241
        !
        CALL update_bar_type( barra,'westpp', 1 )
        !
242
     ENDDO ! ib
Marco Govoni's avatar
Marco Govoni committed
243
     !
244
  ENDDO ! iks
Marco Govoni's avatar
Marco Govoni committed
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
  !
  CALL stop_bar_type( barra, 'westpp' )
  !
  CALL mp_sum( sigma_exx, intra_bgrp_comm )
  CALL mp_sum( sigma_sxx, inter_image_comm )
  !
  sigma_sxx = sigma_exx + sigma_sxx 
  !
  DEALLOCATE( pertg ) 
  IF( noncolin ) THEN 
    DEALLOCATE( pertr_nc ) 
  ELSE
    DEALLOCATE( pertr ) 
  ENDIF
  IF( gamma_only ) THEN
     DEALLOCATE( dproj )
  ELSE
     DEALLOCATE( zproj )
263
264
     DEALLOCATE( evckmq )
     DEALLOCATE( phase )
Marco Govoni's avatar
Marco Govoni committed
265
266
267
268
  ENDIF  
  !
  ! Output it per k-point
  !
Govoni's avatar
Govoni committed
269
270
271
272
273
274
275
276
277
278
279
280
  IF(mpime==root) THEN
     CALL json%initialize()
     CALL json%load_file(filename=TRIM(logfile))
  ENDIF
  !
  ! STDOUT
  !
  WRITE(stdout,"(5X)")  
  CALL io_push_bar()
  WRITE(stdout,"(5X,a,1X,a,1X,a,1X,a,1X,a,1X,a)") &
  & 'K     ', 'B     ', '      Eks [eV]', '       Sx [eV]', '      Sxx [eV]', '        Sxx/Sx'
  CALL io_push_bar()
Govoni's avatar
Govoni committed
281
  !
Govoni's avatar
Govoni committed
282
  ALLOCATE(out_tab(westpp_range(2)-westpp_range(1)+1,5))
Govoni's avatar
Govoni committed
283
  !
284
  DO iks=1,k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
285
286
287
288
289
     DO ib = westpp_range(1), westpp_range(2)
        out_tab( ib - westpp_range(1) + 1, 1) = REAL( ib, KIND=DP) 
        out_tab( ib - westpp_range(1) + 1, 2) = et(ib,iks) * rytoev
        out_tab( ib - westpp_range(1) + 1, 3) = sigma_exx(ib,iks) * rytoev
        out_tab( ib - westpp_range(1) + 1, 4) = sigma_sxx(ib,iks) * rytoev
Govoni's avatar
Govoni committed
290
        out_tab( ib - westpp_range(1) + 1, 5) = sigma_sxx(ib,iks) / sigma_exx(ib,iks)
291
        WRITE(stdout,"(5X,i5.5,1X,i6.6,1X,1f14.6,1X,1f14.6,1X,1f14.6,1X,1f14.6)") &
Govoni's avatar
Govoni committed
292
        & iks, ib, et(ib,iks)*rytoev, sigma_exx(ib,iks)*rytoev, sigma_sxx(ib,iks)*rytoev, sigma_sxx(ib,iks)/sigma_exx(ib,iks)
Marco Govoni's avatar
Marco Govoni committed
293
     ENDDO
294
295
     IF (k_grid%nps>1.AND.iks<k_grid%nps) CALL io_push_bar()
     WRITE(myglobk,'(i5.5)') iks_l2g(iks)
Govoni's avatar
Govoni committed
296
297
298
299
     !
     !CALL serial_table_output(mpime==root,4000,'sxx_K'//myglobk,out_tab,&
     !& westpp_range(2)-westpp_range(1)+1,5,&
     !& (/'      band','    E0[eV]','    Sx[eV]','   Sxx[eV]','    Sxx/Sx'/))
Marco Govoni's avatar
Marco Govoni committed
300
     !
Govoni's avatar
Govoni committed
301
302
     IF(mpime==root) THEN 
        !
Govoni's avatar
Govoni committed
303
304
305
306
307
        CALL json%add('output.S.K'//TRIM(myglobk)//'.bandmap',out_tab(:,1))
        CALL json%add('output.S.K'//TRIM(myglobk)//'.eks',out_tab(:,2))
        CALL json%add('output.S.K'//TRIM(myglobk)//'.sx',out_tab(:,3))
        CALL json%add('output.S.K'//TRIM(myglobk)//'.sxx',out_tab(:,4))
        CALL json%add('output.S.K'//TRIM(myglobk)//'.fraction',out_tab(:,5))
Govoni's avatar
Govoni committed
308
309
        !
     ENDIF
Marco Govoni's avatar
Marco Govoni committed
310
311
  ENDDO
  DEALLOCATE(out_tab)
Govoni's avatar
Govoni committed
312
  CALL io_push_bar()
Marco Govoni's avatar
Marco Govoni committed
313
  !
Govoni's avatar
Govoni committed
314
315
  IF( mpime == root ) THEN 
     !
Govoni's avatar
Govoni committed
316
     OPEN( NEWUNIT=iunit, FILE=TRIM(logfile) )
Govoni's avatar
Govoni committed
317
318
319
320
321
322
323
     CALL json%print_file( iunit )
     CLOSE( iunit )
     !
     CALL json%destroy()
     !
  ENDIF
  !
Marco Govoni's avatar
Marco Govoni committed
324
325
326
327
  DEALLOCATE( sigma_exx )
  DEALLOCATE( sigma_sxx ) 
  !
END SUBROUTINE