solve_wfreq.f90 46.9 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
42
43
44
  USE mp_global,            ONLY : my_image_id,nimage,inter_image_comm,intra_bgrp_comm
  USE mp_world,             ONLY : mpime
  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
57
58
59
60
61
62
63
64
65
  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
  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
  USE distribution_center,  ONLY : pert,macropert,ifr,rfr
  USE class_idistribute,    ONLY : idistribute 
  USE wfreq_restart,        ONLY : solvewfreq_restart_write,solvewfreq_restart_read,bks_type
66
67
  USE class_bz_grid,        ONLY : bz_grid
  USE types_bz_grid,        ONLY : k_grid, q_grid
68
  USE chi_invert,           ONLY : chi_invert_real, chi_invert_complex
69
70
  USE class_coulomb,        ONLY : coulomb
  USE types_coulomb,        ONLY : pot3D
Marco Govoni's avatar
Marco Govoni committed
71
72
73
74
75
76
77
78
79
80
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  LOGICAL,INTENT(IN) :: l_read_restart,l_generate_plot
  !
  ! Workspace
  !
  INTEGER :: i1,i2,i3,im,ip,ig,glob_ip,ir,iv,iks,ipol,m
81
  CHARACTER(LEN=512)    :: fname
Marco Govoni's avatar
Marco Govoni committed
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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
  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(:,:)
  COMPLEX(DP),ALLOCATABLE :: phi(:,:)
  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 
     bks%max_ks        = nks 
     bks%min_ks        = 1 
  ENDIF
  !
  barra_load = 0
  DO iks = 1, nks
     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
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
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
  ! LOOP 
  !
  DO iks = 1, nks   ! KPOINT-SPIN
     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
     !
     IF(nks>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)
     !
     nbndval = nbnd_occ(iks)
     !
     mwo = - wk(iks) / omega
     zmwo = CMPLX( - wk(iks) / omega, 0._DP, KIND=DP)
     !
     bks%max_band = nbndval
     bks%min_band = 1
     !
     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
           !
           ! PHI 
           !
           ALLOCATE(phi(npwx*npol,3))
           ALLOCATE(phi_tmp(npwx*npol,3))
           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
           !
           DEALLOCATE( eprec )
           DEALLOCATE( e )
           DEALLOCATE( phi_tmp )
           !
        ENDIF
        !
        ! PSIC
        !
!        IF(gamma_only) THEN
           CALL single_invfft_gamma(dffts,npw,npwx,evc(1,iv),psic,'Wave') 
!        ELSEIF(noncolin) THEN
!           CALL SINGLEBAND_invfft_k(npw,evc(1     ,iv),npwx,psic_nc(1,1),dffts%nnr,.TRUE.)
!           CALL SINGLEBAND_invfft_k(npw,evc(1+npwx,iv),npwx,psic_nc(1,2),dffts%nnr,.TRUE.)
!        ELSE
!           CALL SINGLEBAND_invfft_k(npw,evc(1,iv),npwx,psic,dffts%nnr,.TRUE.)
!        ENDIF
        !
        ! ZEROS
        !
        dvpsi = 0._DP   
        !
        ! Read PDEP
        !
292
        ALLOCATE( pertg(npwqx) )
Marco Govoni's avatar
Marco Govoni committed
293
294
295
296
297
298
299
300
301
302
303
304
305
        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
              !
              WRITE(my_label_b,'(i6.6)') glob_ip
306
              fname = TRIM( wstat_save_dir ) // "/E"//TRIM(ADJUSTL(my_label_b))//".json"
Marco Govoni's avatar
Marco Govoni committed
307
308
309
              CALL pdep_read_G_and_distribute(fname,pertg)
              !
              ! Multiply by sqvc
310
              DO ig = 1, npwq
311
                 pertg(ig) = pot3D%sqvc(ig) * pertg(ig)
Marco Govoni's avatar
Marco Govoni committed
312
313
314
315
              ENDDO
              !
              ! Bring it to R-space
