solve_wfreq.f90 47.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
! 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_wfreq(l_read_restart,l_generate_plot)
  !-----------------------------------------------------------------------
  !
  USE control_flags,        ONLY : gamma_only 
  !
  IMPLICIT NONE
  !
  ! I/O 
  !
  LOGICAL,INTENT(IN) :: l_read_restart,l_generate_plot
  !
  IF( gamma_only ) THEN 
    CALL solve_wfreq_gamma( l_read_restart,l_generate_plot )
  ELSE
    CALL solve_wfreq_k( l_read_restart,l_generate_plot )
  ENDIF
  !
END SUBROUTINE 
!
!-----------------------------------------------------------------------
SUBROUTINE solve_wfreq_gamma(l_read_restart,l_generate_plot)
  !-----------------------------------------------------------------------
  !
  USE kinds,                ONLY : DP 
38
39
  USE westcom,              ONLY : west_prefix,n_pdep_eigen_to_use,n_lanczos,npwq,l_macropol,iks_l2g,d_epsm1_ifr,z_epsm1_rfr,&
                                 & l_enable_lanczos,nbnd_occ,iuwfc,lrwfc,wfreq_eta,imfreq_list,refreq_list,tr2_dfpt,&
40
                                 & z_head_rfr,d_head_ifr,o_restart_time,l_skip_nl_part_of_hcomr,npwqx,fftdriver, wstat_save_dir
Marco Govoni's avatar
Marco Govoni committed
41
  USE mp_global,            ONLY : my_image_id,nimage,inter_image_comm,intra_bgrp_comm
Marco Govoni's avatar
Marco Govoni committed
42
  USE mp_world,             ONLY : mpime,root
Marco Govoni's avatar
Marco Govoni committed
43
44
  USE mp,                   ONLY : mp_bcast,mp_barrier,mp_sum
  USE io_global,            ONLY : stdout,ionode
45
  USE gvect,                ONLY : g,ngm,gstart
Marco Govoni's avatar
Marco Govoni committed
46
47
48
49
50
51
52
53
54
55
56
  USE gvecw,                ONLY : gcutw
  USE cell_base,            ONLY : tpiba2,bg,omega
  USE fft_base,             ONLY : dffts
  USE constants,            ONLY : tpi,fpi,e2
  USE pwcom,                ONLY : npw,npwx,et,nks,current_spin,isk,xk,nbnd,lsda,igk_k,g2kin,current_k,wk,ngk
  USE wavefunctions_module, ONLY : evc,psic,psic_nc
  USE io_files,             ONLY : tmp_dir,nwordwfc,iunwfc
  USE fft_at_gamma,         ONLY : single_invfft_gamma,single_fwfft_gamma
!  USE fft_at_k,             ONLY : SINGLEBAND_INVFFT_k,SINGLEBAND_FWFFT_k
  USE becmod,               ONLY : becp,allocate_bec_type,deallocate_bec_type
  USE uspp,                 ONLY : vkb,nkb
57
  USE pdep_db,              ONLY : generate_pdep_fname
Marco Govoni's avatar
Marco Govoni committed
58
59
60
61
62
63
  USE pdep_io,              ONLY : pdep_read_G_and_distribute 
  USE io_push,              ONLY : io_push_title
!  USE control_flags,        ONLY : gamma_only
  USE noncollin_module,     ONLY : noncolin,npol 
  USE buffers,              ONLY : get_buffer
  USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
64
  USE distribution_center,  ONLY : pert,macropert,ifr,rfr,occband
Marco Govoni's avatar
Marco Govoni committed
65
66
  USE class_idistribute,    ONLY : idistribute 
  USE wfreq_restart,        ONLY : solvewfreq_restart_write,solvewfreq_restart_read,bks_type
67
  USE types_bz_grid,        ONLY : k_grid
68
  USE chi_invert,           ONLY : chi_invert_real, chi_invert_complex
69
  USE types_coulomb,        ONLY : pot3D
Marco Govoni's avatar
Marco Govoni committed
70
71
72
73
74
75
76
77
78
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  LOGICAL,INTENT(IN) :: l_read_restart,l_generate_plot
  !
  ! Workspace
  !
79
  INTEGER :: i1,i2,i3,im,ip,ig,glob_ip,ir,iv,ivloc,iks,ipol,m
80
81
  CHARACTER(LEN=25) :: filepot
  CHARACTER(LEN=:),ALLOCATABLE :: fname
Marco Govoni's avatar
Marco Govoni committed
82
83
84
85
86
87
  CHARACTER(LEN=6)      :: my_label_b
  COMPLEX(DP),ALLOCATABLE :: auxr(:)
  INTEGER :: nbndval
  REAL(DP),ALLOCATABLE :: diago( :, : ), subdiago( :, :), bnorm(:), braket(:, :, :)
  COMPLEX(DP),ALLOCATABLE :: q_s( :, :, : )
  COMPLEX(DP),ALLOCATABLE :: dvpsi(:,:)
88
  COMPLEX(DP),ALLOCATABLE :: phi(:,:), phis(:,:,:)
Marco Govoni's avatar
Marco Govoni committed
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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
  COMPLEX(DP),ALLOCATABLE :: phi_tmp(:,:)
  COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:)
  COMPLEX(DP) :: zkonstant
  TYPE(bar_type) :: barra
  INTEGER :: barra_load
  REAL(DP),ALLOCATABLE :: ps_r(:,:)
  TYPE(idistribute) :: mypara
  REAL(DP),ALLOCATABLE :: overlap(:,:)
  REAL(DP) :: mwo, ecv, dfactor, frequency, dhead
  COMPLEX(DP) :: zmwo, zfactor,zm,zp, zhead
  INTEGER :: glob_jp,ic,ifreq,il
  REAL(DP),ALLOCATABLE :: dmatilda(:,:), dlambda(:,:)
  COMPLEX(DP),ALLOCATABLE :: zmatilda(:,:), zlambda(:,:)
  REAL(DP),ALLOCATABLE :: dmati(:,:,:)
  COMPLEX(DP),ALLOCATABLE :: zmatr(:,:,:)
  LOGICAL :: l_iks_skip, l_iv_skip
  REAL(DP) :: time_spent(2)
  REAL(DP),EXTERNAL :: get_clock
  TYPE(bks_type) :: bks
  REAL(DP),ALLOCATABLE :: eprec(:)
  INTEGER :: ierr
  REAL(DP),ALLOCATABLE :: e(:)
  !
  CALL io_push_title("(W)-Lanczos")
  ! 
  ! DISTRIBUTING...
  !
  mypara = idistribute()
  IF(l_macropol) THEN
     CALL macropert%copyto(mypara)
  ELSE
     CALL pert%copyto(mypara)
  ENDIF
  !
  ! This is to reduce memory
  !
  CALL deallocate_bec_type( becp )
  IF(l_macropol) THEN 
     CALL allocate_bec_type ( nkb, MAX(mypara%nloc,3), becp ) ! I just need 2 becp at a time
  ELSE
     CALL allocate_bec_type ( nkb, mypara%nloc, becp ) ! I just need 2 becp at a time
  ENDIF
  !
  ! ALLOCATE dmati, zmatr, where chi0 is stored
  !
  ALLOCATE( dmati( mypara%nglob, mypara%nloc, ifr%nloc) )
  ALLOCATE( zmatr( mypara%nglob, mypara%nloc, rfr%nloc) )
  dmati = 0._DP
  zmatr = 0._DP
  !
  IF(l_read_restart) THEN
     CALL solvewfreq_restart_read( bks, dmati, zmatr, mypara%nglob, mypara%nloc )
  ELSE
     bks%lastdone_ks   = 0 
     bks%lastdone_band = 0 
     bks%old_ks        = 0 
     bks%old_band      = 0 
146
     bks%max_ks        = k_grid%nps 
Marco Govoni's avatar
Marco Govoni committed
147
148
149
150
     bks%min_ks        = 1 
  ENDIF
  !
  barra_load = 0
151
  DO iks = 1, k_grid%nps 
