solve_qp.f90 42.5 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
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
! 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
!
!-----------------------------------------------------------------------
SUBROUTINE solve_qp(l_secant,l_generate_plot)
  !-----------------------------------------------------------------------
  !
  USE control_flags,        ONLY : gamma_only 
  !
  IMPLICIT NONE
  !
  ! I/O 
  !
  LOGICAL,INTENT(IN) :: l_secant,l_generate_plot
  !
  IF( gamma_only ) THEN 
    CALL solve_qp_gamma( l_secant, l_generate_plot )
  ELSE
    CALL solve_qp_k( l_secant, l_generate_plot )
  ENDIF
  !
END SUBROUTINE 
!
!-----------------------------------------------------------------------
SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
  !-----------------------------------------------------------------------
  !
  ! ... This subroutine solves the DBS problem for GAMMA, at non-zero freqeuncies. 
  ! ... Perturbations are distributed according to the POT mpi_communicator
  !
  USE kinds,                ONLY : DP 
  USE westcom,              ONLY : n_pdep_eigen_to_use,n_lanczos,qp_bandrange,iks_l2g,imfreq_list_integrate,&
                                 & n_secant_maxiter,nbnd_occ,trev_secant,l_enable_lanczos,imfreq_list,n_imfreq,&
                                 & d_epsm1_ifr,z_epsm1_rfr,l_macropol,n_spectralf,ecut_spectralf, &
                                 & d_head_ifr,d_body1_ifr,d_body2_ifr,d_diago,z_head_rfr,z_body_rfr, &
                                 & sigma_z,sigma_eqplin,sigma_eqpsec,sigma_sc_eks,sigma_sc_eqplin,sigma_sc_eqpsec,sigma_diff, &
                                 & sigma_spectralf,sigma_freq
  USE mp_global,            ONLY : inter_image_comm,intra_bgrp_comm
  USE mp_world,             ONLY : mpime,root
  USE mp,                   ONLY : mp_barrier,mp_sum
  USE io_global,            ONLY : stdout, ionode
  USE pwcom,                ONLY : et,nks,current_spin,isk,xk,nbnd,lsda,g2kin,current_k
  USE io_push,              ONLY : io_push_title,io_push_bar
  USE constants,            ONLY : rytoev,pi
  USE west_io,              ONLY : serial_table_output
  USE distribution_center,  ONLY : pert,ifr,rfr,aband
  USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
  USE wfreq_io,             ONLY : readin_overlap,readin_solvegfreq,readin_solvehf
Govoni's avatar
Govoni committed
58
  USE wfreq_db,             ONLY : wfreq_db_write
59
  USE types_bz_grid,        ONLY : k_grid
Marco Govoni's avatar
Marco Govoni committed
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  LOGICAL,INTENT(IN) :: l_secant,l_generate_plot 
  !
  ! Workspace
  !
  REAL(DP),ALLOCATABLE :: sigma_hf(:,:)
  COMPLEX(DP),ALLOCATABLE :: sigma_cor_in (:,:)
  COMPLEX(DP),ALLOCATABLE :: sigma_cor_out(:,:)
  REAL(DP),ALLOCATABLE :: z_in(:,:)
  REAL(DP),ALLOCATABLE :: qp_energy(:,:)
  COMPLEX(DP),ALLOCATABLE :: sc(:,:,:) 
  REAL(DP),ALLOCATABLE :: en(:,:,:) 
  LOGICAL,ALLOCATABLE :: l_conv(:,:)
  REAL(DP),PARAMETER   :: eshift = 0.007349862_DP ! = 0.1 eV
  INTEGER :: k, ib, iks, ifixed, ip, glob_ip, ifreq, il, im, glob_im, glob_jp, glob_ifreq
  REAL(DP),ALLOCATABLE :: out_tab(:,:)
  CHARACTER(LEN=5) :: myglobk 
  CHARACTER(LEN=6) :: cifixed
  INTEGER :: notconv
  INTEGER,EXTERNAL :: get_nbndval
  INTEGER :: nbndval 
  REAL(DP),ALLOCATABLE :: dtemp(:)
  COMPLEX(DP),ALLOCATABLE :: ztemp(:)
  REAL(DP),ALLOCATABLE :: diago(:,:)
  REAL(DP),ALLOCATABLE :: braket(:,:,:)
  REAL(DP),ALLOCATABLE :: overlap(:,:)
  TYPE(bar_type) :: barra
  INTEGER :: barra_load
  CHARACTER(LEN=126) :: fname
  CHARACTER(LEN=5) :: ib_label, iks_label
  INTEGER,ALLOCATABLE :: un(:,:)
  REAL(DP) :: summed_sf
  ! 
  CALL start_clock('solve_qp')
  !
  ALLOCATE( imfreq_list_integrate( 2, ifr%nloc ) )
  ALLOCATE( dtemp( n_imfreq ) ) 
  !
  !
  dtemp = 0._DP
  DO ifreq = 1, ifr%nloc 
     glob_ifreq = ifr%l2g(ifreq)
     dtemp( glob_ifreq ) = imfreq_list( ifreq )
  ENDDO
  CALL mp_sum( dtemp, intra_bgrp_comm ) 
  !
  DO ifreq = 1, ifr%nloc 
     glob_ifreq = ifr%l2g(ifreq)
     IF( glob_ifreq == 1 ) THEN
        imfreq_list_integrate(1,ifreq) = 0._DP
     ELSE
        imfreq_list_integrate(1,ifreq) = ( dtemp(glob_ifreq) + dtemp(glob_ifreq-1) ) * 0.5_DP
     ENDIF
     !
     IF( glob_ifreq == n_imfreq ) THEN
        imfreq_list_integrate(2,ifreq) = dtemp(n_imfreq)
     ELSE
        imfreq_list_integrate(2,ifreq) = ( dtemp(glob_ifreq) + dtemp(glob_ifreq+1) ) * 0.5_DP
     ENDIF
     !
  ENDDO
  !
  DEALLOCATE( dtemp )
  !
  ! TEMP
  ! 