!              IF(gamma_only) THEN
316
                 CALL single_invfft_gamma(dffts,npwq,npwqx,pertg(1),pertr,TRIM(fftdriver))
Marco Govoni's avatar
Marco Govoni committed
317
318
319
320
321
                 DO ir=1,dffts%nnr 
                    pertr(ir)=psic(ir)*pertr(ir)
                 ENDDO
                 CALL single_fwfft_gamma(dffts,npw,npwx,pertr,dvpsi(1,ip),'Wave')
!              ELSEIF(noncolin) THEN
322
!                 CALL SINGLEBAND_invfft_k(npwq,pertg(1),npwx,pertr,dffts%nnr,.FALSE.)
Marco Govoni's avatar
Marco Govoni committed
323
324
325
326
!                 DO ir=1,dffts%nnr 
!                    pertr(ir)=psic_nc(ir,1)*pertr(ir)
!                 ENDDO
!                 CALL SINGLEBAND_fwfft_k(npw,pertr,dffts%nnr,dvpsi(1,ip),npwx,.TRUE.)
327
!                 CALL SINGLEBAND_invfft_k(npwq,pertg(1),npwx,pertr,dffts%nnr,.FALSE.)
Marco Govoni's avatar
Marco Govoni committed
328
329
330
331
332
!                 DO ir=1,dffts%nnr 
!                    pertr(ir)=psic_nc(ir,2)*pertr(ir)
!                 ENDDO
!                 CALL SINGLEBAND_fwfft_k(npw,pertr,dffts%nnr,dvpsi(1+npwx,ip),npwx,.TRUE.)
!              ELSE
333
!                 CALL SINGLEBAND_invfft_k(npwq,pertg(1),npwx,pertr,dffts%nnr,.FALSE.)
Marco Govoni's avatar
Marco Govoni committed
334
335
336
337
338
339
340
341
342
343
344
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
453
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
536
537
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
!                 DO ir=1,dffts%nnr 
!                    pertr(ir)=psic(ir)*pertr(ir)
!                 ENDDO
!                 CALL SINGLEBAND_fwfft_k(npw,pertr,dffts%nnr,dvpsi(1,ip),npwx,.TRUE.)
!              ENDIF 
              !
           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 ) )
           CALL get_brak_hyper_parallel(dvpsi,mypara%nloc,n_lanczos,q_s,braket,mypara%nloc,mypara%nlocx,mypara%nglob)
           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
     !
     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
565
  !
Marco Govoni's avatar
Marco Govoni committed
566
567
568
569
570
571
572
573
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
  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 
614
  USE westcom,              ONLY : west_prefix,n_pdep_eigen_to_use,n_lanczos,npwq,l_macropol,iks_l2g,z_epsm1_ifr_q,&
615
616
                                 &  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,&
617
                                 & ngq, igq_q
Marco Govoni's avatar
Marco Govoni committed
618
619
620
621
  USE mp_global,            ONLY : my_image_id,nimage,inter_image_comm,intra_bgrp_comm
  USE mp_world,             ONLY : mpime
  USE mp,                   ONLY : mp_bcast,mp_barrier,mp_sum
  USE io_global,            ONLY : stdout,ionode
622
  USE gvect,                ONLY : g,ngm,gstart
Marco Govoni's avatar
Marco Govoni committed
623
624
625
626
627
628
629
630
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
  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
  USE distribution_center,  ONLY : pert,macropert,ifr,rfr
  USE class_idistribute,    ONLY : idistribute 
642
  USE wfreq_restart,        ONLY : solvewfreq_restart_write,solvewfreq_restart_read,bksq_type
643
  USE class_bz_grid,        ONLY : bz_grid
644
  USE types_bz_grid,        ONLY : k_grid, q_grid, compute_phase
645
  USE chi_invert,           ONLY : chi_invert_complex
646
647
  USE class_coulomb,        ONLY : coulomb
  USE types_coulomb,        ONLY : pot3D