Marco Govoni's avatar
Marco Govoni committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
     IF(iks<bks%lastdone_ks) CYCLE
     DO iv = 1, nbnd_occ(iks)
        IF(iks==bks%lastdone_ks .AND. iv <= bks%lastdone_band ) CYCLE
        barra_load = barra_load + 1
     ENDDO
  ENDDO
  ! 
  IF( barra_load == 0 ) THEN 
     CALL start_bar_type ( barra, 'wlanczos', 1 )
     CALL update_bar_type( barra, 'wlanczos', 1 ) 
  ELSE
     CALL start_bar_type ( barra, 'wlanczos', barra_load )
  ENDIF
  !
166
  CALL pot3D%init('Wave',.FALSE.,'default')
167
  !
Marco Govoni's avatar
Marco Govoni committed
168
169
  ! LOOP 
  !
170
  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
Marco Govoni's avatar
Marco Govoni committed
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
     IF(iks<bks%lastdone_ks) CYCLE
     !
     ! ... 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
     !
186
     IF(k_grid%nps>1) THEN
Marco Govoni's avatar
Marco Govoni committed
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
        !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)
     !
     nbndval = nbnd_occ(iks)
     !
214
215
     mwo = - k_grid%weight(iks) / omega
     zmwo = CMPLX( - k_grid%weight(iks) / omega, 0._DP, KIND=DP)
Marco Govoni's avatar
Marco Govoni committed
216
217
218
219
     !
     bks%max_band = nbndval
     bks%min_band = 1
     !
220
221
     ! Parallel macropol
     IF(l_macropol) THEN
Marco Govoni's avatar
Marco Govoni committed
222
        !
223
        ! PHI 
Marco Govoni's avatar
Marco Govoni committed
224
        !
Marco Govoni's avatar
Marco Govoni committed
225
        CALL occband%init(nbndval,'i','occband',.FALSE.)
226
227
228
229
230
231
232
233
234
        !
        ALLOCATE(phis(npwx*npol,3,occband%nloc))
        !
        phis = 0._DP
        !
        ALLOCATE(phi(npwx*npol,3))
        ALLOCATE(phi_tmp(npwx*npol,3))
        !
        DO ivloc = 1, occband%nloc
Marco Govoni's avatar
Marco Govoni committed
235
           !
236
           iv = occband%l2g(ivloc)
Marco Govoni's avatar
Marco Govoni committed
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
           !
           CALL commutator_Hx_psi (iks, 1, 1, evc(1,iv), phi_tmp(1,1), l_skip_nl_part_of_hcomr)
           CALL commutator_Hx_psi (iks, 1, 2, evc(1,iv), phi_tmp(1,2), l_skip_nl_part_of_hcomr)
           CALL commutator_Hx_psi (iks, 1, 3, evc(1,iv), phi_tmp(1,3), l_skip_nl_part_of_hcomr)
           phi = 0._DP
           DO i1 = 1, 3
              DO i2 = 1, 3 
                 zkonstant = CMPLX( bg(i1,i2) , 0._DP, KIND=DP )  
                 CALL ZAXPY(npwx*npol,zkonstant,phi_tmp(1,i2),1,phi(1,i1),1) 
              ENDDO
           ENDDO
           !
           phi_tmp = phi
           CALL apply_alpha_pc_to_m_wfcs(nbndval,3,phi_tmp,(1._DP,0._DP))
           !
           ALLOCATE( eprec(3) )
           ALLOCATE( e(3) )
           !
           CALL set_eprec( 1, evc(1,iv), eprec(1))
           eprec(2) = eprec(1)
           eprec(3) = eprec(1)
           e(1) = et(iv,iks)
           e(2) = et(iv,iks)
           e(3) = et(iv,iks)
           !
           CALL precondition_m_wfcts( 3, phi_tmp, phi, eprec ) 
           !
           CALL linsolve_sternheimer_m_wfcts (nbndval, 3, phi_tmp, phi, e, eprec, tr2_dfpt, ierr ) 
           !
           IF(ierr/=0) THEN 
              WRITE(stdout, '(7X,"** WARNING : MACROPOL not converged, ierr = ",i8)') ierr
           ENDIF
           !
270
271
           phis(:,:,ivloc) = phi(:,:)
           !
Marco Govoni's avatar
Marco Govoni committed
272
273
           DEALLOCATE( eprec )
           DEALLOCATE( e )
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
           !
        END DO
        !
        DEALLOCATE( phi, phi_tmp )
        !
     ENDIF ! macropol
     !
     ALLOCATE(dvpsi(npwx*npol,mypara%nlocx))
     !
     time_spent(1) = get_clock( 'wlanczos' )
     !
     ! LOOP over band states
     !
     DO iv = 1, nbndval
        IF(iks==bks%lastdone_ks .AND. iv <= bks%lastdone_band ) CYCLE
        !
        ! MACROPOL CASE
        !
        IF(l_macropol) THEN
           !
           ALLOCATE(phi(npwx*npol,3))
           phi = 0._DP
           !
           DO ivloc = 1, occband%nloc
               IF( occband%l2g(ivloc) == iv ) THEN
                   phi(:,:) = phis(:,:,ivloc)
               ENDIF
           END DO
           !
           CALL mp_sum(phi, inter_image_comm)
Marco Govoni's avatar
Marco Govoni committed
304
305
306
307
308
           !
        ENDIF
        !
        ! PSIC
        !
309
        CALL single_invfft_gamma(dffts,npw,npwx,evc(1,iv),psic,'Wave') 
Marco Govoni's avatar
Marco Govoni committed
310
311
312
313
314
315
316
        !
        ! ZEROS
        !
        dvpsi = 0._DP   
        !
        ! Read PDEP
        !
317
        ALLOCATE( pertg(npwqx) )
Marco Govoni's avatar
Marco Govoni committed
318
319
320
321
322
323
324
325
326
327
328
329
        ALLOCATE( pertr( dffts%nnr ) )
        !
        DO ip=1,mypara%nloc
           !
           glob_ip = mypara%l2g(ip)
           !
           ! Decide whether read dbs E or dhpi 
           !
           IF(glob_ip<=n_pdep_eigen_to_use) THEN
              !
              ! Exhume dbs eigenvalue
              !
330
331
              CALL generate_pdep_fname( filepot, glob_ip ) 
              fname = TRIM( wstat_save_dir ) // "/"// filepot
Marco Govoni's avatar
Marco Govoni committed
332
333
334
              CALL pdep_read_G_and_distribute(fname,pertg)
              !
              ! Multiply by sqvc
335
              DO ig = 1, npwq
336
                 pertg(ig) = pot3D%sqvc(ig) * pertg(ig)
Marco Govoni's avatar
Marco Govoni committed
337
338
339
              ENDDO
              !
              ! Bring it to R-space
340
341
342
343
344
              CALL single_invfft_gamma(dffts,npwq,npwqx,pertg(1),pertr,TRIM(fftdriver))
              DO ir=1,dffts%nnr 
                 pertr(ir)=psic(ir)*pertr(ir)
              ENDDO
              CALL single_fwfft_gamma(dffts,npw,npwx,pertr,dvpsi(1,ip),'Wave')
