dfpt_module.f90 24.8 KB
Newer Older
1
!
Marco Govoni's avatar
Marco Govoni committed
2
! Copyright (C) 2015-2021 M. Govoni
3
4
5
6
7
8
9
! 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.
!
10
! Contributors to this file:
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
! Marco Govoni
!
!-----------------------------------------------------------------------
MODULE dfpt_module
  !-----------------------------------------------------------------------
  !
  IMPLICIT NONE
  !
  CONTAINS
    !
    !-------------------------------------------------------------------------
    SUBROUTINE dfpt (m,dvg,dng,tr2,iq)
      !-----------------------------------------------------------------------
      !
      USE kinds,                 ONLY : DP
      USE io_global,             ONLY : stdout
27
28
29
      USE wvfct,                 ONLY : nbnd,et
      USE fft_base,              ONLY : dffts
      USE gvect,                 ONLY : gstart
30
31
      USE wavefunctions_module,  ONLY : evc,psic
      USE mp,                    ONLY : mp_sum,mp_barrier,mp_bcast
32
      USE mp_global,             ONLY : inter_image_comm,inter_pool_comm,my_image_id,inter_bgrp_comm
33
34
35
36
37
      USE fft_at_k,              ONLY : single_fwfft_k,single_invfft_k
      USE fft_at_gamma,          ONLY : single_fwfft_gamma,single_invfft_gamma,double_fwfft_gamma,double_invfft_gamma
      USE buffers,               ONLY : get_buffer
      USE noncollin_module,      ONLY : noncolin,npol
      USE bar,                   ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
38
      USE pwcom,                 ONLY : current_spin,nelup,neldw,isk,xk,npw,npwx,lsda,&
39
                                      & current_k,ngk,igk_k
40
41
42
      USE cell_base,             ONLY : omega
      USE control_flags,         ONLY : gamma_only
      USE uspp,                  ONLY : nkb, vkb
43
44
      USE westcom,               ONLY : nbnd_occ,iuwfc,lrwfc,npwqx,npwq,igq_q,fftdriver
      USE io_push,               ONLY : io_push_title
45
      USE mp_world,              ONLY : world_comm
46
      USE types_bz_grid,         ONLY : k_grid, q_grid, compute_phase
47
48
      USE class_idistribute,     ONLY : idistribute
      USE distribution_center,   ONLY : occband
49
50
51
52
53
54
55
56
57
58
59
60
61
      !
      IMPLICIT NONE
      !
      ! I/O
      !
      INTEGER, INTENT(IN), OPTIONAL :: iq
      INTEGER, INTENT(IN) :: m
      COMPLEX(DP), INTENT(IN) :: dvg(npwqx,m)
      COMPLEX(DP), INTENT(OUT) :: dng(npwqx,m)
      REAL(DP),INTENT(IN) :: tr2
      !
      ! Workspace
      !
62
      INTEGER :: ipert, ig, ir, ibnd, ibnd2, lbnd, iks, ikqs, ikq, ik, is
63
64
65
66
      INTEGER :: nbndval, ierr
      INTEGER :: npwkq
      !
      REAL(DP) :: g0(3)
Marco Govoni's avatar
Marco Govoni committed
67
      REAL(DP) :: anorm
68
      REAL(DP), ALLOCATABLE :: eprec(:)
69
70
      REAL(DP), ALLOCATABLE :: eprec_loc(:)
      REAL(DP), ALLOCATABLE :: et_loc(:)
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
      !
      COMPLEX(DP), ALLOCATABLE :: dvpsi(:,:),dpsi(:,:)
      COMPLEX(DP), ALLOCATABLE :: aux_r(:),aux_g(:)
      COMPLEX(DP), ALLOCATABLE :: dpsic(:)
      !
      COMPLEX(DP), ALLOCATABLE :: evckmq(:,:)
      COMPLEX(DP), ALLOCATABLE :: phase(:)
      !
      TYPE(bar_type) :: barra
      !
      LOGICAL :: l_dost
      !
      CHARACTER(LEN=512) :: title
      !
      CALL mp_barrier( world_comm )
      !
      CALL report_dynamical_memory( )
      !
      l_dost = ( tr2 >= 0._DP )
      !
      IF( l_dost ) THEN
         WRITE(title,'(a,es14.6)') "Sternheimer eq. solver... with threshold = ", tr2
      ELSE
         WRITE(title,'(a,es14.6)') "Sternheimer eq. solver... with lite-solver"
      ENDIF
      CALL io_push_title(TRIM(ADJUSTL(title)))
      !