130
131
  ALLOCATE( d_body1_ifr( aband%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
  ALLOCATE( z_body_rfr( aband%nloc, rfr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
Marco Govoni's avatar
Marco Govoni committed
132
  IF( l_enable_lanczos ) THEN
133
134
     ALLOCATE( d_body2_ifr( n_lanczos, pert%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
     ALLOCATE( d_diago( n_lanczos, pert%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
Marco Govoni's avatar
Marco Govoni committed
135
136
137
138
139
140
141
142
143
144
145
146
147
  ENDIF
  !
  d_body1_ifr = 0._DP
  z_body_rfr = 0._DP
  IF( l_enable_lanczos ) THEN
     d_body2_ifr = 0._DP
     d_diago = 0._DP
  ENDIF
  !
  ! d_body1_ifr, d_body2_ifr, z_diago_rfr, d_diago
  !
  CALL io_push_title("Collecting results from W and G")
  !
148
  barra_load = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
Marco Govoni's avatar
Marco Govoni committed
149
150
151
152
  CALL start_bar_type( barra, 'coll_gw', barra_load )
  !
  ! LOOP 
  !
153
  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
Marco Govoni's avatar
Marco Govoni committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
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
201
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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
     !
     nbndval = nbnd_occ(iks)
     !
     DO ib = qp_bandrange(1), qp_bandrange(2)
        !
        IF(ALLOCATED(overlap)) DEALLOCATE(overlap) 
        ALLOCATE(overlap(pert%nglob, nbnd ) )
        !
        CALL readin_overlap( 'g', iks_l2g(iks), ib, overlap, pert%nglob, nbnd )
        !
        ! ------
        ! d_body1_ifr
        ! ------
        !
        ALLOCATE( dtemp(nbnd) )
        !
        DO ifreq = 1, ifr%nloc 
           !
           dtemp = 0._DP
           !
           DO im = 1, nbnd 
              !
              DO glob_jp = 1, n_pdep_eigen_to_use
                 DO ip = 1, pert%nloc
                    glob_ip = pert%l2g(ip)
                    dtemp( im ) = dtemp( im ) + overlap(glob_jp,im)*overlap(glob_ip,im)*d_epsm1_ifr(glob_jp,ip,ifreq)
                 ENDDO
              ENDDO 
              !
           ENDDO ! im
           !
           CALL mp_sum( dtemp, inter_image_comm ) 
           !
           DO im = 1, aband%nloc
              glob_im = aband%l2g(im)
              d_body1_ifr(im,ifreq,ib,iks) = dtemp(glob_im)
           ENDDO
           !
        ENDDO ! ifreq
        !
        DEALLOCATE( dtemp ) 
        !
        ! -----
        ! z_body_rfr
        ! -----
        !
        ALLOCATE( ztemp(nbnd) )
        !
        DO ifreq = 1, rfr%nloc 
           !
           ztemp = 0._DP
           !
           DO im = 1, nbnd 
              !
              DO glob_jp = 1, n_pdep_eigen_to_use
                 DO ip = 1, pert%nloc
                    glob_ip = pert%l2g(ip)
                    ztemp( im ) = ztemp( im ) + overlap(glob_jp,im)*overlap(glob_ip,im)*z_epsm1_rfr(glob_jp,ip,ifreq) 
                 ENDDO
              ENDDO 
              !
           ENDDO ! im
           !
           CALL mp_sum( ztemp, inter_image_comm ) 
           !
           DO im = 1, aband%nloc
              glob_im = aband%l2g(im)
              z_body_rfr(im,ifreq,ib,iks) = ztemp(glob_im)
           ENDDO
           !
        ENDDO ! ifreq
        !
        DEALLOCATE( ztemp ) 
        !
        ! -----------------------------
        ! LANCZOS part : d_diago, d_body2_ifr
        ! -----------------------------
        !
        IF( l_enable_lanczos ) THEN
           !
           ALLOCATE( braket(pert%nglob,n_lanczos,pert%nloc) )
           ALLOCATE( diago(n_lanczos,pert%nloc) )
           CALL readin_solvegfreq( iks_l2g(iks), ib, diago, braket, pert%nloc, pert%nglob, pert%myoffset )
           !
           DO ip = 1, pert%nloc
              DO il = 1, n_lanczos
                 d_diago(il,ip,ib,iks) = diago(il,ip) 
              ENDDO
           ENDDO
           !
           DO ifreq = 1,ifr%nloc
              DO ip = 1, pert%nloc
                 DO il = 1, n_lanczos
                    DO glob_jp = 1, pert%nglob
                       d_body2_ifr(il,ip,ifreq,ib,iks) = d_body2_ifr(il,ip,ifreq,ib,iks) + &
                       & braket(glob_jp,il,ip)*d_epsm1_ifr(glob_jp,ip,ifreq)
                    ENDDO
                 ENDDO
              ENDDO
           ENDDO
           !
           DEALLOCATE(diago)
           DEALLOCATE(braket)
           !
        ENDIF
        !
        CALL update_bar_type( barra, 'coll_gw', 1 )
        !
     ENDDO ! ibnd
     !
  ENDDO ! iks
  !
  CALL stop_bar_type( barra, 'coll_gw' )
  !
  !
  DEALLOCATE( d_epsm1_ifr )
  DEALLOCATE( z_epsm1_rfr )
  !
  ! Get Sigma_X
  !
274
275
  ALLOCATE( sigma_hf (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
  CALL readin_solvehf( sigma_hf(qp_bandrange(1),1), qp_bandrange(2)-qp_bandrange(1)+1, k_grid%nps )
Marco Govoni's avatar
Marco Govoni committed
276
277
278
279
280
281
282
  !
  ! For CORR
  !
  IF( l_secant ) THEN 
     !
     CALL io_push_title("(Q)uasiparticle energies")
     !
283
284
285
286
287
288
289
     ALLOCATE( en(qp_bandrange(1):qp_bandrange(2),k_grid%nps,2) ) 
     ALLOCATE( sc(qp_bandrange(1):qp_bandrange(2),k_grid%nps,2) ) 
     ALLOCATE( l_conv(qp_bandrange(1):qp_bandrange(2),k_grid%nps) ) 
     ALLOCATE( sigma_cor_in (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
     ALLOCATE( sigma_cor_out(qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
     ALLOCATE( z_in (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
     ALLOCATE( qp_energy (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
Marco Govoni's avatar
Marco Govoni committed
290
291
292
     !
     ! 1st step of secant solver : E_KS - 0.5 * eshift
     !
293
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
294
295
296
297
298
299
        DO ib = qp_bandrange(1), qp_bandrange(2)
           en(ib,iks,1) = et(ib,iks) - eshift*0.5_DP
        ENDDO
     ENDDO
     CALL calc_corr_gamma( sc(:,:,1), en(:,:,1), .TRUE.)  
     !
300
     !DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
301
302
303
304
305
306
307
     !   DO ib = qp_bandrange(1), qp_bandrange(2)
     !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,1) * rytoev, sc(ib,iks,1) * rytoev
     !   ENDDO
     !ENDDO
     !
     ! 1st step of secant solver : E_KS + 0.5 * eshift
     !
308
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
309
310
311
312
313
        DO ib = qp_bandrange(1), qp_bandrange(2)
           en(ib,iks,2) = et(ib,iks) + eshift*0.5_DP
        ENDDO
     ENDDO
     CALL calc_corr_gamma( sc(:,:,2), en(:,:,2), .TRUE.) 
314
     !DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
315
316
317
318
319
320
321
     !   DO ib = qp_bandrange(1), qp_bandrange(2)
     !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,2) * rytoev, sc(ib,iks,2) * rytoev
     !   ENDDO
     !ENDDO
     !
     ! Stage sigma_corr_in
     !
322
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
323
324
325
326
327
328
329
        DO ib = qp_bandrange(1), qp_bandrange(2)
           sigma_cor_in(ib,iks) = ( sc(ib,iks,2) + sc(ib,iks,1) ) * 0.5_DP 
        ENDDO
     ENDDO
     !
     ! Stage z_in
     !
330
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
331
332
333
334
335
336
337
338
        DO ib = qp_bandrange(1), qp_bandrange(2)
           z_in(ib,iks) = 1._DP / ( 1._DP-REAL( sc(ib,iks,2)-sc(ib,iks,1), KIND=DP ) / eshift  ) 
        ENDDO
     ENDDO
     !
     ! en 1 = EKS, sc 1 = sigma_corr_in
     ! en 2 = EKS + Z * ( sigmax - vxc + sigmacorrin)
     !
339
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
340
341
342
343
344
345
346
347
348
349
350
        DO ib = qp_bandrange(1), qp_bandrange(2)
           en(ib,iks,1) = et(ib,iks) 
           sc(ib,iks,1) = sigma_cor_in(ib,iks)
           en(ib,iks,2) = et(ib,iks) + z_in(ib,iks) * ( sigma_hf(ib,iks) + REAL( sigma_cor_in(ib,iks),KIND=DP) ) 
        ENDDO
     ENDDO
     sigma_z           = z_in
     sigma_eqplin(:,:) = en(:,:,2)
     sigma_sc_eks      = sigma_cor_in
     CALL output_eqp_report(0,en(:,:,1),en(:,:,2),sigma_cor_in(:,:))
     !
Govoni's avatar
Govoni committed
351
     ! nth step of the secant solver
Marco Govoni's avatar
Marco Govoni committed
352
353
     !
     l_conv = .FALSE.
354
     notconv = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
Marco Govoni's avatar
Marco Govoni committed
355
356
357
358
359
360
     DO ifixed = 1, n_secant_maxiter
        ! 
        WRITE(cifixed,"(i6.6)") ifixed
        !CALL io_push_title("Iter : "//cifixed)
        !
        CALL calc_corr_gamma( sc(:,:,2), en(:,:,2), .TRUE.)  
361
        !DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
362
363
364
365
366
367
368
        !   DO ib = qp_bandrange(1), qp_bandrange(2)
        !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,2) * rytoev, sc(ib,iks,2) * rytoev
        !   ENDDO
        !ENDDO
        !
        ! Resulting new energy:
        !
369
        DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
370
371
372
373
374
375
376
377
378
379
380
           DO ib = qp_bandrange(1), qp_bandrange(2)
               IF( .NOT. l_conv(ib,iks) ) THEN
                  qp_energy(ib,iks) = en(ib,iks,2) + &
                         & ( et(ib,iks) + sigma_hf(ib,iks) + REAL(sc(ib,iks,2),KIND=DP) - en(ib,iks,2) ) / &
                         & ( 1._DP - REAL( sc(ib,iks,2) - sc(ib,iks,1), KIND=DP ) / ( en(ib,iks,2) - en(ib,iks,1) ) )
               ENDIF
           ENDDO
        ENDDO
        !
        ! Estimate l_conv
        !
381
        DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
382
383
384
385
386
           DO ib = qp_bandrange(1), qp_bandrange(2)
              l_conv(ib,iks) = ( ABS( qp_energy(ib,iks) - en(ib,iks,2) ) < trev_secant ) 
           ENDDO
        ENDDO
        !
387
388
389
390
391
       !WRITE(stdout,"(5X)")  
       !CALL io_push_bar()
       !WRITE(stdout,"(5X,'Iter: ',i6.6)") ifixed
       !WRITE(stdout,"(5X,a,1X,a,1X,a,1X,a)") 'K     ', 'B     ', 'QP energ. [eV]', 'conv'
       !CALL io_push_bar()
392
       !DO iks = 1, k_grid%nps
393
394
395
396
397
       !   DO ib = qp_bandrange(1), qp_bandrange(2)
       !      WRITE(stdout,"(5X,i6.6,1X,i6.6,1X,1f14.6,4X,l)") iks, ib, qp_energy(ib,iks) * rytoev, l_conv(ib,iks) 
       !   ENDDO
       !ENDDO
       !CALL io_push_bar()
Govoni's avatar
Govoni committed
398
        !
Marco Govoni's avatar
Marco Govoni committed
399
400
401
        ! Count the number of notconverged QP energies
        !
        notconv = 0 
402
        DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
           DO ib = qp_bandrange(1), qp_bandrange(2)
              IF( .NOT.l_conv(ib,iks) ) notconv = notconv + 1
           ENDDO
        ENDDO
        !
        IF( ifixed == 1 ) THEN
           sigma_sc_eqplin(:,:) = sc(:,:,2)  
        ENDIF 
        !
        sigma_eqpsec = qp_energy 
        sigma_sc_eqpsec = sc(:,:,2)
        sigma_diff(:,:) = qp_energy(:,:) - en(:,:,2)
        !
        IF( notconv == 0 ) THEN 
           ! 
           CALL io_push_title("CONVERGENCE ACHIEVED !!!")
419
           CALL output_eqp_report(-1,en(:,:,2),qp_energy,sc(:,:,2))
Marco Govoni's avatar
Marco Govoni committed
420
421
422
           EXIT
           !
        ELSE
Govoni's avatar
Govoni committed
423
424
425
           !CALL io_push_bar()
           !WRITE(stdout,'(5x,"Number of unconverged QP. = ",i12)') notconv
           !CALL io_push_bar()
Marco Govoni's avatar
Marco Govoni committed
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
           CALL output_eqp_report(ifixed,en(:,:,2),qp_energy,sc(:,:,2))
           !
        ENDIF
        !
        ! Iterate
        !
        en(:,:,1) = en(:,:,2)    
        sc(:,:,1) = sc(:,:,2) 
        en(:,:,2) = qp_energy(:,:)   
        !
     ENDDO
     !
     sigma_cor_out(:,:) = sc(:,:,2)
     DEALLOCATE( en, sc, l_conv )
     !
     IF( notconv .NE. 0 ) THEN 
        CALL io_push_title("CONVERGENCE **NOT** ACHIEVED !!!")
     ENDIF
     !
     ! Output it per k-point
     !
447
448
     ALLOCATE(out_tab(qp_bandrange(2)-qp_bandrange(1)+1,7))
     !
449
     DO iks=1,k_grid%nps
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
        DO ib = qp_bandrange(1), qp_bandrange(2)
           out_tab( ib - qp_bandrange(1) + 1, 1) = REAL( ib, KIND=DP) 
           out_tab( ib - qp_bandrange(1) + 1, 2) = et(ib,iks) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 3) = (et(ib,iks)+sigma_hf(ib,iks)) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 4) = qp_energy(ib,iks) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 5) = (qp_energy(ib,iks) - et(ib,iks) ) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 6) = REAL(  sigma_cor_out(ib,iks), KIND=DP ) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 7) = AIMAG( sigma_cor_out(ib,iks) ) * rytoev
        ENDDO
        WRITE(myglobk,'(i5.5)') iks_l2g(iks)
        !
        CALL serial_table_output(mpime==root,4000,'eqp_K'//myglobk,out_tab,&
           & qp_bandrange(2)-qp_bandrange(1)+1,7,&
           & (/'      band','    E0[eV]','   EHF[eV]','   Eqp[eV]','Eqp-E0[eV]','Sc_Eqp[eV]',' Width[eV]'/))
     ENDDO
     !
     DEALLOCATE( out_tab )  
Marco Govoni's avatar
Marco Govoni committed
467
468
469
470
471
472
473
474
475
476
477
478
479
     DEALLOCATE( sigma_cor_in )
     DEALLOCATE( sigma_cor_out )
     DEALLOCATE( z_in )
     DEALLOCATE( qp_energy )
     !
  ENDIF
  !
  IF( l_generate_plot ) THEN
     !
     CALL io_push_title("(P)lotting the QP corrections...")
     !
     CALL start_bar_type( barra, 'qplot', n_spectralf )
     !
480
481
     ALLOCATE( en(qp_bandrange(1):qp_bandrange(2),k_grid%nps,1) ) 
     ALLOCATE( sc(qp_bandrange(1):qp_bandrange(2),k_grid%nps,1) )
Govoni's avatar
Govoni committed
482
483
484
485
     !ALLOCATE( un(qp_bandrange(1):qp_bandrange(2),nks) )
     !
   ! IF( mpime == 0 ) THEN 
   !    k = 0 
486
   !    DO iks=1,k_grid%nps
Govoni's avatar
Govoni committed
487
488
489
490
491
492
493
494
495
496
497
498
499
500
   !       WRITE(iks_label,"(i5.5)") iks
   !       DO ib = qp_bandrange(1), qp_bandrange(2)
   !          WRITE(ib_label,"(i5.5)") ib
   !          fname = "o-eqp_K"//TRIM(iks_label)//"_B"//TRIM(ib_label)//".tab"
   !          k = k + 1
   !          un(ib,iks) = 10000+k
   !          OPEN( UNIT=un(ib,iks), FILE=TRIM(fname))
   !          WRITE( un(ib,iks),"(a)") '#        E[eV]         E-EKS  Re{S(E)}-Vxc      Im{S(E)}          A(E)'
   !       ENDDO  
   !       fname = "o-eqp_K"//TRIM(iks_label)//"_summed.tab"
   !       OPEN( UNIT=10000-iks, FILE=TRIM(fname))
   !       WRITE( 10000-iks,"(a)") '#        E[eV]          A(E)'
   !    ENDDO
   ! ENDIF 
Marco Govoni's avatar
Marco Govoni committed
501
502
503
504
505
506
507
508
509
     !
     DO glob_ifreq = 1, n_spectralf
        sigma_freq(glob_ifreq) = (ecut_spectralf(2)-ecut_spectralf(1)) / REAL(n_spectralf-1,KIND=DP) * REAL(glob_ifreq-1,KIND=DP) &
                               & + ecut_spectralf(1)
     ENDDO 
     !
     DO glob_ifreq = 1, n_spectralf
        en = (ecut_spectralf(2)-ecut_spectralf(1)) / REAL(n_spectralf-1,KIND=DP) * REAL(glob_ifreq-1,KIND=DP) + ecut_spectralf(1)
        CALL calc_corr_gamma( sc(:,:,1), en(:,:,1), .FALSE.) 
510
        DO iks=1,k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
511
512
513
514
           DO ib = qp_bandrange(1), qp_bandrange(2)
              sigma_spectralf(glob_ifreq,ib,iks) = sc(ib,iks,1) 
           ENDDO
        ENDDO
Govoni's avatar
Govoni committed
515
       !IF( mpime == 0 ) THEN 
516
       !   DO iks=1,k_grid%nps
Govoni's avatar
Govoni committed
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
       !      summed_sf=0._DP
       !      DO ib = qp_bandrange(1), qp_bandrange(2)
       !         WRITE( un(ib,iks),"(5f14.6)") en(ib,iks,1)*rytoev, &
       !         & (en(ib,iks,1)-et(ib,iks))*rytoev, &
       !         & (DBLE(sc(ib,iks,1))+sigma_hf(ib,iks))*rytoev, &
       !         & AIMAG(sc(ib,iks,1))*rytoev, &
       !         & ABS(AIMAG(sc(ib,iks,1)))/((en(ib,iks,1)-et(ib,iks)-sigma_hf(ib,iks)-DBLE(sc(ib,iks,1)))**2+&
       !         &AIMAG(sc(ib,iks,1))**2)/pi/rytoev
       !         IF( ib <= nbnd_occ(iks) ) THEN 
       !            summed_sf=summed_sf+ABS(AIMAG(sc(ib,iks,1)))/((en(ib,iks,1)-et(ib,iks)-sigma_hf(ib,iks)-DBLE(sc(ib,iks,1)))**2&
       !            &+AIMAG(sc(ib,iks,1))**2)/pi
       !         ENDIF
       !      ENDDO
       !      WRITE( 10000-iks,"(2f14.6)") en(qp_bandrange(1),iks,1)*rytoev,summed_sf/rytoev
       !   ENDDO
       !ENDIF
Marco Govoni's avatar
Marco Govoni committed
533
534
535
        CALL update_bar_type(barra,'qplot',1)
     ENDDO
     !
Govoni's avatar
Govoni committed
536
    !IF( mpime == 0 ) THEN 
537
    !   DO iks=1,k_grid%nps
Govoni's avatar
Govoni committed
538
539
540
541
542
543
    !      DO ib = qp_bandrange(1), qp_bandrange(2)
    !         CLOSE( un(ib,iks) )
    !      ENDDO 
    !      CLOSE(10000-iks)
    !   ENDDO
    !ENDIF 
Marco Govoni's avatar
Marco Govoni committed
544
545
546
547
548
     !
     CALL stop_bar_type(barra,'qplot')
     !
     DEALLOCATE( en )
     DEALLOCATE( sc )  
Govoni's avatar
Govoni committed
549
     !DEALLOCATE( un )  
Marco Govoni's avatar
Marco Govoni committed
550
     !
Govoni's avatar
Govoni committed
551
     !CALL io_push_title('Done, take a look at the o-eqp_K*_B*.tab file(s) .')
Marco Govoni's avatar
Marco Govoni committed
552
553
554
555
556
     ! 
  ENDIF 
  !
  DEALLOCATE( sigma_hf )
  !
Govoni's avatar
Govoni committed
557
  CALL wfreq_db_write( )
Marco Govoni's avatar
Marco Govoni committed
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
  !
  CALL stop_clock( "solve_qp" )
  !
  !
END SUBROUTINE
!
!-----------------------------------------------------------------------
SUBROUTINE solve_qp_k(l_secant,l_generate_plot)
  !-----------------------------------------------------------------------
  !
  ! ... This subroutine solves the DBS problem for GAMMA, at non-zero freqeuncies. 
  ! ... Perturbations are distributed according to the POT mpi_communicator
  !
  USE kinds,                ONLY : DP 
  USE westcom,              ONLY : n_pdep_eigen_to_use,n_lanczos,qp_bandrange,iks_l2g,imfreq_list_integrate,&
                                 & n_secant_maxiter,nbnd_occ,trev_secant,l_enable_lanczos,imfreq_list,n_imfreq,&
574
575
                                 & z_epsm1_ifr_q,z_epsm1_rfr_q,l_macropol,n_spectralf,ecut_spectralf, &
                                 & z_head_ifr,z_body1_ifr_q,z_body2_ifr_q,d_diago_q,z_head_rfr,z_body_rfr_q, &
Marco Govoni's avatar
Marco Govoni committed
576
577
578
579
580
581
582
583
584
585
586
587
588
                                 & sigma_z,sigma_eqplin,sigma_eqpsec,sigma_sc_eks,sigma_sc_eqplin,sigma_sc_eqpsec,sigma_diff, &
                                 & sigma_spectralf,sigma_freq
  USE mp_global,            ONLY : inter_image_comm,intra_bgrp_comm
  USE mp_world,             ONLY : mpime,root
  USE mp,                   ONLY : mp_barrier,mp_sum
  USE io_global,            ONLY : stdout, ionode
  USE pwcom,                ONLY : et,nks,current_spin,isk,xk,nbnd,lsda,g2kin,current_k
  USE io_push,              ONLY : io_push_title,io_push_bar
  USE constants,            ONLY : rytoev,pi
  USE west_io,              ONLY : serial_table_output
  USE distribution_center,  ONLY : pert,ifr,rfr,aband
  USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
  USE wfreq_io,             ONLY : readin_overlap,readin_solvegfreq,readin_solvehf
Govoni's avatar
Govoni committed
589
  USE wfreq_db,             ONLY : wfreq_db_write
590
  USE types_bz_grid,        ONLY : k_grid, q_grid
Marco Govoni's avatar
Marco Govoni committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  LOGICAL,INTENT(IN) :: l_secant,l_generate_plot 
  !
  ! Workspace
  !
  REAL(DP),ALLOCATABLE :: sigma_hf(:,:)
  COMPLEX(DP),ALLOCATABLE :: sigma_cor_in (:,:)
  COMPLEX(DP),ALLOCATABLE :: sigma_cor_out(:,:)
  REAL(DP),ALLOCATABLE :: z_in(:,:)
  REAL(DP),ALLOCATABLE :: qp_energy(:,:)
  COMPLEX(DP),ALLOCATABLE :: sc(:,:,:) 
  REAL(DP),ALLOCATABLE :: en(:,:,:) 
  LOGICAL,ALLOCATABLE :: l_conv(:,:)
  REAL(DP),PARAMETER   :: eshift = 0.007349862_DP ! = 0.1 eV
609
610
  INTEGER :: k, ib, iks, ik, ikks, ikk, iq, ifixed, ip, glob_ip, ifreq, il, im, glob_im, glob_jp, glob_ifreq
  REAL(DP) :: q(3), g0(3)
Marco Govoni's avatar
Marco Govoni committed
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
  REAL(DP),ALLOCATABLE :: out_tab(:,:)
  CHARACTER(LEN=5) :: myglobk 
  CHARACTER(LEN=6) :: cifixed
  INTEGER :: notconv
  INTEGER,EXTERNAL :: get_nbndval
  INTEGER :: nbndval 
  REAL(DP),ALLOCATABLE :: dtemp(:)
  COMPLEX(DP),ALLOCATABLE :: ztemp(:)
  REAL(DP),ALLOCATABLE :: diago(:,:)
  COMPLEX(DP),ALLOCATABLE :: braket(:,:,:)
  COMPLEX(DP),ALLOCATABLE :: overlap(:,:)
  TYPE(bar_type) :: barra
  INTEGER :: barra_load
  CHARACTER(LEN=126) :: fname
  CHARACTER(LEN=5) :: ib_label, iks_label
  INTEGER,ALLOCATABLE :: un(:,:)
  REAL(DP) :: summed_sf
628
  !TYPE(bz_grid) :: k1_grid, q_grid_aux
Marco Govoni's avatar
Marco Govoni committed
629
630
631
632
633
634
  ! 
  CALL start_clock('solve_qp')
  !
  ALLOCATE( imfreq_list_integrate( 2, ifr%nloc ) )
  ALLOCATE( dtemp( n_imfreq ) ) 
  !
635
636
637
638
  !k1_grid = bz_grid()
  !CALL k1_grid%init('K')
  !q_grid_aux = bz_grid()
  !CALL q_grid_aux%init_q( k_grid, k1_grid )
Marco Govoni's avatar
Marco Govoni committed
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
  !
  dtemp = 0._DP
  DO ifreq = 1, ifr%nloc 
     glob_ifreq = ifr%l2g(ifreq)
     dtemp( glob_ifreq ) = imfreq_list( ifreq )
  ENDDO
  CALL mp_sum( dtemp, intra_bgrp_comm ) 
  !
  DO ifreq = 1, ifr%nloc 
     glob_ifreq = ifr%l2g(ifreq)
     IF( glob_ifreq == 1 ) THEN
        imfreq_list_integrate(1,ifreq) = 0._DP
     ELSE
        imfreq_list_integrate(1,ifreq) = ( dtemp(glob_ifreq) + dtemp(glob_ifreq-1) ) * 0.5_DP
     ENDIF
     !
     IF( glob_ifreq == n_imfreq ) THEN
        imfreq_list_integrate(2,ifreq) = dtemp(n_imfreq)
     ELSE
        imfreq_list_integrate(2,ifreq) = ( dtemp(glob_ifreq) + dtemp(glob_ifreq+1) ) * 0.5_DP
     ENDIF
     !
  ENDDO
  !
  DEALLOCATE( dtemp )
  !
  ! TEMP
  ! 
667
668
  ALLOCATE( z_body1_ifr_q( aband%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps, q_grid%nps ) )
  ALLOCATE( z_body_rfr_q( aband%nloc, rfr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps, q_grid%nps ) )
Marco Govoni's avatar
Marco Govoni committed
669
  IF( l_enable_lanczos ) THEN
670
671
     ALLOCATE( z_body2_ifr_q( n_lanczos, pert%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps, q_grid%nps ) )
     ALLOCATE( d_diago_q( n_lanczos, pert%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps, q_grid%nps ) )
Marco Govoni's avatar
Marco Govoni committed
672
673
  ENDIF
  !
674
675
  z_body1_ifr_q = 0._DP
  z_body_rfr_q = 0._DP
Marco Govoni's avatar
Marco Govoni committed
676
  IF( l_enable_lanczos ) THEN
677
678
     z_body2_ifr_q = 0._DP
     d_diago_q = 0._DP
Marco Govoni's avatar
Marco Govoni committed
679
680
  ENDIF
  !
681
  ! d_body1_ifr_q, d_body2_ifr_q, z_diago_rfr_q, d_diago_q
Marco Govoni's avatar
Marco Govoni committed
682
683
684
  !
  CALL io_push_title("Collecting results from W and G")
  !
685
  barra_load = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 ) * q_grid%nps
Marco Govoni's avatar
Marco Govoni committed
686
687
688
689
  CALL start_bar_type( barra, 'coll_gw', barra_load )
  !
  ! LOOP 
  !
690
691
692
  ! ... Outer k-point loop (wfc matrix element): iks
  ! ... Inner k-point loop (wfc summed over k'): ikks
  ! ... BEWARE: iks and ikks are switched w.r.t. solve_gfreq_k
693
694
  !
  DO iks = 1, k_grid%nps   ! KPOINT-SPIN (MATRIX ELEMENT)
695
     !
696
     ik = k_grid%ip(iks)
Marco Govoni's avatar
Marco Govoni committed
697
698
699
700
701
     !
     nbndval = nbnd_occ(iks)
     !
     DO ib = qp_bandrange(1), qp_bandrange(2)
        !
702
        DO ikks = 1, k_grid%nps   ! KPOINT-SPIN (INTEGRAL OVER K')
703
           !
704
           ikk = k_grid%ip(ikks)
Marco Govoni's avatar
Marco Govoni committed
705
           !
706
           CALL q_grid%find( k_grid%p_cart(:,ik) - k_grid%p_cart(:,ikk), 1, 'cart', iq, g0 )
Marco Govoni's avatar
Marco Govoni committed
707
           !
708
709
           IF(ALLOCATED(overlap)) DEALLOCATE(overlap) 
           ALLOCATE(overlap(pert%nglob, nbnd ) )
Marco Govoni's avatar
Marco Govoni committed
710
           !
711
           CALL readin_overlap( 'g', iks_l2g(iks), iks_l2g(ikks), ib, overlap, pert%nglob, nbnd )
Marco Govoni's avatar
Marco Govoni committed
712
           !
713
714
715
           ! ------
           ! z_body1_ifr_q
           ! ------
Marco Govoni's avatar
Marco Govoni committed
716
           !
717
           ALLOCATE( ztemp(nbnd) )
Marco Govoni's avatar
Marco Govoni committed
718
           !
719
           DO ifreq = 1, ifr%nloc 
Marco Govoni's avatar
Marco Govoni committed
720
              !
721
              ztemp = 0._DP
Marco Govoni's avatar
Marco Govoni committed
722
              !
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
              DO im = 1, nbnd 
                 !
                 DO glob_jp = 1, n_pdep_eigen_to_use
                    DO ip = 1, pert%nloc
                       glob_ip = pert%l2g(ip)
                       ztemp( im ) = ztemp( im ) + DCONJG(overlap(glob_jp,im))*overlap(glob_ip,im)* &
                                                   & z_epsm1_ifr_q(glob_jp,ip,ifreq,iq)
                    ENDDO
                 ENDDO 
                 !
              ENDDO ! im
              !
              CALL mp_sum( ztemp, inter_image_comm ) 
              !
              DO im = 1, aband%nloc
                 glob_im = aband%l2g(im)
                 z_body1_ifr_q(im,ifreq,ib,iks,iq) = ztemp(glob_im)
              ENDDO
              !
              !
           ENDDO ! ifreq
Marco Govoni's avatar
Marco Govoni committed
744
           !
745
           DEALLOCATE( ztemp ) 
Marco Govoni's avatar
Marco Govoni committed
746
           !
747
748
749
           ! -----
           ! z_body_rfr_q
           ! -----
Marco Govoni's avatar
Marco Govoni committed
750
           !
751
           ALLOCATE( ztemp(nbnd) )
Marco Govoni's avatar
Marco Govoni committed
752
           !
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
           DO ifreq = 1, rfr%nloc 
              !
              ztemp = 0._DP
              !
              DO im = 1, nbnd 
                 !
                 DO glob_jp = 1, n_pdep_eigen_to_use
                    DO ip = 1, pert%nloc
                       glob_ip = pert%l2g(ip)
                       ztemp( im ) = ztemp( im ) + DCONJG(overlap(glob_jp,im))*overlap(glob_ip,im) * &
                       & z_epsm1_rfr_q(glob_jp,ip,ifreq,iq) 
                    ENDDO
                 ENDDO 
                 !
              ENDDO ! im
              !
              CALL mp_sum( ztemp, inter_image_comm ) 
              !
              DO im = 1, aband%nloc
                 glob_im = aband%l2g(im)
                 z_body_rfr_q(im,ifreq,ib,iks,iq) = ztemp(glob_im)
Marco Govoni's avatar
Marco Govoni committed
774
              ENDDO
775
776
              !
           ENDDO ! ifreq
Marco Govoni's avatar
Marco Govoni committed
777
           !
778
779
780
781
782
783
784
785
786
787
788
789
           DEALLOCATE( ztemp ) 
           !
           ! -----------------------------
           ! LANCZOS part : d_diago_q, z_body2_ifr_q
           ! -----------------------------
           !
           IF( l_enable_lanczos ) THEN
              !
              ALLOCATE( braket(pert%nglob,n_lanczos,pert%nloc) )
              ALLOCATE( diago(n_lanczos,pert%nloc) )
              CALL readin_solvegfreq( iks_l2g(iks), iks_l2g(ikks), ib, diago, braket, pert%nloc, pert%nglob, pert%myoffset )
              !
Marco Govoni's avatar
Marco Govoni committed
790
791
              DO ip = 1, pert%nloc
                 DO il = 1, n_lanczos
792
793
794
795
796
797
798
799
800
801
802
                    d_diago_q(il,ip,ib,iks,iq) = diago(il,ip) 
                 ENDDO
              ENDDO
              !
              DO ifreq = 1,ifr%nloc
                 DO ip = 1, pert%nloc
                    DO il = 1, n_lanczos
                       DO glob_jp = 1, pert%nglob
                          z_body2_ifr_q(il,ip,ifreq,ib,iks,iq) = z_body2_ifr_q(il,ip,ifreq,ib,iks,iq) + &
                          & braket(glob_jp,il,ip)*z_epsm1_ifr_q(glob_jp,ip,ifreq,iq)
                       ENDDO
Marco Govoni's avatar
Marco Govoni committed
803
804
805
                    ENDDO
                 ENDDO
              ENDDO
806
807
808
809
810
              !
              DEALLOCATE(diago)
              DEALLOCATE(braket)
              !
           ENDIF
Marco Govoni's avatar
Marco Govoni committed
811
           !
812
           CALL update_bar_type( barra, 'coll_gw', 1 )
Marco Govoni's avatar
Marco Govoni committed
813
           !
814
        ENDDO ! ibnd
Marco Govoni's avatar
Marco Govoni committed
815
        !
816
     ENDDO ! ikks
Marco Govoni's avatar
Marco Govoni committed
817
818
819
820
821
     !
  ENDDO ! iks
  !
  CALL stop_bar_type( barra, 'coll_gw' )
  !
822
823
  DEALLOCATE( z_epsm1_ifr_q )
  DEALLOCATE( z_epsm1_rfr_q )
Marco Govoni's avatar
Marco Govoni committed
824
825
826
  !
  ! Get Sigma_X
  !
827
828
  ALLOCATE( sigma_hf (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
  CALL readin_solvehf( sigma_hf(qp_bandrange(1),1), qp_bandrange(2)-qp_bandrange(1)+1, k_grid%nps )
Marco Govoni's avatar
Marco Govoni committed
829
830
831
832
833
834
835
  !
  ! For CORR
  !
  IF( l_secant ) THEN 
     !
     CALL io_push_title("(Q)uasiparticle energies")
     !
836
837
838
839
840
841
842
     ALLOCATE( en(qp_bandrange(1):qp_bandrange(2),k_grid%nps,2) ) 
     ALLOCATE( sc(qp_bandrange(1):qp_bandrange(2),k_grid%nps,2) ) 
     ALLOCATE( l_conv(qp_bandrange(1):qp_bandrange(2),k_grid%nps) ) 
     ALLOCATE( sigma_cor_in (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
     ALLOCATE( sigma_cor_out(qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
     ALLOCATE( z_in (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
     ALLOCATE( qp_energy (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
Marco Govoni's avatar
Marco Govoni committed
843
844
845
     !
     ! 1st step of secant solver : E_KS - 0.5 * eshift
     !
846
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
847
848
849
850
851
852
853
854
855
856
857
858
859
860
        DO ib = qp_bandrange(1), qp_bandrange(2)
           en(ib,iks,1) = et(ib,iks) - eshift*0.5_DP
        ENDDO
     ENDDO
     CALL calc_corr_k( sc(:,:,1), en(:,:,1), .TRUE.)  
     !
     !DO iks = 1, nks
     !   DO ib = qp_bandrange(1), qp_bandrange(2)
     !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,1) * rytoev, sc(ib,iks,1) * rytoev
     !   ENDDO
     !ENDDO
     !
     ! 1st step of secant solver : E_KS + 0.5 * eshift
     !
861
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
862
863
864
865
866
867
868
869
870
871
872
873
874
        DO ib = qp_bandrange(1), qp_bandrange(2)
           en(ib,iks,2) = et(ib,iks) + eshift*0.5_DP
        ENDDO
     ENDDO
     CALL calc_corr_k( sc(:,:,2), en(:,:,2), .TRUE.) 
     !DO iks = 1, nks
     !   DO ib = qp_bandrange(1), qp_bandrange(2)
     !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,2) * rytoev, sc(ib,iks,2) * rytoev
     !   ENDDO
     !ENDDO
     !
     ! Stage sigma_corr_in
     !
875
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
876
877
878
879
880
881
882
        DO ib = qp_bandrange(1), qp_bandrange(2)
           sigma_cor_in(ib,iks) = ( sc(ib,iks,2) + sc(ib,iks,1) ) * 0.5_DP 
        ENDDO
     ENDDO
     !
     ! Stage z_in
     !
883
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
884
885
886
887
888
889
890
891
        DO ib = qp_bandrange(1), qp_bandrange(2)
           z_in(ib,iks) = 1._DP / ( 1._DP-REAL( sc(ib,iks,2)-sc(ib,iks,1), KIND=DP ) / eshift  ) 
        ENDDO
     ENDDO
     !
     ! en 1 = EKS, sc 1 = sigma_corr_in
     ! en 2 = EKS + Z * ( sigmax - vxc + sigmacorrin)
     !
892
     DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
893
894
895
896
897
898
899
900
901
902
903
        DO ib = qp_bandrange(1), qp_bandrange(2)
           en(ib,iks,1) = et(ib,iks) 
           sc(ib,iks,1) = sigma_cor_in(ib,iks)
           en(ib,iks,2) = et(ib,iks) + z_in(ib,iks) * ( sigma_hf(ib,iks) + REAL( sigma_cor_in(ib,iks),KIND=DP) ) 
        ENDDO
     ENDDO
     sigma_z           = z_in
     sigma_eqplin(:,:) = en(:,:,2)
     sigma_sc_eks      = sigma_cor_in
     CALL output_eqp_report(0,en(:,:,1),en(:,:,2),sigma_cor_in(:,:))
     !
904
905
906
907
908
    !WRITE(stdout,"(5X)")  
    !CALL io_push_bar()
    !WRITE(stdout,"(5X,'Iter: ',i6.6)") 0 
    !WRITE(stdout,"(5X,a,1X,a,1X,a,1X,a)") 'K     ', 'B     ', 'QP energ. [eV]', 'conv'
    !CALL io_push_bar()
909
    !DO iks = 1, k_grid%nps
910
911
912
913
914
    !   DO ib = qp_bandrange(1), qp_bandrange(2)
    !      WRITE(stdout,"(5X,i6.6,1X,i6.6,1X,1f14.6,4X,l)") iks, ib, en(ib,iks,2) * rytoev, .false. 
    !   ENDDO
    !ENDDO
    !CALL io_push_bar()
Marco Govoni's avatar
Marco Govoni committed
915
     !
Govoni's avatar
Govoni committed
916
     ! nth step of the secant solver
Marco Govoni's avatar
Marco Govoni committed
917
918
     !
     l_conv = .FALSE.
919
     notconv =  k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
Marco Govoni's avatar
Marco Govoni committed
920
921
922
923
924
925
     DO ifixed = 1, n_secant_maxiter
        ! 
        WRITE(cifixed,"(i6.6)") ifixed
        !CALL io_push_title("Iter : "//cifixed)
        !
        CALL calc_corr_k( sc(:,:,2), en(:,:,2), .TRUE.)  
926
        !DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
927
928
929
930
931
932
933
        !   DO ib = qp_bandrange(1), qp_bandrange(2)
        !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,2) * rytoev, sc(ib,iks,2) * rytoev
        !   ENDDO
        !ENDDO
        !
        ! Resulting new energy:
        !
934
        DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
935
936
937
938
939
940
941
942
943
944
945
           DO ib = qp_bandrange(1), qp_bandrange(2)
               IF( .NOT. l_conv(ib,iks) ) THEN
                  qp_energy(ib,iks) = en(ib,iks,2) + &
                         & ( et(ib,iks) + sigma_hf(ib,iks) + REAL(sc(ib,iks,2),KIND=DP) - en(ib,iks,2) ) / &
                         & ( 1._DP - REAL( sc(ib,iks,2) - sc(ib,iks,1), KIND=DP ) / ( en(ib,iks,2) - en(ib,iks,1) ) )
               ENDIF
           ENDDO
        ENDDO
        !
        ! Estimate l_conv
        !
946
        DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
947
948
949
950
951
           DO ib = qp_bandrange(1), qp_bandrange(2)
              l_conv(ib,iks) = ( ABS( qp_energy(ib,iks) - en(ib,iks,2) ) < trev_secant ) 
           ENDDO
        ENDDO
        !
952
953
954
955
956
       !WRITE(stdout,"(5X)")  
       !CALL io_push_bar()
       !WRITE(stdout,"(5X,'Iter: ',i6.6)") ifixed
       !WRITE(stdout,"(5X,a,1X,a,1X,a,1X,a)") 'K     ', 'B     ', 'QP energ. [eV]', 'conv'
       !CALL io_push_bar()
957
       !DO iks = 1, k_grid%nps
958
959
960
961
962
       !   DO ib = qp_bandrange(1), qp_bandrange(2)
       !      WRITE(stdout,"(5X,i6.6,1X,i6.6,1X,1f14.6,4X,l)") iks, ib, qp_energy(ib,iks) * rytoev, l_conv(ib,iks) 
       !   ENDDO
       !ENDDO
       !CALL io_push_bar()
963
964
965
       !
       ! Count the number of notconverged QP energies
       !
Marco Govoni's avatar
Marco Govoni committed
966
        notconv = 0 
967
        DO iks = 1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
           DO ib = qp_bandrange(1), qp_bandrange(2)
              IF( .NOT.l_conv(ib,iks) ) notconv = notconv + 1
           ENDDO
        ENDDO
        !
        IF( ifixed == 1 ) THEN
           sigma_sc_eqplin(:,:) = sc(:,:,2)  
        ENDIF 
        !
        sigma_eqpsec = qp_energy 
        sigma_sc_eqpsec = sc(:,:,2)
        sigma_diff(:,:) = qp_energy(:,:) - en(:,:,2)
        !
        IF( notconv == 0 ) THEN 
           ! 
           CALL io_push_title("CONVERGENCE ACHIEVED !!!")
984
           CALL output_eqp_report(-1,en(:,:,2),qp_energy,sc(:,:,2))
Marco Govoni's avatar
Marco Govoni committed
985
986
           EXIT
        ELSE
Govoni's avatar
Govoni committed
987
988
989
           !CALL io_push_bar()
           !WRITE(stdout,'(5x,"Number of unconverged QP. = ",i12)') notconv
           !CALL io_push_bar()
Marco Govoni's avatar
Marco Govoni committed
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
           CALL output_eqp_report(ifixed,en(:,:,2),qp_energy,sc(:,:,2))
        ENDIF
        !
        ! Iterate
        !
        en(:,:,1) = en(:,:,2)    
        sc(:,:,1) = sc(:,:,2) 
        en(:,:,2) = qp_energy(:,:)   
        !
     ENDDO
     !
     sigma_cor_out(:,:) = sc(:,:,2)
     DEALLOCATE( en, sc, l_conv )
     !
     IF( notconv .NE. 0 ) THEN 
        CALL io_push_title("CONVERGENCE **NOT** ACHIEVED !!!")
     ENDIF
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
     !
     ! Output it per k-point
     !
     ALLOCATE(out_tab(qp_bandrange(2)-qp_bandrange(1)+1,7))
     !
     DO iks=1,nks
        DO ib = qp_bandrange(1), qp_bandrange(2)
           out_tab( ib - qp_bandrange(1) + 1, 1) = REAL( ib, KIND=DP) 
           out_tab( ib - qp_bandrange(1) + 1, 2) = et(ib,iks) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 3) = (et(ib,iks)+sigma_hf(ib,iks)) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 4) = qp_energy(ib,iks) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 5) = (qp_energy(ib,iks) - et(ib,iks) ) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 6) = REAL(  sigma_cor_out(ib,iks), KIND=DP ) * rytoev
           out_tab( ib - qp_bandrange(1) + 1, 7) = AIMAG( sigma_cor_out(ib,iks) ) * rytoev
        ENDDO
        WRITE(myglobk,'(i5.5)') iks_l2g(iks)
        !
        CALL serial_table_output(mpime==root,4000,'eqp_K'//myglobk,out_tab,&
           & qp_bandrange(2)-qp_bandrange(1)+1,7,&
           & (/'      band','    E0[eV]','   EHF[eV]','   Eqp[eV]','Eqp-E0[eV]','Sc_Eqp[eV]',' Width[eV]'/))
     ENDDO
     !
     DEALLOCATE( out_tab )  
Marco Govoni's avatar
Marco Govoni committed
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
     DEALLOCATE( sigma_cor_in )
     DEALLOCATE( sigma_cor_out )
     DEALLOCATE( z_in )
     DEALLOCATE( qp_energy )
     !
  ENDIF
  !
  IF( l_generate_plot ) THEN
     !
     CALL io_push_title("(P)lotting the QP corrections...")
     !
     CALL start_bar_type( barra, 'qplot', n_spectralf )
     !
1043
1044
     ALLOCATE( en(qp_bandrange(1):qp_bandrange(2),k_grid%nps,1) ) 
     ALLOCATE( sc(qp_bandrange(1):qp_bandrange(2),k_grid%nps,1) )
Govoni's avatar
Govoni committed
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
     !ALLOCATE( un(qp_bandrange(1):qp_bandrange(2),nks) )
     !
    !IF( mpime == 0 ) THEN 
    !   k = 0 
    !   DO iks=1,nks
    !      WRITE(iks_label,"(i5.5)") iks
    !      DO ib = qp_bandrange(1), qp_bandrange(2)
    !         WRITE(ib_label,"(i5.5)") ib
    !         fname = "o-eqp_K"//TRIM(iks_label)//"_B"//TRIM(ib_label)//".tab"
    !         k = k + 1
    !         un(ib,iks) = 10000+k
    !         OPEN( UNIT=un(ib,iks), FILE=TRIM(fname))
    !         WRITE( un(ib,iks),"(a)") '#        E[eV]         E-EKS  Re{S(E)}-Vxc      Im{S(E)}          A(E)'
    !      ENDDO  
    !      fname = "o-eqp_K"//TRIM(iks_label)//"_summed.tab"
    !      OPEN( UNIT=10000-iks, FILE=TRIM(fname))
    !      WRITE( 10000-iks,"(a)") '#        E[eV]          A(E)'
    !   ENDDO
    !ENDIF
Marco Govoni's avatar
Marco Govoni committed
1064
1065
1066
1067
1068
1069
1070
1071
1072
     !
     DO glob_ifreq = 1, n_spectralf
        sigma_freq(glob_ifreq) = (ecut_spectralf(2)-ecut_spectralf(1)) / REAL(n_spectralf-1,KIND=DP) * REAL(glob_ifreq-1,KIND=DP) &
                               & + ecut_spectralf(1)
     ENDDO 
     !
     DO glob_ifreq = 1, n_spectralf
        en = (ecut_spectralf(2)-ecut_spectralf(1)) / REAL(n_spectralf-1,KIND=DP) * REAL(glob_ifreq-1,KIND=DP) + ecut_spectralf(1) 
        CALL calc_corr_k( sc(:,:,1), en(:,:,1), .FALSE.)
1073
        DO iks=1, k_grid%nps
Marco Govoni's avatar
Marco Govoni committed
1074
1075
1076
1077
           DO ib = qp_bandrange(1), qp_bandrange(2)
              sigma_spectralf(glob_ifreq,ib,iks) = sc(ib,iks,1) 
           ENDDO
        ENDDO
Govoni's avatar
Govoni committed
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
       !IF( mpime == 0 ) THEN 
       !   DO iks=1,nks
       !      summed_sf=0._DP
       !      DO ib = qp_bandrange(1), qp_bandrange(2)
       !         WRITE( un(ib,iks),"(5f14.6)") en(ib,iks,1)*rytoev, &
       !         & (en(ib,iks,1)-et(ib,iks))*rytoev, &
       !         & (DBLE(sc(ib,iks,1))+sigma_hf(ib,iks))*rytoev, &
       !         & AIMAG(sc(ib,iks,1))*rytoev, &
       !         & ABS(AIMAG(sc(ib,iks,1)))/((en(ib,iks,1)-et(ib,iks)-sigma_hf(ib,iks)-DBLE(sc(ib,iks,1)))**2+&
       !         &AIMAG(sc(ib,iks,1))**2)/pi/rytoev
       !         IF( ib <= nbnd_occ(iks) ) THEN 
       !            summed_sf=summed_sf+ABS(AIMAG(sc(ib,iks,1)))/((en(ib,iks,1)-et(ib,iks)-sigma_hf(ib,iks)-DBLE(sc(ib,iks,1)))**2&
       !            &+AIMAG(sc(ib,iks,1))**2)/pi
       !         ENDIF
       !      ENDDO
       !      WRITE( 10000-iks,"(2f14.6)") en(qp_bandrange(1),iks,1)*rytoev,summed_sf/rytoev
       !   ENDDO
       !ENDIF
Marco Govoni's avatar
Marco Govoni committed
1096
1097
1098
        CALL update_bar_type(barra,'qplot',1)
     ENDDO
     !
Govoni's avatar
Govoni committed
1099
1100
1101
1102
1103
1104
1105
1106
    !IF( mpime == 0 ) THEN 
    !   DO iks=1,nks
    !      DO ib = qp_bandrange(1), qp_bandrange(2)
    !         CLOSE( un(ib,iks) )
    !      ENDDO 
    !      CLOSE(10000-iks)
    !   ENDDO
    !ENDIF 
Marco Govoni's avatar
Marco Govoni committed
1107
1108
1109
1110
1111
     !
     CALL stop_bar_type(barra,'qplot')
     !
     DEALLOCATE( en )
     DEALLOCATE( sc )  
Govoni's avatar
Govoni committed
1112
     !DEALLOCATE( un )
Marco Govoni's avatar
Marco Govoni committed
1113
     !
Govoni's avatar
Govoni committed
1114
     !CALL io_push_title('Done, take a look at the o-eqp_K*_B*.tab file(s) .')
Marco Govoni's avatar
Marco Govoni committed
1115
1116
1117
1118
1119
     ! 
  ENDIF 
  !
  DEALLOCATE( sigma_hf )
  !
Govoni's avatar
Govoni committed
1120
  CALL wfreq_db_write( )
Marco Govoni's avatar
Marco Govoni committed
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
  !
  CALL stop_clock( "solve_qp" )
  !
  !
END SUBROUTINE
!
!
SUBROUTINE output_eqp_report(iteration,en1,en2,sc1)
  !
  USE kinds,                ONLY : DP
Govoni's avatar
Govoni committed
1131
  USE westcom,              ONLY : qp_bandrange,trev_secant,logfile,iks_l2g
Marco Govoni's avatar
Marco Govoni committed
1132
1133
  USE pwcom,                ONLY : nks,et
  USE constants,            ONLY : rytoev
Govoni's avatar
Govoni committed
1134
  !USE west_io,              ONLY : serial_table_output
Marco Govoni's avatar
Marco Govoni committed
1135
  USE mp_world,             ONLY : mpime,root
1136
1137
  USE io_global,            ONLY : stdout
  USE io_push,              ONLY : io_push_title,io_push_bar
Govoni's avatar
Govoni committed
1138
  USE json_module,          ONLY : json_file
1139
  USE types_bz_grid,        ONLY : k_grid
Marco Govoni's avatar
Marco Govoni committed
1140
1141
1142
1143
1144
1145
  ! 
  IMPLICIT NONE
  !
  ! I/O
  !
  INTEGER,INTENT(IN) :: iteration
1146
1147
1148
  REAL(DP),INTENT(IN) :: en1(qp_bandrange(1):qp_bandrange(2),k_grid%nps)
  REAL(DP),INTENT(IN) :: en2(qp_bandrange(1):qp_bandrange(2),k_grid%nps)
  COMPLEX(DP),INTENT(IN) :: sc1(qp_bandrange(1):qp_bandrange(2),k_grid%nps)
Marco Govoni's avatar
Marco Govoni committed
1149
1150
1151
1152
1153
  !
  ! Workspace
  !
  CHARACTER(LEN=9) :: prefisso
  INTEGER :: contatore
1154
  REAL(DP) :: out_tabella(k_grid%nps*(qp_bandrange(2)-qp_bandrange(1)+1),7)
Marco Govoni's avatar
Marco Govoni committed
1155
  INTEGER :: ib, iks
1156
  CHARACTER(LEN=4) :: symb
Govoni's avatar
Govoni committed
1157
1158
1159
1160
  TYPE(json_file) :: json
  CHARACTER(LEN=6) :: my_label_k,my_label_b,citr
  INTEGER :: secitr, iunit
  LOGICAL :: found
1161
1162
1163
1164
1165
1166
1167
1168
  !
  ! STDOUT
  !
  WRITE(stdout,"(5X)")  
  CALL io_push_bar()
  IF( iteration >= 0 ) WRITE(stdout,"(5X,'Iter: ',i6.6)") iteration 
  WRITE(stdout,"(5X,a,1X,a,1X,a,1X,a)") 'K     ', 'B     ', 'QP energ. [eV]', 'conv'
  CALL io_push_bar()
1169
  DO iks = 1, k_grid%nps
1170
1171
1172
1173
1174
     DO ib = qp_bandrange(1), qp_bandrange(2)
        symb='  no'
        IF( (iteration .NE. 0) .AND. (ABS(en2(ib,iks)-en1(ib,iks)) < trev_secant) ) symb=' yes'
        WRITE(stdout,"(5X,i6.6,1X,i6.6,1X,1f14.6,1X,a)") iks, ib, en2(ib,iks) * rytoev, symb
     ENDDO
1175
     IF (k_grid%nps>1.AND.iks<k_grid%nps)  WRITE(stdout,"(5X, 33(a))") '---------------------------------'
1176
1177
1178
1179
  ENDDO
  CALL io_push_bar()
  !
  ! LOGFILE 
Marco Govoni's avatar
Marco Govoni committed
1180
  !
Govoni's avatar
Govoni committed
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
 !contatore=0
 !DO iks=1,nks
 !   DO ib = qp_bandrange(1), qp_bandrange(2)
 !      contatore = contatore + 1
 !      out_tabella(contatore,1) = REAL(iks,KIND=DP)
 !      out_tabella(contatore,2) = REAL(ib,KIND=DP)
 !      out_tabella(contatore,3) = et(ib,iks)*rytoev
 !      out_tabella(contatore,4) = en1(ib,iks)*rytoev
 !      out_tabella(contatore,5) = REAL( sc1(ib,iks), KIND=DP) *rytoev
 !      out_tabella(contatore,6) = en2(ib,iks)*rytoev
 !      out_tabella(contatore,7) = (en2(ib,iks)-en1(ib,iks))*rytoev
 !   ENDDO
 !ENDDO
 !!
 !IF(iteration>=0) THEN 
 !   WRITE(prefisso,"('itr_',i5.5)") iteration
 !ELSE
 !   prefisso='converged'
 !ENDIF
 !CALL serial_table_output(mpime==root,4000,'eqp.'//TRIM(ADJUSTL(prefisso)),out_tabella,&
 !& nks*(qp_bandrange(2)-qp_bandrange(1)+1),7,&
 !& (/'       iks','        ib','   Eks[eV]','   Ein[eV]','Sc_Ein[eV]','  Eout[eV]',' Diff.[eV]'/))
  !
  IF( mpime == root ) THEN 
     !
     CALL json%initialize()
     !
     CALL json%load_file(filename=TRIM(logfile))
     !
Govoni's avatar
Govoni committed
1210
     CALL json%get('exec.Q.secitr',secitr, found )
Govoni's avatar
Govoni committed
1211
1212
1213
1214
1215
1216
1217
1218
1219
     !
     IF( found ) THEN 
        secitr = secitr+1
     ELSE 
        secitr = 1
     ENDIF
     !
     WRITE(citr,'(i6)') secitr
     !
Govoni's avatar
Govoni committed
1220
     CALL json%update('exec.Q.secitr', secitr, found )
Govoni's avatar
Govoni committed
1221
     !
1222
     DO iks = 1, k_grid%nps
Govoni's avatar
Govoni committed
1223
1224
1225
        WRITE( my_label_k, '(i6.6)') iks_l2g(iks)
        DO ib = qp_bandrange(1), qp_bandrange(2) 
           WRITE( my_label_b, '(i6.6)') ib
Govoni's avatar
Govoni committed
1226
           CALL json%add('exec.Q.K'//TRIM(my_label_k)//'.B'//TRIM(my_label_b)//'.ein('//TRIM(ADJUSTL(citr))//')',&
Govoni's avatar
Govoni committed
1227
           & en1(ib,iks)*rytoev)
Govoni's avatar
Govoni committed
1228
           CALL json%add('exec.Q.K'//TRIM(my_label_k)//'.B'//TRIM(my_label_b)//'.eout('//TRIM(ADJUSTL(citr))//')',&
Govoni's avatar
Govoni committed
1229
           & en2(ib,iks)*rytoev)
Govoni's avatar
Govoni committed
1230
           CALL json%add('exec.Q.K'//TRIM(my_label_k)//'.B'//TRIM(my_label_b)//'.sc_ein('//TRIM(ADJUSTL(citr))//')',&
Govoni's avatar
Govoni committed
1231
1232
           & REAL(sc1(ib,iks))*rytoev)
        ENDDO
Marco Govoni's avatar
Marco Govoni committed
1233
     ENDDO
Govoni's avatar
Govoni committed
1234
1235
1236
1237
1238
1239
1240
1241
     !
     OPEN( NEWUNIT=iunit,FILE=TRIM(logfile) )
     CALL json%print_file( iunit )
     CLOSE( iunit )
     !
     CALL json%destroy()
     !
  ENDIF 
Marco Govoni's avatar
Marco Govoni committed
1242
1243
  !
END SUBROUTINE