Marco Govoni's avatar
Marco Govoni committed
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
              !
           ELSE
              !
              ipol = glob_ip-n_pdep_eigen_to_use
              !
              dvpsi(:,ip) = phi(:,ipol) * DSQRT(fpi * e2)
              !
              !
           ENDIF
           !
        ENDDO ! pert
        ! 
        DEALLOCATE(pertr)
        DEALLOCATE(pertg)
        IF(l_macropol) THEN
           DEALLOCATE(phi)
        ENDIF
        !
        CALL apply_alpha_pc_to_m_wfcs(nbndval,mypara%nloc,dvpsi,(1._DP,0._DP))
        !
        ! OVERLAP( glob_ip, im=1:nbnd ) = < psi_im iks | dvpsi_glob_ip >
        !
        IF(ALLOCATED(ps_r)) DEALLOCATE(ps_r)
        ALLOCATE(ps_r(nbnd-nbndval,mypara%nloc))
        CALL glbrak_gamma(evc(1,nbndval+1),dvpsi(1,1),ps_r,npw,npwx,nbnd-nbndval,mypara%nloc,nbnd-nbndval,npol)
        CALL mp_sum(ps_r,intra_bgrp_comm) 
        !
        IF(ALLOCATED(overlap)) DEALLOCATE(overlap)
        ALLOCATE(overlap(mypara%nglob, nbndval+1:nbnd ) )
        overlap = 0._DP
        DO ic = nbndval+1, nbnd
           DO ip = 1, mypara%nloc
              glob_ip = mypara%l2g(ip)
              overlap(glob_ip,ic) = ps_r(ic-nbndval,ip)
           ENDDO
        ENDDO
        !
        DEALLOCATE(ps_r)
        CALL mp_sum(overlap,inter_image_comm) 
        !
        ! Update dmati with cond
        !
        DO ifreq = 1, ifr%nloc
           !
           frequency = imfreq_list( ifreq ) 
           !
           DO ic = nbndval+1, nbnd 
              !
              ecv = et(ic,iks)-et(iv,iks)
              dfactor = mwo * 2._DP * ecv / ( ecv**2 + frequency**2 )
              !
              DO ip = 1, mypara%nloc
                 glob_ip = mypara%l2g(ip)
                 DO glob_jp = 1, mypara%nglob
                    !
                    dmati(glob_jp,ip,ifreq) = dmati(glob_jp,ip,ifreq) + overlap( glob_jp, ic) * overlap( glob_ip, ic) * dfactor
                    !
                 ENDDO
              ENDDO
              !
           ENDDO ! ic
        ENDDO ! ifreq
        !
        ! Update zmatr with cond
        !
        DO ifreq = 1, rfr%nloc
           !
           frequency = refreq_list( ifreq ) 
           !
           DO ic = nbndval+1, nbnd 
              !
              ecv = et(ic,iks)-et(iv,iks)
              zp = CMPLX( ecv + frequency, - wfreq_eta, KIND=DP )
              zm = CMPLX( ecv - frequency, - wfreq_eta, KIND=DP )
              zfactor = zmwo / zp + zmwo / zm 
              !
              DO ip = 1, mypara%nloc
                 glob_ip = mypara%l2g(ip)
                 DO glob_jp = 1, mypara%nglob
                    !
                    zmatr(glob_jp,ip,ifreq) = zmatr(glob_jp,ip,ifreq) + CMPLX( overlap( glob_jp, ic) * &
                    & overlap( glob_ip, ic), 0._DP, KIND=DP ) * zfactor
                    !
                 ENDDO
              ENDDO
              !
           ENDDO ! ic
        ENDDO ! ifreq 
        !
        DEALLOCATE(overlap)
        !
        ! Apply Pc, to be sure
        !
        CALL apply_alpha_pc_to_m_wfcs(nbnd,mypara%nloc,dvpsi,(1._DP,0._DP))
        !
        ! Now dvpsi is distributed according to eigen_distr (image), I need to use it for lanczos
        ! In the gamma_only case I need to process 2 dvpsi at a time (+ the odd last one, eventually), otherwise 1 at a time.
        !
        IF( l_enable_lanczos ) THEN  
           !
           ALLOCATE( bnorm    (                             mypara%nloc ) )
           ALLOCATE( diago    (               n_lanczos   , mypara%nloc ) )
           ALLOCATE( subdiago (               n_lanczos-1 , mypara%nloc ) )
           ALLOCATE( q_s      ( npwx*npol   , mypara%nloc , n_lanczos   ) )  ! WARNING ORDER INVERTED TO SMOOTHEN LANCZOS ALGORITHM 
           !
           CALL solve_deflated_lanczos_w_full_ortho ( nbnd, mypara%nloc, n_lanczos, dvpsi, diago, subdiago, q_s, bnorm)
           !
           ALLOCATE( braket   ( mypara%nglob, n_lanczos   , mypara%nloc ) )
Victor Yu's avatar
Victor Yu committed
453
           CALL get_brak_hyper_parallel(dvpsi,mypara%nloc,n_lanczos,q_s,braket,mypara)
Marco Govoni's avatar
Marco Govoni committed
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
           DEALLOCATE( q_s )
           !
           DO ip = 1, mypara%nloc
              CALL diago_lanczos( bnorm(ip), diago( :, ip), subdiago( :, ip), braket(:,:,ip), mypara%nglob )
           ENDDO
           !
           DEALLOCATE( bnorm )
           DEALLOCATE( subdiago )
           !
           ! Update dmati with lanczos
           !
           DO ifreq = 1, ifr%nloc
              !
              frequency = imfreq_list( ifreq ) 
              !
              DO il = 1, n_lanczos 
                 !
                 DO ip = 1, mypara%nloc
                    glob_ip = mypara%l2g(ip)
                    ecv = diago( il, ip ) - et(iv,iks) 
                    dfactor = mwo * 2._DP * ecv / ( ecv**2 + frequency**2 )
                    DO glob_jp = 1, mypara%nglob
                       !
                       dmati(glob_jp,ip,ifreq) = dmati(glob_jp,ip,ifreq) + braket( glob_jp, il, ip) * dfactor
                       !
                    ENDDO
                 ENDDO
                 !
              ENDDO ! il
           ENDDO ! ifreq
           !
           ! Update zmatr with lanczos
           !
           DO ifreq = 1, rfr%nloc
              !
              frequency = refreq_list( ifreq ) 
              !
              DO il = 1, n_lanczos 
                 !
                 DO ip = 1, mypara%nloc
                    glob_ip = mypara%l2g(ip)
                    ecv = diago( il, ip ) - et(iv,iks) 
                    zp = CMPLX( ecv + frequency, - wfreq_eta, KIND=DP )
                    zm = CMPLX( ecv - frequency, - wfreq_eta, KIND=DP )
                    zfactor = zmwo / zp + zmwo / zm
                    DO glob_jp = 1, mypara%nglob
                       !
                       zmatr(glob_jp,ip,ifreq) = zmatr(glob_jp,ip,ifreq) + &
                       & CMPLX( braket( glob_jp, il, ip), 0._DP, KIND=DP ) * zfactor
                       !
                    ENDDO
                 ENDDO
                 !
              ENDDO ! il
           ENDDO ! ifreq
           !
           ! MPI-IO 
           !
           !CALL writeout_solvewfreq( iks_l2g(iks), iv, diago, braket, io_comm, mypara%nloc, mypara%nglob, mypara%myoffset )
           !
           DEALLOCATE( diago ) 
           DEALLOCATE( braket )
           !
        ENDIF ! l_enable_lanczos
        !
        time_spent(2) = get_clock( 'wlanczos' ) 
        !
        IF( o_restart_time >= 0._DP ) THEN 
           IF( (time_spent(2)-time_spent(1)) > o_restart_time*60._DP .OR. iv == nbndval ) THEN 
              bks%lastdone_ks=iks
              bks%lastdone_band=iv
              CALL solvewfreq_restart_write(bks,dmati,zmatr,mypara%nglob,mypara%nloc)
              bks%old_ks = iks 
              bks%old_band = iv
              time_spent(1) = get_clock( 'wlanczos' )
           ENDIF
        ENDIF
        !
        CALL update_bar_type( barra, 'wlanczos', 1 )
        !
     ENDDO ! BANDS
     !
mahe's avatar
mahe committed
536
537
     IF(l_macropol) DEALLOCATE(phis)
     !
Marco Govoni's avatar
Marco Govoni committed
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
     DEALLOCATE(dvpsi)
     !
  ENDDO ! KPOINT-SPIN 
  !
  CALL stop_bar_type( barra, 'wlanczos' )
  !
  ! EPS-1 imfreq
  !
  ALLOCATE( dmatilda( mypara%nglob, mypara%nglob ) )
  ALLOCATE( dlambda( n_pdep_eigen_to_use, n_pdep_eigen_to_use) ) 
  ALLOCATE( d_epsm1_ifr( pert%nglob, pert%nloc, ifr%nloc) )
  d_epsm1_ifr = 0._DP
  IF(l_macropol) ALLOCATE( d_head_ifr( ifr%nloc) )
  !
  CALL start_clock('chi0_diag')
  !
  DO ifreq = 1, ifr%nloc
     !
     dmatilda = 0._DP
     DO ip = 1, mypara%nloc
        glob_ip = mypara%l2g(ip) 
        dmatilda( :, glob_ip) = dmati( :, ip, ifreq )
     ENDDO
     !
     CALL mp_sum( dmatilda, inter_image_comm )
     ! 
     CALL chi_invert_real( dmatilda, dhead, dlambda, mypara%nglob)
     !
     DO ip = 1, pert%nloc
        glob_ip = pert%l2g(ip)
        d_epsm1_ifr(1:n_pdep_eigen_to_use,ip,ifreq) = dlambda( 1:n_pdep_eigen_to_use, glob_ip)
     ENDDO 
     IF( l_macropol ) d_head_ifr( ifreq) = dhead
     !
  ENDDO