98
99
      occband = idistribute()
      !
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
      dng=0.0_DP
      !
      CALL start_bar_type( barra, 'dfpt', MAX(m,1) * k_grid%nps )
      !
      IF (.NOT.gamma_only) THEN
         ALLOCATE( evckmq(npwx*npol,nbnd) )
         ALLOCATE( phase(dffts%nnr) )
      ENDIF
      !
      DO iks = 1, k_grid%nps  ! KPOINT-SPIN LOOP
         !
         ik = k_grid%ip(iks)
         is = k_grid%is(iks)
         !
         ! ... Set k-point, spin, kinetic energy, needed by Hpsi
         !
         current_k = iks
         IF ( lsda ) current_spin = isk(iks)
118
         CALL g2_kin(iks)
119
120
121
122
123
124
125
         !
         ! ... More stuff needed by the hamiltonian: nonlocal projectors
         !
         IF ( nkb > 0 ) CALL init_us_2( ngk(iks), igk_k(1,iks), k_grid%p_cart(1,ik), vkb )
         !
         nbndval = nbnd_occ(iks)
         !
126
127
         CALL occband%init( nbndval, 'b', 'occband', .FALSE. )
         !
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
         ! ... Number of G vectors for PW expansion of wfs at k
         !
         npw = ngk(iks)
         !
         ! ... Read wavefuctions at k in G space, for all bands, and store them in evc
         !
         IF(k_grid%nps>1) THEN
            IF ( my_image_id == 0 ) CALL get_buffer( evc, lrwfc, iuwfc, iks )
            CALL mp_bcast( evc, 0, inter_image_comm )
         ENDIF
         !
         IF (gamma_only) THEN
            !
            ikqs = iks
            g0 = 0.0_DP
            !
         ELSE
            !
            ! ... Find G0 and compute phase
            !
148
            !CALL k_grid%find( k_grid%p_cart(:,ik) - q_grid%p_cart(:,iq), is, 'cart', ikqs, g0 )   !M
149
150
            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
151
            !
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
            CALL compute_phase( g0, 'cart', phase )
            !
            ! ... Number of G vectors for PW expansion of wfs at [k-q]
            !
            npwkq = ngk(ikqs)
            !
            ! ... Set wavefunctions at [k-q] in G space, for all bands, and store them in evckmq
            !
            IF ( my_image_id == 0 ) CALL get_buffer( evckmq, lrwfc, iuwfc, ikqs )
            CALL mp_bcast( evckmq, 0, inter_image_comm )
            !
         ENDIF
         !
         !
         ALLOCATE(eprec(nbndval))
167
168
         ALLOCATE(eprec_loc(occband%nloc))
         ALLOCATE(et_loc(occband%nloc))
169
170
         CALL set_eprec(nbndval,evc(1,1),eprec)
         !
171
172
173
174
175
176
177
178
         DO lbnd = 1,occband%nloc
            ibnd = occband%l2g(lbnd)
            eprec_loc(lbnd) = eprec(ibnd)
            et_loc(lbnd) = et(ibnd,ikqs)
         ENDDO
         !
         ALLOCATE(dvpsi(npwx*npol,occband%nloc))
         ALLOCATE(dpsi(npwx*npol,occband%nloc))
179
180
181
         !
         DO ipert = 1, m
            !
Marco Govoni's avatar
Marco Govoni committed
182
183
            ALLOCATE( aux_g(npwqx) ); aux_g = 0._DP
            ALLOCATE( aux_r(dffts%nnr) ); aux_r = 0._DP