Marco Govoni's avatar
Marco Govoni committed
648
649
650
651
652
653
654
655
656
  !
  IMPLICIT NONE
  !
  ! I/O
  !
  LOGICAL,INTENT(IN) :: l_read_restart,l_generate_plot
  !
  ! Workspace
  !
657
  INTEGER :: i1,i2,i3,im,ip,ig,glob_ip,ir,iv,iks,ik,is,iq,ikqs,ipol,m
658
  CHARACTER(LEN=512)    :: fname
Marco Govoni's avatar
Marco Govoni committed
659
  CHARACTER(LEN=6)      :: my_label_b
660
  CHARACTER(LEN=5)      :: my_label_q
Marco Govoni's avatar
Marco Govoni committed
661
662
663
664
665
666
667
668
669
  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(:,:)
  COMPLEX(DP),ALLOCATABLE :: phi(:,:)
  COMPLEX(DP),ALLOCATABLE :: phi_tmp(:,:)
  COMPLEX(DP),ALLOCATABLE :: pertg(:),pertr(:)
670
671
672
673
  COMPLEX(DP),ALLOCATABLE :: evckpq(:,:)
  COMPLEX(DP),ALLOCATABLE :: psick(:),psick_nc(:,:)
  COMPLEX(DP),ALLOCATABLE :: phase(:)
  INTEGER :: npwkq
Marco Govoni's avatar
Marco Govoni committed
674
675
676
677
678
679
680
681
682
683
  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(:,:)
684
685
  COMPLEX(DP),ALLOCATABLE :: zmati_q(:,:,:,:)
  COMPLEX(DP),ALLOCATABLE :: zmatr_q(:,:,:,:)
Marco Govoni's avatar
Marco Govoni committed
686
  LOGICAL :: l_iks_skip, l_iv_skip
687
  LOGICAL :: l_gammaq
Marco Govoni's avatar
Marco Govoni committed
688
689
  REAL(DP) :: time_spent(2)
  REAL(DP),EXTERNAL :: get_clock
690
  TYPE(bksq_type) :: bksq
Marco Govoni's avatar
Marco Govoni committed
691
692
693
  REAL(DP),ALLOCATABLE :: eprec(:)
  INTEGER :: ierr
  REAL(DP),ALLOCATABLE :: e(:)
694
  REAL(DP) :: kpq(3), g0(3) 
Marco Govoni's avatar
Marco Govoni committed
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
  !
  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
  !
716
  ! ALLOCATE zmati_q, zmatr_q, where chi0 is stored
Marco Govoni's avatar
Marco Govoni committed
717
  !
718
719
  ALLOCATE( zmati_q( mypara%nglob, mypara%nloc, ifr%nloc, q_grid%np) )
  ALLOCATE( zmatr_q( mypara%nglob, mypara%nloc, rfr%nloc, q_grid%np) )
720
721
722
723
724
725
726
727
728
729
  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
730
731
  !
  IF(l_read_restart) THEN
732
     CALL solvewfreq_restart_read( bksq, zmati_q, zmatr_q, mypara%nglob, mypara%nloc )
Marco Govoni's avatar
Marco Govoni committed
733
  ELSE
734
735
736
737
738
739
740
741
742
743
     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
744
745
746
  ENDIF
  !
  barra_load = 0
747
  DO iq = 1, q_grid%np
748
     IF (iq<bksq%lastdone_q) CYCLE
749
     DO iks = 1, k_grid%nps
750
        IF(iq==bksq%lastdone_q .AND. iks<bksq%lastdone_ks) CYCLE
751
        DO iv = 1, nbnd_occ(iks)
752
           IF(iq==bksq%lastdone_q .AND. iks==bksq%lastdone_ks .AND. iv <= bksq%lastdone_band ) CYCLE
753
754
           barra_load = barra_load + 1
        ENDDO
Marco Govoni's avatar
Marco Govoni committed
755
756
757
758
759
760
761
762
763
764
765
766
     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 
  !