573
  !
Marco Govoni's avatar
Marco Govoni committed
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
  CALL stop_clock('chi0_diag')
  !
  DEALLOCATE( dlambda )
  DEALLOCATE( dmatilda )
  DEALLOCATE(dmati)
  !
  ! EPS-1 refreq
  !
  ALLOCATE( zmatilda( mypara%nglob, mypara%nglob ) )
  ALLOCATE( zlambda( n_pdep_eigen_to_use, n_pdep_eigen_to_use ) ) 
  ALLOCATE( z_epsm1_rfr( pert%nglob, pert%nloc, rfr%nloc) )
  z_epsm1_rfr = 0._DP
  IF(l_macropol) ALLOCATE( z_head_rfr( rfr%nloc) )
  !
  DO ifreq = 1, rfr%nloc
     !
     zmatilda = 0._DP
     DO ip = 1, mypara%nloc
        glob_ip = mypara%l2g(ip) 
        zmatilda( :, glob_ip) = zmatr( :, ip, ifreq )
     ENDDO
     !
     CALL mp_sum( zmatilda, inter_image_comm ) 
     CALL chi_invert_complex( zmatilda, zhead, zlambda, mypara%nglob)
     !
     DO ip = 1, pert%nloc
        glob_ip = pert%l2g(ip)
        z_epsm1_rfr(1:n_pdep_eigen_to_use,ip,ifreq) = zlambda( 1:n_pdep_eigen_to_use, glob_ip)
     ENDDO 
     IF( l_macropol ) z_head_rfr( ifreq) = zhead
     !
  ENDDO
  !
  DEALLOCATE( zlambda )
  DEALLOCATE( zmatilda )
  DEALLOCATE( zmatr ) 
  !
  IF(l_generate_plot) THEN
     CALL output_eps_head( )
  ENDIF
  !
END SUBROUTINE 
!
!-----------------------------------------------------------------------
SUBROUTINE solve_wfreq_k(l_read_restart,l_generate_plot)
  !-----------------------------------------------------------------------
  !
  USE kinds,                ONLY : DP 
622
  USE westcom,              ONLY : west_prefix,n_pdep_eigen_to_use,n_lanczos,npwq,l_macropol,iks_l2g,z_epsm1_ifr_q,&
623
624
                                 &  z_epsm1_rfr_q,l_enable_lanczos,nbnd_occ,iuwfc,lrwfc,wfreq_eta,imfreq_list,refreq_list,tr2_dfpt,&
                                 & z_head_rfr,z_head_ifr,o_restart_time,l_skip_nl_part_of_hcomr,npwqx,fftdriver, wstat_save_dir,&
625
                                 & ngq, igq_q
Marco Govoni's avatar
Marco Govoni committed
626
  USE mp_global,            ONLY : my_image_id,nimage,inter_image_comm,intra_bgrp_comm
Marco Govoni's avatar
Marco Govoni committed
627
  USE mp_world,             ONLY : mpime,root
Marco Govoni's avatar
Marco Govoni committed
628
629
  USE mp,                   ONLY : mp_bcast,mp_barrier,mp_sum
  USE io_global,            ONLY : stdout,ionode
630
  USE gvect,                ONLY : g,ngm,gstart
Marco Govoni's avatar
Marco Govoni committed
631
632
633
634
635
636
637
638
639
640
641
  USE gvecw,                ONLY : gcutw
  USE cell_base,            ONLY : tpiba2,bg,omega
  USE fft_base,             ONLY : dffts
  USE constants,            ONLY : tpi,fpi,e2
  USE pwcom,                ONLY : npw,npwx,et,nks,current_spin,isk,xk,nbnd,lsda,igk_k,g2kin,current_k,wk,ngk
  USE wavefunctions_module, ONLY : evc,psic,psic_nc
  USE io_files,             ONLY : tmp_dir,nwordwfc,iunwfc
!  USE fft_at_gamma,         ONLY : DOUBLEBAND_INVFFT,SINGLEBAND_INVFFT,DOUBLEBAND_FWFFT,SINGLEBAND_FWFFT
  USE fft_at_k,             ONLY : single_invfft_k,single_fwfft_k
  USE becmod,               ONLY : becp,allocate_bec_type,deallocate_bec_type
  USE uspp,                 ONLY : vkb,nkb
642
  USE pdep_db,              ONLY : generate_pdep_fname
Marco Govoni's avatar
Marco Govoni committed
643
644
645
646
647
648
  USE pdep_io,              ONLY : pdep_read_G_and_distribute 
  USE io_push,              ONLY : io_push_title
!  USE control_flags,        ONLY : gamma_only
  USE noncollin_module,     ONLY : noncolin,npol 
  USE buffers,              ONLY : get_buffer
  USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
mahe's avatar
mahe committed
649
  USE distribution_center,  ONLY : pert,macropert,ifr,rfr,occband
Marco Govoni's avatar
Marco Govoni committed
650
  USE class_idistribute,    ONLY : idistribute 
651
  USE wfreq_restart,        ONLY : solvewfreq_restart_write,solvewfreq_restart_read,bksq_type
652
  USE types_bz_grid,        ONLY : k_grid, q_grid, compute_phase
653
  USE chi_invert,           ONLY : chi_invert_complex
654
  USE types_coulomb,        ONLY : pot3D
Marco Govoni's avatar
Marco Govoni committed
655
656
657
658
659
660
661
662
663
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  LOGICAL,INTENT(IN) :: l_read_restart,l_generate_plot
  !
  ! Workspace
  !
mahe's avatar
mahe committed
664
  INTEGER :: i1,i2,i3,im,ip,ig,glob_ip,ir,iv,ivloc,iks,ik,is,iq,ikqs,ikq,ipol,m
665
666
  CHARACTER(LEN=25) :: filepot 
  CHARACTER(LEN=:),ALLOCATABLE    :: fname
Marco Govoni's avatar
Marco Govoni committed
667
  CHARACTER(LEN=6)      :: my_label_b
668
  CHARACTER(LEN=5)      :: my_label_q
Marco Govoni's avatar
Marco Govoni committed
669
670
671
672
673
674
  COMPLEX(DP),ALLOCATABLE :: auxr(:)
  INTEGER :: nbndval
  REAL(DP),ALLOCATABLE :: diago( :, : ), subdiago( :, :), bnorm(:)
  COMPLEX(DP),ALLOCATABLE :: braket(:, :, :)
  COMPLEX(DP),ALLOCATABLE :: q_s( :, :, : )
  COMPLEX(DP),ALLOCATABLE :: dvpsi(:,:)
mahe's avatar
mahe committed
675
  COMPLEX(DP),ALLOCATABLE :: phi(:,:), phis(:,:,:)
Marco Govoni's avatar
Marco Govoni committed
676
677
  COMPLEX(DP),ALLOCATABLE :: phi_tmp(:,:)
  COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:)
678
679
680
681
  COMPLEX(DP),ALLOCATABLE :: evckpq(:,:)
  COMPLEX(DP),ALLOCATABLE :: psick(:),psick_nc(:,:)
  COMPLEX(DP),ALLOCATABLE :: phase(:)
  INTEGER :: npwkq
Marco Govoni's avatar
Marco Govoni committed
682
683
684
685
686
687
688
689
690
691
  COMPLEX(DP) :: zkonstant
  TYPE(bar_type) :: barra
  INTEGER :: barra_load
  COMPLEX(DP),ALLOCATABLE :: ps_c(:,:)
  TYPE(idistribute) :: mypara
  COMPLEX(DP),ALLOCATABLE :: overlap(:,:)
  REAL(DP) :: mwo, ecv, dfactor, frequency
  COMPLEX(DP) :: zmwo, zfactor,zm,zp, zhead
  INTEGER :: glob_jp,ic,ifreq,il
  COMPLEX(DP),ALLOCATABLE :: zmatilda(:,:), zlambda(:,:)
692
693
  COMPLEX(DP),ALLOCATABLE :: zmati_q(:,:,:,:)
  COMPLEX(DP),ALLOCATABLE :: zmatr_q(:,:,:,:)
Marco Govoni's avatar
Marco Govoni committed
694
  LOGICAL :: l_iks_skip, l_iv_skip