184
            !
185
186
            DO CONCURRENT (ig = 1:npwq)
               aux_g(ig) = dvg(ig,ipert)
187
188
189
190
191
            ENDDO
            !
            ! ... inverse Fourier transform of the perturbation: (q+)G ---> R
            !
            IF (gamma_only) THEN
192
               CALL single_invfft_gamma(dffts,npwq,npwqx,aux_g,aux_r,TRIM(fftdriver))
193
194
195
196
197
198
199
200
201
202
203
204
            ELSE
               CALL single_invfft_k(dffts,npwq,npwqx,aux_g,aux_r,'Wave',igq_q(1,iq))
            ENDIF
            !
            ! The perturbation is in aux_r
            !
            dvpsi=0._DP
            dpsi =0._DP
            !
            IF(gamma_only) THEN
               !
               ! double bands @ gamma
205
206
207
208
               DO lbnd = 1,occband%nloc-MOD(occband%nloc,2),2
                  !
                  ibnd = occband%l2g(lbnd)
                  ibnd2 = occband%l2g(lbnd+1)
209
                  !
210
                  CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),evc(1,ibnd2),psic,'Wave')
Marco Govoni's avatar
Marco Govoni committed
211
                  DO CONCURRENT (ir=1:dffts%nnr)
212
213
                     psic(ir) = psic(ir) * REAL(aux_r(ir),KIND=DP)
                  ENDDO
214
                  CALL double_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,lbnd),dvpsi(1,lbnd+1),'Wave')
215
216
                  !
               ENDDO
217
               !
218
               ! single band @ gamma
219
220
221
222
               IF( MOD(occband%nloc,2) == 1 ) THEN
                  !
                  lbnd = occband%nloc
                  ibnd = occband%l2g(lbnd)
223
224
                  !
                  CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),psic,'Wave')
Marco Govoni's avatar
Marco Govoni committed
225
                  DO CONCURRENT (ir=1:dffts%nnr)
226
227
                     psic(ir) = CMPLX( REAL(psic(ir),KIND=DP) * REAL(aux_r(ir),KIND=DP), 0._DP, KIND=DP)
                  ENDDO
228
                  CALL single_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,lbnd),'Wave')
229
230
231
232
233
                  !
               ENDIF
               !
            ELSE
               !
234
235
236
               DO lbnd = 1,occband%nloc
                  !
                  ibnd = occband%l2g(lbnd)
237
238
239
240
241
242
243
244
                  !
                  ! ... inverse Fourier transform of wfs at [k-q]: (k-q+)G ---> R
                  !
                  CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1,ibnd),psic,'Wave',igk_k(1,ikqs))
                  !
                  ! ... construct right-hand-side term of Sternheimer equation:
                  ! ... product of wavefunction at [k-q], phase and perturbation in real space
                  !
Marco Govoni's avatar
Marco Govoni committed
245
                  DO CONCURRENT (ir = 1:dffts%nnr)
246
247
248
249
250
251
                     psic(ir) = psic(ir) * phase(ir) * aux_r(ir)
                  ENDDO
                  !
                  ! Fourier transform product of wf at [k-q], phase and
                  ! perturbation of wavevector q: R ---> (k+)G
                  !
252
                  CALL single_fwfft_k(dffts,npw,npwx,psic,dvpsi(1,lbnd),'Wave',igk_k(1,iks))
253
254
255
256
257
                  !
                  ! dv|psi> is in dvpsi
                  !
               ENDDO
               !
258
               IF (noncolin) THEN
259
260
261
                  DO lbnd = 1,occband%nloc
                     !
                     ibnd = occband%l2g(lbnd)
262
263
264
                     !
                     CALL single_invfft_k(dffts,npwkq,npwx,evckmq(npwx+1,ibnd),psic,'Wave',igk_k(1,ikqs))
                     !
Marco Govoni's avatar
Marco Govoni committed
265
                     DO CONCURRENT (ir = 1:dffts%nnr)
266
267
268
                        psic(ir) = psic(ir) * phase(ir) * aux_r(ir)
                     ENDDO
                     !