767
  DO iq = 1, q_grid%np   ! Q-POINT
768
     IF (iq<bksq%lastdone_q) CYCLE
Marco Govoni's avatar
Marco Govoni committed
769
     !
770
     npwq = ngq(iq)
771
     l_gammaq = q_grid%l_pIsGamma(iq)
Marco Govoni's avatar
Marco Govoni committed
772
     !
773
     CALL pot3D%init('Wave',.TRUE.,'default',iq)
Marco Govoni's avatar
Marco Govoni committed
774
     !
775
     DO iks = 1, k_grid%nps   ! KPOINT-SPIN
776
        !
777
        ik = k_grid%ip(iks) 
778
779
        is = k_grid%is(iks) 
        !
780
        IF(iq==bksq%lastdone_q .AND. iks<bksq%lastdone_ks) CYCLE
Marco Govoni's avatar
Marco Govoni committed
781
        !
782
        ! ... Set k-point, spin, kinetic energy, needed by Hpsi
Marco Govoni's avatar
Marco Govoni committed
783
        !
784
785
786
        current_k = iks
        IF ( lsda ) current_spin = isk(iks)
        call g2_kin( iks )
Marco Govoni's avatar
Marco Govoni committed
787
        !
788
        ! ... More stuff needed by the hamiltonian: nonlocal projectors
Marco Govoni's avatar
Marco Govoni committed
789
        !
790
        IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), k_grid%p_cart(1,ik), vkb )
791
        npw = ngk(iks)
Marco Govoni's avatar
Marco Govoni committed
792
        !
793
        ! ... read in wavefunctions from the previous iteration
Marco Govoni's avatar
Marco Govoni committed
794
        !
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
        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
820
        !
821
822
        !ikqs = kpq_grid%index_kq(iks,iq)
        !npwkq = ngk(ikqs)
Marco Govoni's avatar
Marco Govoni committed
823
        !
824
        CALL k_grid%find( k_grid%p_cart(:,ik) + q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )
Marco Govoni's avatar
Marco Govoni committed
825
        !
826
827
        npwkq = ngk(ikqs)
        !
828
        CALL compute_phase( g0, 'cart', phase )
829
830
831
        !
        ! Set wavefunctions at [k+q] in G space, for all bands,
        ! and store them in evckpq
Marco Govoni's avatar
Marco Govoni committed
832
        ! 
833
834
        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
835
        !
836
        nbndval = nbnd_occ(iks)
Marco Govoni's avatar
Marco Govoni committed
837
        !
838
839
        mwo = - k_grid%weight(iks) / omega
        zmwo = CMPLX( - k_grid%weight(iks) / omega, 0._DP, KIND=DP)
Marco Govoni's avatar
Marco Govoni committed
840
        !
841
842
        bksq%max_band = nbndval
        bksq%min_band = 1
Marco Govoni's avatar
Marco Govoni committed
843
        !
844
        ALLOCATE(dvpsi(npwx*npol,mypara%nlocx)) 
Marco Govoni's avatar
Marco Govoni committed
845
        !
846
        time_spent(1) = get_clock( 'wlanczos' ) 
Marco Govoni's avatar
Marco Govoni committed
847
        !
848
        ! LOOP over band states 
Marco Govoni's avatar
Marco Govoni committed
849
        !
850
        DO iv = 1, nbndval
851
           IF(iq==bksq%lastdone_q .AND. iks==bksq%lastdone_ks .AND. iv <= bksq%lastdone_band ) CYCLE
Marco Govoni's avatar
Marco Govoni committed
852
           !
853
           ! MACROPOL CASE
Marco Govoni's avatar
Marco Govoni committed
854
           !
855
           IF(l_macropol .AND. l_gammaq) THEN
Marco Govoni's avatar
Marco Govoni committed
856
              !
857
              ! PHI 
Marco Govoni's avatar
Marco Govoni committed
858
              !