695
  LOGICAL :: l_gammaq
Marco Govoni's avatar
Marco Govoni committed
696
697
  REAL(DP) :: time_spent(2)
  REAL(DP),EXTERNAL :: get_clock
698
  TYPE(bksq_type) :: bksq
Marco Govoni's avatar
Marco Govoni committed
699
700
701
  REAL(DP),ALLOCATABLE :: eprec(:)
  INTEGER :: ierr
  REAL(DP),ALLOCATABLE :: e(:)
702
  REAL(DP) :: kpq(3), g0(3) 
Marco Govoni's avatar
Marco Govoni committed
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
  !
  CALL io_push_title("(W)-Lanczos")
  ! 
  ! DISTRIBUTING...
  !
  mypara = idistribute()
  IF(l_macropol) THEN
     CALL macropert%copyto(mypara)
  ELSE
     CALL pert%copyto(mypara)
  ENDIF
  !
  ! This is to reduce memory
  !
  CALL deallocate_bec_type( becp ) 
  IF(l_macropol) THEN
     CALL allocate_bec_type ( nkb, MAX(mypara%nloc,3), becp ) ! I just need 2 becp at a time
  ELSE
     CALL allocate_bec_type ( nkb, mypara%nloc, becp ) ! I just need 2 becp at a time
  ENDIF
  !
724
  ! ALLOCATE zmati_q, zmatr_q, where chi0 is stored
Marco Govoni's avatar
Marco Govoni committed
725
  !
726
727
  ALLOCATE( zmati_q( mypara%nglob, mypara%nloc, ifr%nloc, q_grid%np) )
  ALLOCATE( zmatr_q( mypara%nglob, mypara%nloc, rfr%nloc, q_grid%np) )
728
729
730
731
732
733
734
735
736
737
  zmati_q = 0._DP
  zmatr_q = 0._DP
  !
  ALLOCATE( evckpq(npwx*npol,nbnd) )
  IF (noncolin) THEN
     ALLOCATE( psick_nc(dffts%nnr,npol) )
  ELSE
     ALLOCATE( psick(dffts%nnr) )
  ENDIF
  ALLOCATE( phase(dffts%nnr) )
Marco Govoni's avatar
Marco Govoni committed
738
739
  !
  IF(l_read_restart) THEN
740
     CALL solvewfreq_restart_read( bksq, zmati_q, zmatr_q, mypara%nglob, mypara%nloc )
Marco Govoni's avatar
Marco Govoni committed
741
  ELSE
742
743
744
745
746
747
748
749
750
751
     bksq%lastdone_ks   = 0 
     bksq%lastdone_q    = 0 
     bksq%lastdone_band = 0 
     bksq%old_ks        = 0 
     bksq%old_q         = 0 
     bksq%old_band      = 0 
     bksq%max_q         = q_grid%np
     bksq%max_ks        = k_grid%np
     bksq%min_q         = 1 
     bksq%min_ks        = 1 
Marco Govoni's avatar
Marco Govoni committed
752
753
754
  ENDIF
  !
  barra_load = 0
755
  DO iq = 1, q_grid%np
756
     IF (iq<bksq%lastdone_q) CYCLE
757
     DO iks = 1, k_grid%nps
758
        IF(iq==bksq%lastdone_q .AND. iks<bksq%lastdone_ks) CYCLE
759
        DO iv = 1, nbnd_occ(iks)
760
           IF(iq==bksq%lastdone_q .AND. iks==bksq%lastdone_ks .AND. iv <= bksq%lastdone_band ) CYCLE
761
762
           barra_load = barra_load + 1
        ENDDO
Marco Govoni's avatar
Marco Govoni committed
763
764
765
766
767
768
769
770
771
772
773
774
     ENDDO
  ENDDO
  ! 
  IF( barra_load == 0 ) THEN 
     CALL start_bar_type ( barra, 'wlanczos', 1 )
     CALL update_bar_type( barra, 'wlanczos', 1 ) 
  ELSE
     CALL start_bar_type ( barra, 'wlanczos', barra_load )
  ENDIF
  !
  ! LOOP 
  !
775
  DO iq = 1, q_grid%np   ! Q-POINT
776
     IF (iq<bksq%lastdone_q) CYCLE
Marco Govoni's avatar
Marco Govoni committed
777
     !
778
     npwq = ngq(iq)
779
     l_gammaq = q_grid%l_pIsGamma(iq)
Marco Govoni's avatar
Marco Govoni committed
780
     !
781
     CALL pot3D%init('Wave',.TRUE.,'default',iq)
Marco Govoni's avatar
Marco Govoni committed
782
     !
783
     DO iks = 1, k_grid%nps   ! KPOINT-SPIN
784
        !
785
        ik = k_grid%ip(iks) 
786
787
        is = k_grid%is(iks) 
        !
788
        IF(iq==bksq%lastdone_q .AND. iks<bksq%lastdone_ks) CYCLE
Marco Govoni's avatar
Marco Govoni committed
789
        !
790
        ! ... Set k-point, spin, kinetic energy, needed by Hpsi
Marco Govoni's avatar
Marco Govoni committed
791
        !
792
793
794
        current_k = iks
        IF ( lsda ) current_spin = isk(iks)
        call g2_kin( iks )
Marco Govoni's avatar
Marco Govoni committed
795
        !
796
        ! ... More stuff needed by the hamiltonian: nonlocal projectors
Marco Govoni's avatar
Marco Govoni committed
797
        !
798
        IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), k_grid%p_cart(1,ik), vkb )
799
        npw = ngk(iks)
Marco Govoni's avatar
Marco Govoni committed
800
        !
801
        ! ... read in wavefunctions from the previous iteration
Marco Govoni's avatar
Marco Govoni committed
802
        !
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
        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)
Marco Govoni's avatar
Marco Govoni committed
828
        !
829
830
        !ikqs = kpq_grid%index_kq(iks,iq)
        !npwkq = ngk(ikqs)
Marco Govoni's avatar
Marco Govoni committed
831
        !
832
833
834
        !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)                                                                                   
Marco Govoni's avatar
Marco Govoni committed
835
        !
836
837
        npwkq = ngk(ikqs)
        !
838
        CALL compute_phase( g0, 'cart', phase )
839
840
841
        !
        ! Set wavefunctions at [k+q] in G space, for all bands,
        ! and store them in evckpq
Marco Govoni's avatar
Marco Govoni committed
842
        ! 
843
844
        IF ( my_image_id == 0 ) CALL get_buffer( evckpq, lrwfc, iuwfc, ikqs )
        CALL mp_bcast( evckpq, 0, inter_image_comm )
Marco Govoni's avatar
Marco Govoni committed
845
        !
846
        nbndval = nbnd_occ(iks)
Marco Govoni's avatar
Marco Govoni committed
847
        !
848
849
        mwo = - k_grid%weight(iks) / omega
        zmwo = CMPLX( - k_grid%weight(iks) / omega, 0._DP, KIND=DP)
Marco Govoni's avatar
Marco Govoni committed
850
        !
851
852
        bksq%max_band = nbndval
        bksq%min_band = 1
Marco Govoni's avatar
Marco Govoni committed
853
        !
mahe's avatar
mahe committed
854
        ! MACROPOL CASE
Marco Govoni's avatar
Marco Govoni committed
855
        !
mahe's avatar
mahe committed
856
        IF(l_macropol .AND. l_gammaq) THEN
Marco Govoni's avatar
Marco Govoni committed
857
           !
mahe's avatar
mahe committed
858
           ! PHI 
Marco Govoni's avatar
Marco Govoni committed
859
           !
mahe's avatar
mahe committed
860
861
862
863
864
865
866
867
868
869
           CALL occband%init(nbndval,'i','occband',.TRUE.)
           !
           ALLOCATE(phis(npwx*npol,3,occband%nloc))
           !
           phis = 0._DP
           !
           ALLOCATE(phi(npwx*npol,3))
           ALLOCATE(phi_tmp(npwx*npol,3))
           !
           DO ivloc = 1, occband%nloc
Marco Govoni's avatar
Marco Govoni committed
870
              !
mahe's avatar
mahe committed
871
              iv = occband%l2g(ivloc)
Marco Govoni's avatar
Marco Govoni committed
872
              !