269
                     CALL single_fwfft_k(dffts,npw,npwx,psic,dvpsi(npwx+1,lbnd),'Wave',igk_k(1,iks))
270
271
272
273
274
275
276
277
278
279
280
                     !
                  ENDDO
               ENDIF
               !
            ENDIF
            !
            DEALLOCATE( aux_g )
            DEALLOCATE( aux_r )
            !
            ! - P_c | dvpsi >
            !
281
            CALL apply_alpha_pc_to_m_wfcs( nbndval, occband%nloc, dvpsi, (-1._DP,0._DP) )
282
            !
283
            CALL precondition_m_wfcts( occband%nloc, dvpsi, dpsi, eprec_loc )
284
285
            !
            IF( l_dost) THEN
286
               !
287
288
289
290
               ! The Sternheimer operator is (H_k - E_(k-q) + alpha * P_v)
               ! The Hamiltonian is evaluated at the k-point current_k in h_psi
               ! (see also PHonon/PH/cch_psi_all.f90, where H_(k+q) is evaluated)
               !
291
               CALL linsolve_sternheimer_m_wfcts (nbndval, occband%nloc, dvpsi, dpsi, et_loc, eprec_loc, tr2, ierr )
292
293
294
295
296
297
298
299
300
301
302
303
304
305
               !
               IF(ierr/=0) THEN
                  WRITE(stdout, '(7X,"** WARNING : PERT ",i8," iks ",I8," not converged, ierr = ",i8)') ipert,iks,ierr
               ENDIF
               !
            ENDIF
            !
            ALLOCATE( aux_r(dffts%nnr) )
            !
            aux_r=0._DP
            !
            IF(gamma_only) THEN
               !
               ! double band @ gamma
306
307
308
               DO lbnd = 1,occband%nloc
                  !
                  ibnd = occband%l2g(lbnd)
309
                  !
310
                  CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),dpsi(1,lbnd),psic,'Wave')
Marco Govoni's avatar
Marco Govoni committed
311
312
                  DO CONCURRENT (ir=1:dffts%nnr)
                     aux_r(ir) = aux_r(ir) + CMPLX( REAL( psic(ir),KIND=DP) * DIMAG( psic(ir)) , 0.0_DP, KIND=DP)
313
314
315
316
317
318
319
320
                  ENDDO
                  !
               ENDDO
               !
            ELSE
               !
               ALLOCATE( dpsic(dffts%nnr) )
               !
321
322
323
               DO lbnd = 1,occband%nloc
                  !
                  ibnd = occband%l2g(lbnd)
324
325
326
327
328
329
330
                  !
                  ! inverse Fourier transform of wavefunction at [k-q]: (k-q+)G ---> R
                  !
                  CALL single_invfft_k(dffts,npwkq,npwx,evckmq(1,ibnd),psic,'Wave',igk_k(1,ikqs))
                  !
                  ! inverse Fourier transform of perturbed wavefunction: (k+)G ---> R
                  !
331
                  CALL single_invfft_k(dffts,npw,npwx,dpsi(1,lbnd),dpsic,'Wave',igk_k(1,iks))
332
                  !
333
                  DO CONCURRENT (ir = 1: dffts%nnr)
334
335
336
337
338
                     aux_r(ir) = aux_r(ir) + CONJG( psic(ir) * phase(ir) ) * dpsic(ir)
                  ENDDO
                  !
               ENDDO
               !
339
               IF (noncolin) THEN
340
341
342
                  DO lbnd = 1,occband%nloc
                     !
                     ibnd = occband%l2g(lbnd)
343
                     !
Matteo Gerosa's avatar
Matteo Gerosa committed
344
                     CALL single_invfft_k(dffts,npwkq,npwx,evckmq(npwx+1,ibnd),psic,'Wave',igk_k(1,ikqs))
345
                     !
346
                     CALL single_invfft_k(dffts,npw,npwx,dpsi(npwx+1,lbnd),dpsic,'Wave',igk_k(1,iks))