859
860
861
862
863
864
865
866
867
868
              ALLOCATE(phi(npwx*npol,3))
              ALLOCATE(phi_tmp(npwx*npol,3))
              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
869
870
871
                 ENDDO
              ENDDO
              !
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
              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
              !
              DEALLOCATE( eprec )
              DEALLOCATE( e )
              DEALLOCATE( phi_tmp )
              !
              !
           ENDIF
Marco Govoni's avatar
Marco Govoni committed
899
           !
900
           ! PSIC
Marco Govoni's avatar
Marco Govoni committed
901
           !
902
903
           IF(noncolin) THEN
              CALL single_invfft_k(dffts,npwkq,npwx,evckpq(1     ,iv),psick_nc(1,1),'Wave',igk_k(1,ikqs))
904
              CALL single_invfft_k(dffts,npwkq,npwx,evckpq(npwx+1,iv),psick_nc(1,2),'Wave',igk_k(1,ikqs))
905
906
907
908
909
910
911
912
913
914
915
916
917
918
           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
919
              !
920
              glob_ip = mypara%l2g(ip)
Marco Govoni's avatar
Marco Govoni committed
921
              !
922
923
924
925
926
927
928
929
930
931
932
933
934
              ! Decide whether read dbs E or dhpi 
              !
              IF(glob_ip<=n_pdep_eigen_to_use) THEN
                 !
                 ! Exhume dbs eigenvalue
                 !
                 WRITE(my_label_b,'(i6.6)') glob_ip
                 WRITE(my_label_q,'(i5.5)') iq
                 fname = TRIM( wstat_save_dir ) // "/EQ"//TRIM(ADJUSTL(my_label_q))//"_"//TRIM(ADJUSTL(my_label_b))//".json"
                 CALL pdep_read_G_and_distribute(fname,pertg,iq)
                 !
                 ! Multiply by sqvc
                 DO ig = 1, npwq
935
                    pertg(ig) = pot3D%sqvc(ig) * pertg(ig)
Marco Govoni's avatar
Marco Govoni committed
936
                 ENDDO
937
938
939
940
941
942
943
944
                 !
                 ! 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))
945
                    CALL single_invfft_k(dffts,npwq,npwqx,pertg(1),pertr,'Wave',igq_q(1,iq))
946
947
948
                    DO ir=1,dffts%nnr 
                       pertr(ir)=phase(ir)*psick_nc(ir,2)*DCONJG(pertr(ir))
                    ENDDO
949
                    CALL single_fwfft_k(dffts,npw,npwx,pertr,dvpsi(npwx+1,ip),'Wave',igk_k(1,current_k))
950
951
952
953
954
955
956
957
958
959
                 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
                 !
960
961
962
963
                 IF (l_gammaq) THEN 
                    ipol = glob_ip-n_pdep_eigen_to_use
                    dvpsi(:,ip) = phi(:,ipol) * DSQRT(fpi * e2)
                 ENDIF
964
965
                 !
              ENDIF
Marco Govoni's avatar
Marco Govoni committed
966
              !
967
968
969
970
971
972
973
           ENDDO ! pert
           ! 
           DEALLOCATE(pertr)
           DEALLOCATE(pertg)
           IF(l_macropol.AND.l_gammaq) THEN
              DEALLOCATE(phi)
           ENDIF
Marco Govoni's avatar
Marco Govoni committed
974
           !
975
           CALL apply_alpha_pc_to_m_wfcs(nbndval,mypara%nloc,dvpsi,(1._DP,0._DP))
Marco Govoni's avatar
Marco Govoni committed
976
           !
977
           ! OVERLAP( glob_ip, im=1:nbnd ) = < psi_im iks | dvpsi_glob_ip >
Marco Govoni's avatar
Marco Govoni committed
978
           !
979
980
981
982
           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
983
           !
984
985
986
987
988
989
990
991
           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
992
993
           ENDDO
           !
994
995
           DEALLOCATE(ps_c)
           CALL mp_sum(overlap,inter_image_comm) 