873
874
875
876
877
878
879
880
              CALL commutator_Hx_psi (iks, 1, 1, evc(1,iv), phi_tmp(1,1), l_skip_nl_part_of_hcomr)
              CALL commutator_Hx_psi (iks, 1, 2, evc(1,iv), phi_tmp(1,2), l_skip_nl_part_of_hcomr)
              CALL commutator_Hx_psi (iks, 1, 3, evc(1,iv), phi_tmp(1,3), l_skip_nl_part_of_hcomr)
              phi = 0._DP
              DO i1 = 1, 3
                 DO i2 = 1, 3 
                    zkonstant = CMPLX( bg(i1,i2) , 0._DP, KIND=DP )  
                    CALL ZAXPY(npwx*npol,zkonstant,phi_tmp(1,i2),1,phi(1,i1),1) 
Marco Govoni's avatar
Marco Govoni committed
881
882
883
                 ENDDO
              ENDDO
              !
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
              phi_tmp = phi
              CALL apply_alpha_pc_to_m_wfcs(nbndval,3,phi_tmp,(1._DP,0._DP))
              !
              ALLOCATE( eprec(3) )
              ALLOCATE( e(3) )
              !
              CALL set_eprec( 1, evc(1,iv), eprec(1))
              eprec(2) = eprec(1)
              eprec(3) = eprec(1)
              e(1) = et(iv,iks)
              e(2) = et(iv,iks)
              e(3) = et(iv,iks)
              !
              CALL precondition_m_wfcts( 3, phi_tmp, phi, eprec ) 
              !
              CALL linsolve_sternheimer_m_wfcts (nbndval, 3, phi_tmp, phi, e, eprec, tr2_dfpt, ierr ) 
              !
              IF(ierr/=0) THEN 
                 WRITE(stdout, '(7X,"** WARNING : MACROPOL not converged, ierr = ",i8)') ierr
              ENDIF
              !
mahe's avatar
mahe committed
905
906
              phis(:,:,ivloc) = phi(:,:)
              !
907
908
909
              DEALLOCATE( eprec )
              DEALLOCATE( e )
              !
mahe's avatar
mahe committed
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
           END DO
           !
           DEALLOCATE( phi, phi_tmp )
           !
        ENDIF ! macropol
        !
        ALLOCATE(dvpsi(npwx*npol,mypara%nlocx))
        !
        time_spent(1) = get_clock( 'wlanczos' )
        !
        ! LOOP over band states
        !
        DO iv = 1, nbndval
           IF(iq==bksq%lastdone_q .AND. iks==bksq%lastdone_ks .AND. iv <= bksq%lastdone_band ) CYCLE
           !
           ! MACROPOL CASE
           !
           IF(l_macropol .AND. l_gammaq) THEN
              !
              ALLOCATE(phi(npwx*npol,3))
              phi = 0._DP
              !
              DO ivloc = 1, occband%nloc
                 IF( occband%l2g(ivloc) == iv ) THEN
                    phi(:,:) = phis(:,:,ivloc)
                 ENDIF
              END DO
              !
              CALL mp_sum(phi, inter_image_comm)
939
940
              !
           ENDIF
Marco Govoni's avatar
Marco Govoni committed
941
           !
942
           ! PSIC
Marco Govoni's avatar
Marco Govoni committed
943
           !
944
945
           IF(noncolin) THEN
              CALL single_invfft_k(dffts,npwkq,npwx,evckpq(1     ,iv),psick_nc(1,1),'Wave',igk_k(1,ikqs))
946
              CALL single_invfft_k(dffts,npwkq,npwx,evckpq(npwx+1,iv),psick_nc(1,2),'Wave',igk_k(1,ikqs))
947
948
949
950
951
952
953
954
955
956
957
958
959
960
           ELSE
              CALL single_invfft_k(dffts,npwkq,npwx,evckpq(1,iv),psick,'Wave',igk_k(1,ikqs))
           ENDIF
           !
           ! ZEROS
           !
           dvpsi = 0._DP   
           !
           ! Read PDEP
           !
           ALLOCATE( pertg(npwqx) )
           ALLOCATE( pertr( dffts%nnr ) )
           !
           DO ip=1,mypara%nloc
Marco Govoni's avatar
Marco Govoni committed
961
              !
962
              glob_ip = mypara%l2g(ip)
Marco Govoni's avatar
Marco Govoni committed
963
              !
964
965
966
967
968
969
              ! Decide whether read dbs E or dhpi 
              !
              IF(glob_ip<=n_pdep_eigen_to_use) THEN
                 !
                 ! Exhume dbs eigenvalue
                 !
970
971
                 CALL generate_pdep_fname( filepot, glob_ip, iq ) 
                 fname = TRIM( wstat_save_dir ) // "/"// filepot
972
973
974
975
                 CALL pdep_read_G_and_distribute(fname,pertg,iq)
                 !
                 ! Multiply by sqvc
                 DO ig = 1, npwq
976
                    pertg(ig) = pot3D%sqvc(ig) * pertg(ig)
Marco Govoni's avatar
Marco Govoni committed
977
                 ENDDO
978
979
980
981
982
983
984
985
                 !
                 ! Bring it to R-space
                 IF(noncolin) THEN
                    CALL single_invfft_k(dffts,npwq,npwqx,pertg(1),pertr,'Wave',igq_q(1,iq))
                    DO ir=1,dffts%nnr 
                       pertr(ir)=phase(ir)*psick_nc(ir,1)*DCONJG(pertr(ir))
                    ENDDO
                    CALL single_fwfft_k(dffts,npw,npwx,pertr,dvpsi(1,ip),'Wave',igk_k(1,current_k))
986
                    CALL single_invfft_k(dffts,npwq,npwqx,pertg(1),pertr,'Wave',igq_q(1,iq))
987
988
989
                    DO ir=1,dffts%nnr 
                       pertr(ir)=phase(ir)*psick_nc(ir,2)*DCONJG(pertr(ir))
                    ENDDO
990
                    CALL single_fwfft_k(dffts,npw,npwx,pertr,dvpsi(npwx+1,ip),'Wave',igk_k(1,current_k))
991
992
993
994
995
996
997
998
999
1000
                 ELSE
                    CALL single_invfft_k(dffts,npwq,npwqx,pertg(1),pertr,'Wave',igq_q(1,iq))
                    DO ir=1,dffts%nnr 
                       pertr(ir)=phase(ir)*psick(ir)*DCONJG(pertr(ir))
                    ENDDO
                    CALL single_fwfft_k(dffts,npw,npwx,pertr,dvpsi(1,ip),'Wave',igk_k(1,current_k))
                 ENDIF 
                 !
              ELSE
                 !
1001
1002
1003
1004
                 IF (l_gammaq) THEN 
                    ipol = glob_ip-n_pdep_eigen_to_use
                    dvpsi(:,ip) = phi(:,ipol) * DSQRT(fpi * e2)
                 ENDIF
1005
1006
                 !
              ENDIF
Marco Govoni's avatar
Marco Govoni committed
1007
              !
1008
1009
1010
1011
1012
1013
1014
           ENDDO ! pert
           ! 
           DEALLOCATE(pertr)
           DEALLOCATE(pertg)
           IF(l_macropol.AND.l_gammaq) THEN
              DEALLOCATE(phi)
           ENDIF
Marco Govoni's avatar
Marco Govoni committed
1015
           !
1016
           CALL apply_alpha_pc_to_m_wfcs(nbndval,mypara%nloc,dvpsi,(1._DP,0._DP))
Marco Govoni's avatar
Marco Govoni committed
1017
           !
1018
           ! OVERLAP( glob_ip, im=1:nbnd ) = < psi_im iks | dvpsi_glob_ip >
Marco Govoni's avatar
Marco Govoni committed
1019
           !
1020
1021
1022
1023
           IF(ALLOCATED(ps_c)) DEALLOCATE(ps_c)
           ALLOCATE(ps_c(nbnd-nbndval,mypara%nloc))
           CALL glbrak_k(evc(1,nbndval+1),dvpsi(1,1),ps_c,npw,npwx,nbnd-nbndval,mypara%nloc,nbnd-nbndval,npol)
           CALL mp_sum(ps_c,intra_bgrp_comm) 
Marco Govoni's avatar
Marco Govoni committed
1024
           !