347
                     !
Marco Govoni's avatar
Marco Govoni committed
348
                     DO CONCURRENT (ir = 1: dffts%nnr)
349
350
351
352
353
354
355
356
357
358
                        aux_r(ir) = aux_r(ir) + CONJG( psic(ir) * phase(ir) ) * dpsic(ir)
                     ENDDO
                     !
                  ENDDO
               ENDIF
               !
               DEALLOCATE( dpsic )
               !
            ENDIF
            !
359
360
361
362
            ! Sum up aux_r from band groups
            !
            CALL mp_sum(aux_r,inter_bgrp_comm)
            !
363
364
365
            ! The perturbation is in aux_r
            !
            ALLOCATE( aux_g(npwqx) )
366
            !
367
368
369
370
371
372
            IF(gamma_only) THEN
               CALL single_fwfft_gamma(dffts,npwq,npwqx,aux_r,aux_g,TRIM(fftdriver))
            ELSE
               CALL single_fwfft_k(dffts,npwq,npwqx,aux_r,aux_g,'Wave',igq_q(1,iq))
            ENDIF
            !
Marco Govoni's avatar
Marco Govoni committed
373
374
            DO CONCURRENT( ig = 1: npwq) ! pert acts only on body
               dng(ig,ipert) = dng(ig,ipert) + 2._DP * k_grid%weight(iks) * aux_g(ig) / omega
375
376
377
378
379
            ENDDO
            !
            DEALLOCATE( aux_g )
            DEALLOCATE( aux_r )
            !
380
            CALL update_bar_type( barra, 'dfpt', 1 )
381
382
383
384
385
386
            !
         ENDDO ! ipert
         !
         IF( m == 0 ) CALL update_bar_type( barra, 'dfpt', 1 )
         !
         DEALLOCATE( eprec )
387
388
         DEALLOCATE( eprec_loc )
         DEALLOCATE( et_loc )
389
390
391
392
393
394
395
         DEALLOCATE( dpsi )
         DEALLOCATE( dvpsi )
         !
      ENDDO ! K-POINT and SPIN
      !
      IF ( gamma_only ) THEN
         IF ( gstart == 2 ) dng(1,1:m) = CMPLX( 0._DP, 0._DP, KIND=DP )
396
      ELSE