Marco Govoni's avatar
Marco Govoni committed
996
           !
997
           ! Update zmati with cond
Marco Govoni's avatar
Marco Govoni committed
998
999
1000
1001
1002
           !
           DO ifreq = 1, ifr%nloc
              !
              frequency = imfreq_list( ifreq ) 
              !
1003
1004
1005
1006
1007
              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
1008
1009
1010
1011
1012
                 !
                 DO ip = 1, mypara%nloc
                    glob_ip = mypara%l2g(ip)
                    DO glob_jp = 1, mypara%nglob
                       !
1013
1014
                       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
1015
1016
1017
1018
                       !
                    ENDDO
                 ENDDO
                 !
1019
              ENDDO ! ic
Marco Govoni's avatar
Marco Govoni committed
1020
1021
           ENDDO ! ifreq
           !
1022
           ! Update zmatr with cond
Marco Govoni's avatar
Marco Govoni committed
1023
1024
1025
1026
1027
           !
           DO ifreq = 1, rfr%nloc
              !
              frequency = refreq_list( ifreq ) 
              !
1028
1029
1030
1031
1032
1033
              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
1034
1035
1036
1037
1038
                 !
                 DO ip = 1, mypara%nloc
                    glob_ip = mypara%l2g(ip)
                    DO glob_jp = 1, mypara%nglob
                       !
1039
1040
                       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
1041
1042
1043
1044
                       !
                    ENDDO
                 ENDDO
                 !
1045
1046
              ENDDO ! ic
           ENDDO ! ifreq 
Marco Govoni's avatar
Marco Govoni committed
1047
           !
1048
           DEALLOCATE(overlap)
Marco Govoni's avatar
Marco Govoni committed
1049
           !
1050
           ! Apply Pc, to be sure
Marco Govoni's avatar
Marco Govoni committed
1051
           !
1052
           CALL apply_alpha_pc_to_m_wfcs(nbnd,mypara%nloc,dvpsi,(1._DP,0._DP))
Marco Govoni's avatar
Marco Govoni committed
1053
           !
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
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
           ! 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 ) )
              CALL get_brak_hyper_parallel_complex(dvpsi,mypara%nloc,n_lanczos,q_s,braket,mypara%nloc,mypara%nlocx,mypara%nglob)
              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
1137
1138
1139
1140
1141
1142
1143
                 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
1144
1145
                 time_spent(1) = get_clock( 'wlanczos' )
              ENDIF
Marco Govoni's avatar
Marco Govoni committed
1146
           ENDIF
1147
1148
1149
1150
           !
           CALL update_bar_type( barra, 'wlanczos', 1 )
           !
        ENDDO ! BANDS
Marco Govoni's avatar
Marco Govoni committed
1151
        !
1152
        DEALLOCATE(dvpsi)
Marco Govoni's avatar
Marco Govoni committed
1153
        !
1154
     ENDDO ! KPOINT-SPIN 
Marco Govoni's avatar
Marco Govoni committed
1155
     !
1156
  ENDDO ! QPOINT
Marco Govoni's avatar
Marco Govoni committed
1157
1158
1159
  !
  CALL stop_bar_type( barra, 'wlanczos' )
  !
1160
1161
1162
  DEALLOCATE( evckpq )
  DEALLOCATE( phase )
  !
Marco Govoni's avatar
Marco Govoni committed
1163
1164
1165
1166
  ! EPS-1 imfreq
  !
  ALLOCATE( zmatilda( mypara%nglob, mypara%nglob ) )
  ALLOCATE( zlambda( n_pdep_eigen_to_use, n_pdep_eigen_to_use) ) 
1167
  ALLOCATE( z_epsm1_ifr_q( pert%nglob, pert%nloc, ifr%nloc, q_grid%np) )
1168
  z_epsm1_ifr_q = 0._DP