1025
1026
1027
1028
1029
1030
1031
1032
           IF(ALLOCATED(overlap)) DEALLOCATE(overlap)
           ALLOCATE(overlap(mypara%nglob, nbndval+1:nbnd ) )
           overlap = 0._DP
           DO ic = nbndval+1, nbnd
              DO ip = 1, mypara%nloc
                 glob_ip = mypara%l2g(ip)
                 overlap(glob_ip,ic) = ps_c(ic-nbndval,ip)
              ENDDO
Marco Govoni's avatar
Marco Govoni committed
1033
1034
           ENDDO
           !
1035
1036
           DEALLOCATE(ps_c)
           CALL mp_sum(overlap,inter_image_comm) 
Marco Govoni's avatar
Marco Govoni committed
1037
           !
1038
           ! Update zmati with cond
Marco Govoni's avatar
Marco Govoni committed
1039
1040
1041
1042
1043
           !
           DO ifreq = 1, ifr%nloc
              !
              frequency = imfreq_list( ifreq ) 
              !
1044
1045
1046
1047
1048
              DO ic = nbndval+1, nbnd 
                 !
                 ecv = et(ic,iks)-et(iv,ikqs)
                 dfactor = mwo * 2._DP * ecv / ( ecv**2 + frequency**2 )
                 zfactor = CMPLX( dfactor, 0._DP, KIND=DP )
Marco Govoni's avatar
Marco Govoni committed
1049
1050
1051
1052
1053
                 !
                 DO ip = 1, mypara%nloc
                    glob_ip = mypara%l2g(ip)
                    DO glob_jp = 1, mypara%nglob
                       !
1054
1055
                       zmati_q(glob_jp,ip,ifreq,iq) = zmati_q(glob_jp,ip,ifreq,iq) + DCONJG(overlap( glob_ip, ic)) * &
                       & overlap( glob_jp, ic) * zfactor
Marco Govoni's avatar
Marco Govoni committed
1056
1057
1058
1059
                       !
                    ENDDO
                 ENDDO
                 !
1060
              ENDDO ! ic
Marco Govoni's avatar
Marco Govoni committed
1061
1062
           ENDDO ! ifreq
           !
1063
           ! Update zmatr with cond
Marco Govoni's avatar
Marco Govoni committed
1064
1065
1066
1067
1068
           !
           DO ifreq = 1, rfr%nloc
              !
              frequency = refreq_list( ifreq ) 
              !
1069
1070
1071
1072
1073
1074
              DO ic = nbndval+1, nbnd 
                 !
                 ecv = et(ic,iks)-et(iv,ikqs)
                 zp = CMPLX( ecv + frequency, - wfreq_eta, KIND=DP )
                 zm = CMPLX( ecv - frequency, - wfreq_eta, KIND=DP )
                 zfactor = zmwo / zp + zmwo / zm 
Marco Govoni's avatar
Marco Govoni committed
1075
1076
1077
1078
1079
                 !
                 DO ip = 1, mypara%nloc
                    glob_ip = mypara%l2g(ip)
                    DO glob_jp = 1, mypara%nglob
                       !
1080
1081
                       zmatr_q(glob_jp,ip,ifreq,iq) = zmatr_q(glob_jp,ip,ifreq,iq) + DCONJG( overlap( glob_ip, ic) ) * &
                       & overlap( glob_jp, ic) * zfactor
Marco Govoni's avatar
Marco Govoni committed
1082
1083
1084
1085
                       !
                    ENDDO
                 ENDDO
                 !
1086
1087
              ENDDO ! ic
           ENDDO ! ifreq 
Marco Govoni's avatar
Marco Govoni committed
1088
           !
1089
           DEALLOCATE(overlap)
Marco Govoni's avatar
Marco Govoni committed
1090
           !
1091
           ! Apply Pc, to be sure
Marco Govoni's avatar
Marco Govoni committed
1092
           !
1093
           CALL apply_alpha_pc_to_m_wfcs(nbnd,mypara%nloc,dvpsi,(1._DP,0._DP))
Marco Govoni's avatar
Marco Govoni committed
1094
           !
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
           ! Now dvpsi is distributed according to eigen_distr (image), I need to use it for lanczos
           ! In the gamma_only case I need to process 2 dvpsi at a time (+ the odd last one, eventually), otherwise 1 at a time.
           !
           IF( l_enable_lanczos ) THEN  
              !
              ALLOCATE( bnorm    (                             mypara%nloc ) )
              ALLOCATE( diago    (               n_lanczos   , mypara%nloc ) )
              ALLOCATE( subdiago (               n_lanczos-1 , mypara%nloc ) )
              ALLOCATE( q_s      ( npwx*npol   , mypara%nloc , n_lanczos   ) )  ! WARNING ORDER INVERTED TO SMOOTHEN LANCZOS ALGORITHM 
              !
              CALL solve_deflated_lanczos_w_full_ortho ( nbnd, mypara%nloc, n_lanczos, dvpsi, diago, subdiago, q_s, bnorm)
              !
              ALLOCATE( braket   ( mypara%nglob, n_lanczos   , mypara%nloc ) )
Victor Yu's avatar
Victor Yu committed
1108
              CALL get_brak_hyper_parallel_complex(dvpsi,mypara%nloc,n_lanczos,q_s,braket,mypara)
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
              DEALLOCATE( q_s )
              !
              DO ip = 1, mypara%nloc
                 CALL diago_lanczos_complex( bnorm(ip), diago( :, ip), subdiago( :, ip), braket(:,:,ip), mypara%nglob )
              ENDDO
              !
              DEALLOCATE( bnorm )
              DEALLOCATE( subdiago )
              !
              ! Update zmati with lanczos
              !
              DO ifreq = 1, ifr%nloc
                 !
                 frequency = imfreq_list( ifreq ) 
                 !
                 DO il = 1, n_lanczos 
                    !
                    DO ip = 1, mypara%nloc
                       glob_ip = mypara%l2g(ip)
                       ecv = diago( il, ip ) - et(iv,ikqs) 
                       dfactor = mwo * 2._DP * ecv / ( ecv**2 + frequency**2 )
                       zfactor = CMPLX( dfactor, 0._DP, KIND=DP )
                       DO glob_jp = 1, mypara%nglob
                          !
                          zmati_q(glob_jp,ip,ifreq,iq) = zmati_q(glob_jp,ip,ifreq,iq) + DCONJG(braket( glob_jp, il, ip)) * zfactor
                          !
                       ENDDO
                    ENDDO
                    !
                 ENDDO ! il
              ENDDO ! ifreq
              !
              ! Update zmatr with lanczos
              !
              DO ifreq = 1, rfr%nloc
                 !
                 frequency = refreq_list( ifreq ) 
                 !
                 DO il = 1, n_lanczos 
                    !
                    DO ip = 1, mypara%nloc
                       glob_ip = mypara%l2g(ip)
                       ecv = diago( il, ip ) - et(iv,ikqs) 
                       zp = CMPLX( ecv + frequency, - wfreq_eta, KIND=DP )
                       zm = CMPLX( ecv - frequency, - wfreq_eta, KIND=DP )
                       zfactor = zmwo / zp + zmwo / zm
                       DO glob_jp = 1, mypara%nglob
                          !
                          zmatr_q(glob_jp,ip,ifreq,iq) = zmatr_q(glob_jp,ip,ifreq,iq) + DCONJG(braket( glob_jp, il, ip)) * zfactor
                          !
                       ENDDO
                    ENDDO
                    !
                 ENDDO ! il
              ENDDO ! ifreq
              !
              ! MPI-IO 
              !
              !CALL writeout_solvewfreq( iks_l2g(iks), iv, diago, braket, io_comm, mypara%nloc, mypara%nglob, mypara%myoffset )
              !
              DEALLOCATE( diago ) 
              DEALLOCATE( braket )
              !
           ENDIF ! l_enable_lanczos
           !
           time_spent(2) = get_clock( 'wlanczos' ) 
           !
           IF( o_restart_time >= 0._DP ) THEN 
              IF( (time_spent(2)-time_spent(1)) > o_restart_time*60._DP .OR. iv == nbndval) THEN
1178
1179
1180
1181
1182
1183
1184
                 bksq%lastdone_q=iq
                 bksq%lastdone_ks=iks
                 bksq%lastdone_band=iv
                 CALL solvewfreq_restart_write(bksq,zmati_q(:,:,:,:),zmatr_q(:,:,:,:),mypara%nglob,mypara%nloc)
                 bksq%old_q = iq
                 bksq%old_ks = iks 
                 bksq%old_band = iv