397
398
399
400
401
402
403
404
405
406
407
408
         IF ( gstart == 2 .AND. q_grid%l_pIsGamma(iq) ) dng(1,1:m) = CMPLX( 0._DP, 0._DP, KIND=DP )
         DEALLOCATE( evckmq )
         DEALLOCATE( phase )
      ENDIF
      !
      CALL mp_sum(dng,inter_pool_comm)
      !
      CALL mp_barrier( world_comm )
      !
      CALL stop_bar_type( barra, 'dfpt' )
      !
    END SUBROUTINE
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
END MODULE
!!-----------------------------------------------------------------------
!SUBROUTINE dfpt (m,dvg,dng,tr2)
!  !-----------------------------------------------------------------------
!  !
!  USE kinds,                 ONLY : DP
!  USE io_global,             ONLY : stdout
!  USE wvfct,                 ONLY : nbnd,g2kin,et
!  USE fft_base,              ONLY : dfftp,dffts
!  USE gvect,                 ONLY : nl,gstart,g,ngm
!  USE wavefunctions_module,  ONLY : evc,psic
!  USE gvecw,                 ONLY : gcutw
!  USE mp,                    ONLY : mp_sum,mp_barrier,mp_bcast
!  USE mp_global,             ONLY : inter_image_comm,inter_pool_comm,my_image_id
!  USE fft_at_gamma,          ONLY : single_fwfft_gamma,single_invfft_gamma,double_fwfft_gamma,double_invfft_gamma
!  USE fft_at_k,              ONLY : single_fwfft_k,single_invfft_k
!  USE buffers,               ONLY : get_buffer
!  USE noncollin_module,      ONLY : noncolin,npol
!  USE bar,                   ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
!  USE pwcom,                 ONLY : current_spin,wk,nks,nelup,neldw,isk,igk_k,xk,npw,npwx,lsda,nkstot,&
!                                  & current_k,ngk
!  USE cell_base,             ONLY : tpiba2,omega
!  USE control_flags,         ONLY : gamma_only, io_level
!  USE io_files,              ONLY : tmp_dir, nwordwfc, iunwfc, diropn
!  USE uspp,                  ONLY : nkb, vkb, okvan
!  USE constants,             ONLY : e2,fpi
!  USE westcom,               ONLY : npwq,nbnd_occ,iuwfc,lrwfc,npwqx,fftdriver
!  USE io_push,               ONLY : io_push_title
!  USE mp_world,              ONLY : mpime,world_comm
!  USE types_coulomb,         ONLY : pot3D
!  !
!  IMPLICIT NONE
!  !
!  ! I/O
!  !
!  INTEGER,INTENT(IN) :: m
!  COMPLEX(DP),INTENT(IN) :: dvg(npwqx,m)
!  COMPLEX(DP),INTENT(OUT) :: dng(npwqx,m)
!  REAL(DP),INTENT(IN) :: tr2
!  !
!  ! Workspace
!  !
!  INTEGER :: ipert, ig, ir, ibnd, iks
!  INTEGER :: nbndval, ierr
!  !
455
!  REAL(DP) :: anorm, prod
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
!  REAL(DP),ALLOCATABLE :: eprec(:)
!  !
!  COMPLEX(DP),ALLOCATABLE :: dvpsi(:,:),dpsi(:,:)
!  COMPLEX(DP),ALLOCATABLE :: aux_r(:),aux_g(:)
!  COMPLEX(DP),ALLOCATABLE :: dpsic(:)
!  !
!  TYPE(bar_type) :: barra
!  !
!  LOGICAL :: conv_dfpt
!  LOGICAL :: exst,exst_mem
!  LOGICAL :: l_dost
!  !
!  CHARACTER(LEN=512) :: title
!  !
!  CALL mp_barrier( world_comm )
!  !
!  CALL report_dynamical_memory()
!  !
474
!  l_dost = ( tr2 >= 0._DP )
475
!  !
476
477
!  IF( l_dost ) THEN
!     WRITE(title,'(a,es14.6)') "Sternheimer eq. solver... with threshold = ", tr2
478
479
480
481
482
483
484
!  ELSE
!     WRITE(title,'(a,es14.6)') "Sternheimer eq. solver... with lite-solver"
!  ENDIF
!  CALL io_push_title(TRIM(ADJUSTL(title)))
!  !
!  dng=0.0_DP
!  !
485
486
!  CALL start_bar_type( barra, 'dfpt', MAX(m,1) * nks )
!  !IF(nks>1) CALL diropn(iuwfc,'wfc',lrwfc,exst)
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
!  !
!  DO iks = 1, nks  ! KPOINT-SPIN LOOP
!     !
!     ! ... 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
505
!        !lrwfc = nbnd * npwx * npol
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
!        !!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)
!     IF(nbndval==0) THEN
!        CALL update_bar_type( barra,'dfpt', MAX(1,m) )
!        CYCLE
!     ENDIF
!     !
!     ALLOCATE(eprec(nbndval))
!     CALL set_eprec(nbndval,evc(1,1),eprec)
!     !
538
!     ALLOCATE(dvpsi(npwx*npol,nbndval))
539
540
541
542
543
544
545
546
547
548
549
!     ALLOCATE(dpsi(npwx*npol,nbndval))
!     !
!     DO ipert = 1, m
!        !
!        ALLOCATE(aux_g(npwqx))
!        ALLOCATE(aux_r(dffts%nnr))
!        !
!        aux_g = 0._DP
!        aux_r = 0._DP
!        !
!        DO ig = 1, npwq  ! perturbation acts only on body
550
!           aux_g(ig) = dvg(ig,ipert) * pot3D%sqvc(ig)
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
!        ENDDO
!        !
!        IF(gamma_only) THEN
!          CALL single_invfft_gamma(dffts,npwq,npwqx,aux_g,aux_r,TRIM(fftdriver))
!        ELSE
!          CALL single_invfft_k(dffts,npwq,npwqx,aux_g,aux_r,TRIM(fftdriver)) ! no igk
!        ENDIF
!        !
!        ! The perturbation is in aux_r
!        !
!        dvpsi=0._DP
!        dpsi =0._DP
!        !
!        IF(gamma_only) THEN
!           !
!           ! double bands @ gamma
!           DO ibnd=1,nbndval-MOD(nbndval,2),2
!              !
!              CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),evc(1,ibnd+1),psic,'Wave')
!              DO ir=1,dffts%nnr
!                 psic(ir) = psic(ir) * REAL(aux_r(ir),KIND=DP)
!              ENDDO
!              CALL double_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,ibnd),dvpsi(1,ibnd+1),'Wave')
!              !
!           ENDDO
576
!           !
577
578
579
580
581
582
583
584
!           ! single band @ gamma
!           IF( MOD(nbndval,2) == 1 ) THEN
!              ibnd=nbndval
!              !
!              CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),psic,'Wave')
!              DO ir=1,dffts%nnr
!                 psic(ir) = CMPLX( REAL(psic(ir),KIND=DP) * REAL(aux_r(ir),KIND=DP), 0._DP, KIND=DP)
!              ENDDO
585
!              CALL single_fwfft_gamma(dffts,npw,npwx,psic,dvpsi(1,ibnd),'Wave')
586
587
588
589
590
591
592
593
594
595
596
597
!              !
!           ENDIF
!           !
!        ELSE
!           !
!           ! only single bands
!           DO ibnd=1,nbndval
!              !
!              CALL single_invfft_k(dffts,npw,npwx,evc(1,ibnd),psic,'Wave',igk_k(1,current_k))
!              DO ir=1,dffts%nnr
!                 psic(ir) = psic(ir) * aux_r(ir)
!              ENDDO
598
!              CALL single_fwfft_k(dffts,npw,npwx,psic,dvpsi(1,ibnd),'Wave',igk_k(1,current_k))
599
600
601
602
603
604
605
606
607
608
!              !
!           ENDDO
!           !
!           IF(npol==2) THEN
!              DO ibnd=1,nbndval
!                 !
!                 CALL single_invfft_k(dffts,npw,npwx,evc(npwx+1,ibnd),psic,'Wave',igk_k(1,current_k))
!                 DO ir=1,dffts%nnr
!                    psic(ir) = psic(ir) * aux_r(ir)
!                 ENDDO
609
!                 CALL single_fwfft_k(dffts,npw,npwx,psic,dvpsi(npwx+1,ibnd),'Wave',igk_k(1,current_k))
610
611
612
613
!                 !
!              ENDDO
!           ENDIF
!           !
614
!        ENDIF
615
616
617
618
619
620
!        !
!        DEALLOCATE(aux_g)
!        DEALLOCATE(aux_r)
!        !
!        CALL apply_alpha_pc_to_m_wfcs(nbndval,nbndval,dvpsi,(-1._DP,0._DP))
!        !
621
!        CALL precondition_m_wfcts( nbndval, dvpsi, dpsi, eprec )
622
623
624
!        !
!        IF( l_dost) THEN
!           !
625
626
627
!           CALL linsolve_sternheimer_m_wfcts (nbndval, nbndval, dvpsi, dpsi, et(1,iks), eprec, tr2, ierr )
!           !
!           IF(ierr/=0) THEN
628
629
630
631
632
633
634
635
!              WRITE(stdout, '(7X,"** WARNING : PERT ",i8," not converged, ierr = ",i8)') ipert,ierr
!           ENDIF
!           !
!        ENDIF
!        !
!        ALLOCATE(aux_r(dffts%nnr))
!        !
!        aux_r=0._DP
636
!        !
637
638
639
640
641
642
643
!        IF(gamma_only) THEN
!           !
!           ! double band @ gamma
!           DO ibnd=1,nbndval
!              !
!              CALL double_invfft_gamma(dffts,npw,npwx,evc(1,ibnd),dpsi(1,ibnd),psic,'Wave')
!              DO ir=1,dffts%nnr
644
!                 prod =  REAL( psic(ir),KIND=DP) * DIMAG( psic(ir))
645
646
647
648
649
650
651
652
653
654
!                 aux_r(ir) = aux_r(ir) + CMPLX( prod, 0.0_DP, KIND=DP)
!              ENDDO
!              !
!           ENDDO
!           !
!        ELSE
!           !
!           ALLOCATE(dpsic(dffts%nnr))
!           !
!           ! only single bands
655
!           DO ibnd=1,nbndval
656
657
658
!              !
!              CALL single_invfft_k(dffts,npw,npwx,evc(1,ibnd),psic,'Wave',igk_k(1,current_k))
!              CALL single_invfft_k(dffts,npw,npwx,dpsi(1,ibnd),dpsic,'Wave',igk_k(1,current_k))
659
!              DO ir=1,dffts%nnr
660
661
662
663
664
665
!                 aux_r(ir) = aux_r(ir) + DCONJG(psic(ir))*dpsic(ir)
!              ENDDO
!              !
!           ENDDO
!           !
!           IF(npol==2) THEN
666
!              DO ibnd=1,nbndval
667
668
669
!                 !
!                 CALL single_invfft_k(dffts,npw,npwx,evc(npwx+1,ibnd),psic,'Wave',igk_k(1,current_k))
!                 CALL single_invfft_k(dffts,npw,npwx,dpsi(npwx+1,ibnd),dpsic,'Wave',igk_k(1,current_k))
670
!                 DO ir=1,dffts%nnr
671
672
673
674
675
676
677
678
679
680
681
682
683
684
!                    aux_r(ir) = aux_r(ir) + DCONJG(psic(ir))*dpsic(ir)
!                 ENDDO
!                 !
!              ENDDO
!           ENDIF
!           !
!           DEALLOCATE(dpsic)
!           !
!        ENDIF
!        !
!        ! The perturbation is in aux_r
!        !
!        ALLOCATE(aux_g(npwqx))
!        IF(gamma_only) THEN
685
!           CALL single_fwfft_gamma(dffts,npwq,npwqx,aux_r,aux_g,TRIM(fftdriver))
686
687
688
689
690
691
692
693
694
695
696
!        ELSE
!           CALL single_fwfft_k(dffts,npwq,npwqx,aux_r,aux_g,TRIM(fftdriver)) ! no igk
!        ENDIF
!        !
!        DO ig=1,npwq ! pert acts only on body
!           dng(ig,ipert) = dng(ig,ipert) + 2._DP * wk(iks) * aux_g(ig) * pot3D%sqvc(ig) / omega
!        ENDDO
!        !
!        DEALLOCATE(aux_g)
!        DEALLOCATE(aux_r)
!        !
697
!        CALL update_bar_type( barra,'dfpt', 1 )
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
!        !
!     ENDDO ! ipert
!     !
!     IF( m==0 ) CALL update_bar_type( barra,'dfpt', 1 )
!     !
!     DEALLOCATE(eprec)
!     DEALLOCATE(dpsi)
!     DEALLOCATE(dvpsi)
!     !
!  ENDDO ! K-POINT and SPIN
!  !
!  IF( gstart==2 ) dng(1,1:m) = CMPLX( 0._DP, 0._DP, KIND=DP )
!  !
!  CALL mp_sum(dng,inter_pool_comm)
!  !
!  CALL mp_barrier( world_comm )
!  !
!  CALL stop_bar_type( barra, 'dfpt' )
!  !
!!  CALL close_buffer(iuwfc,'delete')
!  !IF( nks > 1 ) THEN
!  !   IF ( exst ) THEN
!  !      CLOSE(unit=iuwfc,status='keep')
!  !   ELSE
!  !      CLOSE(unit=iuwfc,status='delete')
!  !   ENDIF
!  !ENDIF
!  !
!END SUBROUTINE