Marco Govoni's avatar
Marco Govoni committed
1169
1170
1171
1172
  IF(l_macropol) ALLOCATE( z_head_ifr( ifr%nloc) )
  !
  CALL start_clock('chi0_diag')
  !
1173
  DO iq = 1, q_grid%np
Marco Govoni's avatar
Marco Govoni committed
1174
     !
1175
     l_gammaq = q_grid%l_pIsGamma(iq)
Marco Govoni's avatar
Marco Govoni committed
1176
     !
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
     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 )
        ! 
1187
        CALL chi_invert_complex( zmatilda, zhead, zlambda, mypara%nglob, l_gammaq)
1188
1189
1190
1191
1192
1193
1194
1195
        !
        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
1196
1197
     !
  ENDDO
1198
  !
Marco Govoni's avatar
Marco Govoni committed
1199
1200
1201
1202
  CALL stop_clock('chi0_diag')
  !
  DEALLOCATE( zlambda )
  DEALLOCATE( zmatilda )
1203
  DEALLOCATE( zmati_q )
Marco Govoni's avatar
Marco Govoni committed
1204
1205
1206
1207
1208
  !
  ! EPS-1 refreq
  !
  ALLOCATE( zmatilda( mypara%nglob, mypara%nglob ) )
  ALLOCATE( zlambda( n_pdep_eigen_to_use, n_pdep_eigen_to_use ) ) 
1209
  ALLOCATE( z_epsm1_rfr_q( pert%nglob, pert%nloc, rfr%nloc, q_grid%np) )
1210
  z_epsm1_rfr_q = 0._DP
Marco Govoni's avatar
Marco Govoni committed
1211
1212
  IF(l_macropol) ALLOCATE( z_head_rfr( rfr%nloc) )
  !
1213
  DO iq = 1, q_grid%np
Marco Govoni's avatar
Marco Govoni committed
1214
     !
1215
     l_gammaq = q_grid%l_pIsGamma(iq)
Marco Govoni's avatar
Marco Govoni committed
1216
     !
1217
1218
1219
1220
1221
1222
1223
1224
1225
     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 ) 
1226
        CALL chi_invert_complex( zmatilda, zhead, zlambda, mypara%nglob, l_gammaq)
1227
1228
1229
1230
1231
        !
        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 
1232
        IF( l_macropol .AND. l_gammaq ) z_head_rfr( ifreq ) = zhead
1233
1234
        !
     ENDDO
Marco Govoni's avatar
Marco Govoni committed
1235
1236
1237
1238
1239
     !
  ENDDO
  !
  DEALLOCATE( zlambda )
  DEALLOCATE( zmatilda )
1240
  DEALLOCATE( zmatr_q ) 
Marco Govoni's avatar
Marco Govoni committed
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
  !
  IF(l_generate_plot) THEN 
     CALL output_eps_head( )
  ENDIF
  !
END SUBROUTINE 
!
!
SUBROUTINE output_eps_head( )
  !
  USE kinds,                ONLY : DP
1252
  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
1253
  USE constants,            ONLY : rytoev,fpi
Govoni's avatar
Govoni committed
1254
  !USE west_io,              ONLY : serial_table_output
Marco Govoni's avatar
Marco Govoni committed
1255
1256
1257
1258
1259
1260
1261
1262
  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
1263
  USE json_module,          ONLY : json_file
Marco Govoni's avatar
Marco Govoni committed
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
  ! 
  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
1277
1278
  TYPE(json_file) :: json
  INTEGER :: iunit
Marco Govoni's avatar
Marco Govoni committed
1279
  !
1280
  IF(l_macropol) THEN 
Marco Govoni's avatar
Marco Govoni committed
1281
1282
1283
1284
1285
     !
     CALL io_push_title("(O)ptics") 
     !
     CALL start_bar_type ( barra, 'optics', rfr%nloc )
     !
Govoni's avatar
Govoni committed
1286
1287
     time_spent(1) = get_clock( 'optics' ) 
     !