1185
1186
                 time_spent(1) = get_clock( 'wlanczos' )
              ENDIF
Marco Govoni's avatar
Marco Govoni committed
1187
           ENDIF
1188
1189
1190
1191
           !
           CALL update_bar_type( barra, 'wlanczos', 1 )
           !
        ENDDO ! BANDS
Marco Govoni's avatar
Marco Govoni committed
1192
        !
mahe's avatar
mahe committed
1193
1194
        IF(l_macropol .AND. l_gammaq) DEALLOCATE(phis)
        !
1195
        DEALLOCATE(dvpsi)
Marco Govoni's avatar
Marco Govoni committed
1196
        !
1197
     ENDDO ! KPOINT-SPIN 
Marco Govoni's avatar
Marco Govoni committed
1198
     !
1199
  ENDDO ! QPOINT
Marco Govoni's avatar
Marco Govoni committed
1200
1201
1202
  !
  CALL stop_bar_type( barra, 'wlanczos' )
  !
1203
  DEALLOCATE( evckpq )
1204
1205
1206
1207
1208
  IF (noncolin) THEN
     DEALLOCATE( psick_nc )
  ELSE
     DEALLOCATE( psick )
  ENDIF
1209
1210
  DEALLOCATE( phase )
  !
Marco Govoni's avatar
Marco Govoni committed
1211
1212
1213
1214
  ! EPS-1 imfreq
  !
  ALLOCATE( zmatilda( mypara%nglob, mypara%nglob ) )
  ALLOCATE( zlambda( n_pdep_eigen_to_use, n_pdep_eigen_to_use) ) 
1215
  ALLOCATE( z_epsm1_ifr_q( pert%nglob, pert%nloc, ifr%nloc, q_grid%np) )
1216
  z_epsm1_ifr_q = 0._DP
Marco Govoni's avatar
Marco Govoni committed
1217
1218
1219
1220
  IF(l_macropol) ALLOCATE( z_head_ifr( ifr%nloc) )
  !
  CALL start_clock('chi0_diag')
  !
1221
  DO iq = 1, q_grid%np
Marco Govoni's avatar
Marco Govoni committed
1222
     !
1223
     l_gammaq = q_grid%l_pIsGamma(iq)
Marco Govoni's avatar
Marco Govoni committed
1224
     !
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
     DO ifreq = 1, ifr%nloc
        !
        zmatilda = 0._DP
        DO ip = 1, mypara%nloc
           glob_ip = mypara%l2g(ip) 
           zmatilda( :, glob_ip) = zmati_q( :, ip, ifreq, iq )
        ENDDO
        !
        CALL mp_sum( zmatilda, inter_image_comm )
        ! 
1235
        CALL chi_invert_complex( zmatilda, zhead, zlambda, mypara%nglob, l_gammaq)
1236
1237
1238
1239
1240
1241
1242
1243
        !
        DO ip = 1, pert%nloc
           glob_ip = pert%l2g(ip)
           z_epsm1_ifr_q(1:n_pdep_eigen_to_use,ip,ifreq,iq) = zlambda( 1:n_pdep_eigen_to_use, glob_ip)
        ENDDO 
        IF( l_macropol .AND. l_gammaq ) z_head_ifr( ifreq) = zhead
        !
     ENDDO
Marco Govoni's avatar
Marco Govoni committed
1244
1245
     !
  ENDDO
1246
  !
Marco Govoni's avatar
Marco Govoni committed
1247
1248
1249
1250
  CALL stop_clock('chi0_diag')
  !
  DEALLOCATE( zlambda )
  DEALLOCATE( zmatilda )
1251
  DEALLOCATE( zmati_q )
Marco Govoni's avatar
Marco Govoni committed
1252
1253
1254
1255
1256
  !
  ! EPS-1 refreq
  !
  ALLOCATE( zmatilda( mypara%nglob, mypara%nglob ) )
  ALLOCATE( zlambda( n_pdep_eigen_to_use, n_pdep_eigen_to_use ) ) 
1257
  ALLOCATE( z_epsm1_rfr_q( pert%nglob, pert%nloc, rfr%nloc, q_grid%np) )
1258
  z_epsm1_rfr_q = 0._DP
Marco Govoni's avatar
Marco Govoni committed
1259
1260
  IF(l_macropol) ALLOCATE( z_head_rfr( rfr%nloc) )
  !
1261
  DO iq = 1, q_grid%np
Marco Govoni's avatar
Marco Govoni committed
1262
     !
1263
     l_gammaq = q_grid%l_pIsGamma(iq)
Marco Govoni's avatar
Marco Govoni committed
1264
     !
1265
1266
1267
1268
1269
1270
1271
1272
1273
     DO ifreq = 1, rfr%nloc
        !
        zmatilda = 0._DP
        DO ip = 1, mypara%nloc
           glob_ip = mypara%l2g(ip) 
           zmatilda( :, glob_ip) = zmatr_q( :, ip, ifreq, iq )
        ENDDO
        !
        CALL mp_sum( zmatilda, inter_image_comm ) 
1274
        CALL chi_invert_complex( zmatilda, zhead, zlambda, mypara%nglob, l_gammaq)
1275
1276
1277
1278
1279
        !
        DO ip = 1, pert%nloc
           glob_ip = pert%l2g(ip)
           z_epsm1_rfr_q(1:n_pdep_eigen_to_use,ip,ifreq,iq) = zlambda( 1:n_pdep_eigen_to_use, glob_ip)
        ENDDO 
1280
        IF( l_macropol .AND. l_gammaq ) z_head_rfr( ifreq ) = zhead
1281
1282
        !
     ENDDO
Marco Govoni's avatar
Marco Govoni committed
1283
1284
1285
1286
1287
     !
  ENDDO
  !
  DEALLOCATE( zlambda )
  DEALLOCATE( zmatilda )
1288
  DEALLOCATE( zmatr_q ) 
Marco Govoni's avatar
Marco Govoni committed
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
  !
  IF(l_generate_plot) THEN 
     CALL output_eps_head( )
  ENDIF
  !
END SUBROUTINE 
!
!
SUBROUTINE output_eps_head( )
  !
  USE kinds,                ONLY : DP
1300
  USE westcom,              ONLY : d_head_ifr,z_head_ifr,z_head_rfr,refreq_list,l_macropol,imfreq_list,wfreq_save_dir
Marco Govoni's avatar
Marco Govoni committed
1301
  USE constants,            ONLY : rytoev,fpi
Govoni's avatar
Govoni committed
1302
  !USE west_io,              ONLY : serial_table_output
Marco Govoni's avatar
Marco Govoni committed
1303
1304
1305
1306
1307
1308
1309
1310
  USE mp_world,             ONLY : mpime,root
  USE distribution_center,  ONLY : ifr,rfr
  USE mp,                   ONLY : mp_sum
  USE mp_global,            ONLY : intra_bgrp_comm 
  USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
  USE io_push,              ONLY : io_push_title,io_push_bar
  USE io_global,            ONLY : stdout
  USE cell_base,            ONLY : omega
Govoni's avatar
Govoni committed
1311
  USE json_module,          ONLY : json_file
Marco Govoni's avatar
Marco Govoni committed
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
  ! 
  IMPLICIT NONE
  !
  ! Workspace
  !
  CHARACTER(LEN=9) :: prefisso
  REAL(DP),ALLOCATABLE :: out_tabella(:,:)
  INTEGER :: ifreq,glob_ifreq
  REAL(DP) :: lf, rf, ep1, ep2, nn, kk, rr 
  TYPE(bar_type) :: barra
  REAL(DP) :: time_spent(2)
  REAL(DP),EXTERNAL :: get_clock
  CHARACTER(20),EXTERNAL :: human_readable_time
Govoni's avatar
Govoni committed
1325
1326
  TYPE(json_file) :: json
  INTEGER :: iunit
Marco Govoni's avatar
Marco Govoni committed
1327
  !
1328
  IF(l_macropol) THEN 
Marco Govoni's avatar
Marco Govoni committed
1329
1330
1331
1332
1333
     !
     CALL io_push_title("(O)ptics") 
     !
     CALL start_bar_type ( barra, 'optics', rfr%nloc )
     !
Govoni's avatar