Marco Govoni's avatar
Marco Govoni committed
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
     ! head_rfr
     !
     ALLOCATE( out_tabella(rfr%nglob,8) )
     !
     out_tabella = 0._DP
     DO ifreq = 1, rfr%nloc 
        glob_ifreq = rfr%l2g(ifreq)
        rf =  REAL( z_head_rfr( ifreq ), KIND=DP ) + 1._DP
        lf = -AIMAG( z_head_rfr( ifreq ) )
        ep1 = rf / (rf**2+lf**2)
        ep2 = lf / (rf**2+lf**2)
        nn = DSQRT( 0.5_DP * ( ep1 + DSQRT( ep1**2 + ep2**2 ) ) ) 
        kk = DSQRT( 0.5_DP * (-ep1 + DSQRT( ep1**2 + ep2**2 ) ) )
        rr = ((1._DP-nn)**2+kk**2)/((1._DP+nn)**2+kk**2) 
        !
        out_tabella( glob_ifreq, 1 ) = refreq_list( ifreq )*rytoev
        out_tabella( glob_ifreq, 2 ) = ep1
        out_tabella( glob_ifreq, 3 ) = ep2
        out_tabella( glob_ifreq, 4 ) = lf
        out_tabella( glob_ifreq, 5 ) = nn
        out_tabella( glob_ifreq, 6 ) = kk
        out_tabella( glob_ifreq, 7 ) = rr
        out_tabella( glob_ifreq, 8 ) = (ep1-1._DP)*3._DP*omega/fpi/(ep1+2._DP)
        !
        CALL update_bar_type( barra, 'optics', 1 )
        !
     ENDDO
     !
Govoni's avatar
Govoni committed
1316
1317
     CALL stop_bar_type( barra, 'optics' )
     !
Marco Govoni's avatar
Marco Govoni committed
1318
1319
     CALL mp_sum( out_tabella, intra_bgrp_comm ) 
     !
Govoni's avatar
Govoni committed
1320
1321
1322
     !CALL serial_table_output(mpime==root,4000,'optics',out_tabella,&
     !& rfr%nglob,8,&
     !& (/'     E[eV]','      eps1','      eps2','      EELF','         n','         k','      Refl',' pol[au^3]'/))
Marco Govoni's avatar
Marco Govoni committed
1323
     !
Govoni's avatar
Govoni committed
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
     IF( mpime == root ) THEN 
        !
        CALL json%initialize()
        !
        CALL json%add("e",out_tabella(:,1))
        CALL json%add("eps1",out_tabella(:,2))
        CALL json%add("eps2",out_tabella(:,3))
        CALL json%add("EELF",out_tabella(:,4))
        CALL json%add("n",out_tabella(:,5))
        CALL json%add("k",out_tabella(:,6))
        CALL json%add("refl",out_tabella(:,7))
        CALL json%add("pol",out_tabella(:,8))
        !
        OPEN( NEWUNIT=iunit, FILE=TRIM(wfreq_save_dir)//"/optics.json" )
        CALL json%print_file( iunit )
        CLOSE( iunit )
        !
        CALL json%destroy()
        !
     ENDIF
Marco Govoni's avatar
Marco Govoni committed
1344
1345
1346
1347
1348
     !
     time_spent(2) = get_clock( 'optics' ) 
     !
     WRITE(stdout,'(  5x," ")')
     CALL io_push_bar()
Govoni's avatar
Govoni committed
1349
1350
     WRITE(stdout, "(5x, 'File ',a,' written in ',a20)") TRIM(wfreq_save_dir)//"/optics.json",&
     & human_readable_time(time_spent(2)-time_spent(1))
Marco Govoni's avatar
Marco Govoni committed
1351
     CALL io_push_bar()
Govoni's avatar
Govoni committed
1352
1353
1354
     ! 
     DEALLOCATE( out_tabella )
     !
Marco Govoni's avatar
Marco Govoni committed
1355
1356
1357
  ENDIF
  !
END SUBROUTINE