TABLE OF CONTENTS


ABINIT/m_rf2_init [ Modules ]

[ Top ] [ Modules ]

NAME

  m_rf2_init

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group ()
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

TODO

  Can be merged with m_rf2

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 module m_rf2_init
31 
32  use defs_basis
33  use m_errors
34  use m_abicore
35 
36  implicit none
37 
38  private

ABINIT/rf2_init [ Functions ]

[ Top ] [ Functions ]

NAME

 rf2_init

FUNCTION

 Compute terms needed for the 2nd order Sternheimer equation.
 All terms are stored in a rf2_t object.

COPYRIGHT

 Copyright (C) 2015-2018 ABINIT group (LB,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cg(2,mpw*nspinor*mband*nsppol)=planewave coefficients of wavefunctions at k
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
              projected with non-local projectors: cprj=<p_i|Cnk>
  rf2 : the object we want to initialize (see m_rf2.F90 for more information)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eig0_k(mband*nsppol)=GS eigenvalues at k (hartree)
  eig1_k(2*mband*mband*nsppol)=2nd-order eigenvalues at k,q (hartree)
  ffnl1=nonlocal form factors
  ffnl1_test=nonlocal form factors used for tests (i.e. when dtset%nonlinear_info>2)
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  ibg=shift to be applied on the location of data in the array cprj
  icg=shift to be applied on the location of data in the array cg
  idir=direction of the perturbation
  ikpt=number of the k-point
  ipert=type of the perturbation
  isppol=index of current spin component
  mkmem =number of k points trated by this node (GS data).
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw or wfs at k
  nband_k=number of bands at this k point for that spin polarization
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
  rf_hamk_dir2 <type(rf_hamiltonian_type)>= (used only when ipert=natom+11, so q=0)
    same as rf_hamkq, but the direction of the perturbation is different
  occ_k(nband_k)=occupation number for each band (usually 2) for each k.
  rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
  ddk_f<wfk_t>=struct info for DDK file.

OUTPUT

  rf2%RHS_Stern

NOTES

PARENTS

      dfpt_vtowfk

CHILDREN

      cg_zaxpy,dotprod_g,getcprj,pawcprj_alloc,pawcprj_free,pawcprj_get
      rf2_accumulate_bands,rf2_apply_hamiltonian,rf2_getidirs,sqnorm_g
      wfk_read_bks,wrtout,xmpi_allgather,xmpi_barrier

SOURCE

110 subroutine rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,ffnl1,ffnl1_test,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,&
111                      mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f)
112 
113  use defs_basis
114  use defs_datatypes
115  use defs_abitypes
116  use m_xmpi
117  use m_errors
118  use m_wfk
119  use m_hamiltonian
120  use m_cgtools
121  use m_rf2
122  use m_time   , only : timab
123  use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_copy,pawcprj_get,pawcprj_free,pawcprj_output
124  use m_cgprj,   only : getcprj
125 
126 !This section has been created automatically by the script Abilint (TD).
127 !Do not modify the following lines by hand.
128 #undef ABI_FUNC
129 #define ABI_FUNC 'rf2_init'
130 !End of the abilint section
131 
132  implicit none
133 
134 ! *************************************************************************
135 !Arguments -------------------------------
136 !scalars
137 
138  integer,intent(in) :: ibg,icg,idir,ipert,isppol,ikpt
139  integer,intent(in) :: mkmem,mpw,nband_k,nsppol
140  type(datafiles_type),intent(in) :: dtfil
141  type(dataset_type),intent(in) :: dtset
142  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
143  type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq,rf_hamk_dir2
144  type(MPI_type),intent(in) :: mpi_enreg
145 
146 !arrays
147  real(dp),intent(in),target :: cg(2,mpw*gs_hamkq%nspinor*dtset%mband*mkmem*nsppol)
148  real(dp),intent(in) :: eig0_k(dtset%mband)
149  real(dp),intent(inout) :: eig1_k(2*dtset%mband**2) ! Here eig1_k contains 2nd order eigenvalues...
150  real(dp),intent(in) :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
151  real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k)
152  type(pawcprj_type),intent(in) :: cprj(gs_hamkq%natom,gs_hamkq%nspinor*dtset%mband*mkmem*nsppol*gs_hamkq%usecprj)
153  type(rf2_t),intent(inout) :: rf2
154  type(wfk_t),intent(inout) :: ddk_f(4)
155 !
156 !Local variables-------------------------------
157 !scalars
158  integer,parameter :: iorder_cprj=0
159  integer :: choice_cprj,cpopt_cprj,iband,icpgr_loc,idir1,idir2,idir_cprj,ierr
160  integer :: indb,ipert1,ipert2,iproc,jband,kdir1
161  integer :: me,my_nband,natom,ncpgr_loc,nproc_band,debug_mode
162  integer :: size_cprj,size_wf,shift_band1,shift_band2,shift_cprj_band1,shift_cprj_dir1,shift_proc
163  integer :: shift_dir1_lambda,shift_dir2_lambda,shift_dir1,shift_dir1_loc,shift_dir2,shift_jband_lambda
164  logical :: has_cprj_jband,has_dudkprj
165  real(dp) :: doti,dotr,dot2i,dot2r,invocc,tol_final,factor
166  character(len=500) :: msg
167 !arrays
168  integer :: file_index(2)
169  real(dp) :: lambda_ij(2),tsec(2)
170  real(dp),allocatable :: cg_jband(:,:,:),ddk_read(:,:),dudkdk(:,:),dudk_dir2(:,:)
171  real(dp),allocatable :: eig1_read(:),gvnl1(:,:),h_cwave(:,:),s_cwave(:,:),dsusdu_loc(:,:),dsusdu_gather(:,:)
172  real(dp),allocatable,target :: dsusdu(:,:),dudk(:,:),eig1_k_stored(:)
173  real(dp), ABI_CONTIGUOUS pointer :: cwave_dudk(:,:),cwave_i(:,:),cwave_j(:,:),eig1_k_jband(:)
174  real(dp),pointer :: rhs_j(:,:)
175  type(pawcprj_type),target :: cprj_empty(0,0)
176  type(pawcprj_type),allocatable,target :: cprj_jband(:,:),dudkprj(:,:)
177  type(pawcprj_type),pointer :: cprj_dudk(:,:),cprj_j(:,:)
178  type(rf_hamiltonian_type),pointer :: rf_hamk_idir
179 
180 ! *********************************************************************
181 
182  DBG_ENTER("COLL")
183 
184  call timab(514,1,tsec)
185 
186 !my mpi rank :
187  me=mpi_enreg%me_kpt
188 
189  size_wf=gs_hamkq%npw_k*gs_hamkq%nspinor
190  size_cprj=gs_hamkq%nspinor
191  natom = gs_hamkq%natom
192  debug_mode = 0
193  if (dtset%nonlinear_info>2) debug_mode = 1 ! also active a lot of tests
194 
195 !Define some attributes of the rf2 object
196  rf2%nband_k = nband_k
197  rf2%size_wf = size_wf
198  rf2%size_cprj = size_cprj
199 
200  if(ipert<natom+10.or.ipert>natom+11) then
201    write(msg,'(a)') 'ipert must be equal to natom+10 or natom+11 for rf2 calculations.'
202    MSG_BUG(msg)
203  end if
204 
205 !Define perturbations and idirs
206  rf2%iperts(1) = natom+1
207  rf2%iperts(2) = natom+1
208  if (ipert==natom+11)  rf2%iperts(2) = natom+2
209 
210  if (ipert==natom+10.and.idir<=3) then ! One perturbation, one direction
211    rf2%ndir=1
212    rf2%idirs(1)=idir ; rf2%idirs(2)=idir
213  else ! Two perturbations or/and two directions
214    rf2%ndir=2
215    call rf2_getidirs(idir,idir1,idir2)
216    rf2%idirs(1)=idir1
217    rf2%idirs(2)=idir2
218  end if
219 
220 ! **************************************************************************************************
221 ! Get info from ddk files
222 ! **************************************************************************************************
223 
224 !Allocate work spaces
225  ABI_ALLOCATE(eig1_read,(2*nband_k))
226  ABI_ALLOCATE(ddk_read,(2,size_wf))
227  eig1_read(:)=zero
228  ddk_read(:,:)=zero
229 
230 ! "eig1_k_stored" contains dLambda_{nm}/dpert every bands n and m and ndir (=1 or 2) directions
231 ! pert = k_dir (wavevector) or E_dir (electric field)
232  ABI_ALLOCATE(eig1_k_stored,(2*rf2%ndir*nband_k**2))
233  eig1_k_stored=zero
234 
235 ! "dudk" contains du/dpert1 for every bands and ndir (=1 or 2) directions
236  ABI_STAT_ALLOCATE(dudk,(2,rf2%ndir*nband_k*size_wf), ierr)
237  ABI_CHECK(ierr==0, "out of memory in m_rf2 : dudk")
238  dudk=zero
239  has_dudkprj=.false.
240  if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
241    ABI_DATATYPE_ALLOCATE(dudkprj,(natom,rf2%ndir*nband_k*size_cprj))
242    ncpgr_loc=1;if(ipert==natom+10.or.ipert==natom+11) ncpgr_loc=3
243    call pawcprj_alloc(dudkprj,ncpgr_loc,gs_hamkq%dimcprj)
244    choice_cprj=5 ; cpopt_cprj=0
245    has_dudkprj=.true.
246  else
247    ABI_DATATYPE_ALLOCATE(dudkprj,(natom,0))
248  end if
249 
250  if (debug_mode/=0) then
251    write(msg,'(4(a,i2))') 'RF2_INIT : ipert-natom = ',ipert-natom,' , idir = ',idir,&
252    ' , ikpt = ',ikpt,' , isppol = ',isppol
253    call wrtout(std_out,msg,'COLL')
254  end if
255 
256  file_index(1)=1 ! dir1
257  file_index(2)=2 ! dir2
258  if (ipert==natom+11) then ! see dfpt_looppert.F90
259    file_index(1)=3 ! dir1
260    file_index(2)=2 ! dir2
261  end if
262 
263  do kdir1=1,rf2%ndir
264    idir1=rf2%idirs(kdir1)
265    ipert1=rf2%iperts(kdir1)
266    do iband=1,nband_k
267      call wfk_read_bks(ddk_f(file_index(kdir1)),iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read)
268 !    Copy ddk_read in "dudk"
269      shift_band1=(iband-1)*size_wf
270      shift_dir1=(kdir1-1)*nband_k*size_wf
271      dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=ddk_read(:,:)
272 !    Copy eig1_read in "eig1_k_stored"
273      shift_band1=(iband-1)*2*nband_k
274      shift_dir1_lambda=2*(kdir1-1)*nband_k**2
275      eig1_k_stored(1+shift_band1+shift_dir1_lambda:2*nband_k+shift_band1+shift_dir1_lambda)=eig1_read(:)
276 !    Get this dudk projected on NL projectors
277      if (has_dudkprj.and.mpi_enreg%proc_distrb(ikpt,iband,isppol)==me) then
278        shift_cprj_band1=(iband-1)*size_cprj
279        shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj
280        cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1: &
281 &       size_cprj+shift_cprj_band1+shift_cprj_dir1)
282        idir_cprj=0;if (dudkprj(1,1)%ncpgr/=3) idir_cprj=idir1
283        call getcprj(choice_cprj,cpopt_cprj,ddk_read,cprj_dudk,gs_hamkq%ffnl_k,idir_cprj,&
284 &       gs_hamkq%indlmn,gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kpg_k,gs_hamkq%kpt_k,&
285 &       gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,&
286 &       gs_hamkq%nloalg,gs_hamkq%npw_k,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkxred,&
287 &       gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm)
288      end if
289    end do
290  end do
291 
292  ABI_ALLOCATE(dudkdk,(2,0))
293  ABI_ALLOCATE(dudk_dir2,(2,0))
294 
295 !Get dudkdk for ipert==natom+11
296  if(ipert==natom+11) then
297    ABI_DEALLOCATE(dudkdk)
298    ABI_ALLOCATE(dudkdk,(2,nband_k*size_wf))
299    if (idir>3) then
300      ABI_DEALLOCATE(dudk_dir2)
301      ABI_ALLOCATE(dudk_dir2,(2,nband_k*size_wf))
302    end if
303    do iband=1,nband_k
304      call wfk_read_bks(ddk_f(1),iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read)
305      shift_band1=(iband-1)*size_wf
306      dudkdk(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:)
307 !    Check that < u^(0) | u^(2) > = - Re[< u^(1) | u^(1) >]
308      if (debug_mode/=0 .and. gs_hamkq%usepaw==0) then
309 !      Compute < u^(0) | u^(2) >
310        do jband=1,nband_k
311          cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
312          call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,ddk_read,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
313          if (idir<=3 .and. iband==jband .and. abs(occ_k(iband))>tol8) then
314 !          Compute < u^(1) | u^(1) > = Re[< u^(1) | u^(1) >]
315            cwave_dudk => dudk(:,1+shift_band1:size_wf+shift_band1)
316            call sqnorm_g(dot2r,gs_hamkq%istwf_k,size_wf,cwave_dudk,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
317            dotr = dotr + dot2r
318            dotr = sqrt(dotr**2+doti**2)
319            if (dotr > tol7) then
320              write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST dudkdk iband = ',iband,&
321              ' : NOT PASSED. | < u^(0) | u^(2) > + Re[< u^(1) | u^(1) >] | = ',dotr
322              call wrtout(std_out,msg)
323              call wrtout(ab_out,msg)
324 !           else
325 !             write(msg,'(a,i2,a)') 'RF2 TEST dudkdk iband = ',iband,' : OK.'
326 !             call wrtout(std_out,msg)
327            end if
328          end if ! idir<=3
329        end do ! jband
330      end if ! debug_mode
331 !    Read ddk for idir2
332      if (idir>3) then
333        call wfk_read_bks(ddk_f(4),iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read)
334        dudk_dir2(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:)
335      end if
336    end do !iband
337  end if ! ipert=natom+11
338 
339  ABI_DEALLOCATE(ddk_read)
340  ABI_DEALLOCATE(eig1_read)
341 
342 ! **************************************************************************************************
343 ! COMPUTATION OF "dsusdu", A PART OF "A_mn" AND A PART OF "Lambda_mn" (see defs in m_rf2)
344 ! **************************************************************************************************
345 
346 !Allocate work spaces for one band
347  ABI_ALLOCATE(h_cwave,(2,size_wf))
348  ABI_ALLOCATE(s_cwave,(2,size_wf))
349  ABI_ALLOCATE(gvnl1,(2,size_wf))
350  h_cwave(:,:) = zero
351  s_cwave(:,:) = zero
352  gvnl1(:,:) = zero
353 
354 ! "dsusdu" contains dS/dpert_dir |u_band> + S|du_band/dpert1> for every bands and ndir (=1 or 2) directions
355  ABI_STAT_ALLOCATE(dsusdu,(2,rf2%ndir*nband_k*size_wf), ierr)
356  ABI_CHECK(ierr==0, "out of memory in rf2_init : dsusdu")
357  dsusdu=zero
358 
359  ABI_ALLOCATE(rf2%amn,(2,nband_k**2))
360  rf2%amn=zero
361 
362  ABI_ALLOCATE(rf2%lambda_mn,(2,nband_k**2))
363  rf2%lambda_mn(:,:)=zero
364 
365 !Allocate work spaces when debug_mode is activated
366  has_cprj_jband=.false.
367  if (debug_mode/=0) then ! Only for test purposes
368    ABI_ALLOCATE(cg_jband,(2,size_wf*nband_k,2))
369    cg_jband(:,:,1) = cg(:,1+icg:size_wf*nband_k+icg)
370    if (ipert==natom+11) then ! Note the multiplication by "i"
371      if (idir<=3) then
372        cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k) ! for dir1
373        cg_jband(2,:,2) =  dudk(1,1:size_wf*nband_k) ! for dir1
374      else
375        cg_jband(1,:,2) = -dudk_dir2(2,1:size_wf*nband_k) ! for dir2
376        cg_jband(2,:,2) =  dudk_dir2(1,1:size_wf*nband_k) ! for dir2
377      end if
378    end if
379    if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
380      ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,size_cprj*nband_k))
381      has_cprj_jband=.true.
382    else
383      ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,0))
384    end if
385  else
386    ABI_ALLOCATE(cg_jband,(2,0,2))
387    ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,0))
388  end if
389 
390  factor=one
391  if(ipert==natom+10 .and. idir<=3) factor=two ! in order to not compute same terms twice
392 
393  do kdir1=1,rf2%ndir
394 !  First iteration (kdir1=1) :
395 !  pert1 = rf2%iperts(1) along rf2%idirs(1)
396 !  pert2 = rf2%iperts(2) along rf2%idirs(2)
397 !  Second iteration (kdir1=2) :
398 !  pert1 = rf2%iperts(2) along rf2%idirs(2)
399 !  pert2 = rf2%iperts(1) along rf2%idirs(1)
400    idir1=rf2%idirs(kdir1)
401    ipert1=rf2%iperts(kdir1)
402    shift_dir1=(kdir1-1)*nband_k*size_wf
403    shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj
404    shift_dir1_lambda=(kdir1-1)*2*nband_k**2
405    if(ipert==natom+10 .and. idir<=3) then
406      shift_dir2=0
407      idir2 = idir1
408      ipert2 = ipert1
409      rf_hamk_idir => rf_hamkq
410    else
411      shift_dir2=(2-kdir1)*nband_k*size_wf
412      idir2 = rf2%idirs(3-kdir1)
413      ipert2 = rf2%iperts(3-kdir1)
414      if (kdir1==1) rf_hamk_idir => rf_hamkq
415      if (kdir1==2) rf_hamk_idir => rf_hamk_dir2
416    end if
417 
418 !  Load projected WF according to ipert1 and idir1
419    cprj_j => cprj_empty ; cprj_dudk => cprj_empty
420    if (has_cprj_jband) then
421      call pawcprj_free(cprj_jband)
422      ncpgr_loc= 3;if(ipert1==natom+1.or.ipert1==natom+2) ncpgr_loc=1
423      icpgr_loc=-1;if(ipert1==natom+1.or.ipert1==natom+2) icpgr_loc=idir1
424      call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj)
425      call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,&
426 &     isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,&
427 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc)
428    end if
429 
430 !  LOOP OVER BANDS
431    do jband=1,nband_k ! = band n
432 
433 !    Skip bands not treated by current proc
434      if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
435 
436      shift_band1=(jband-1)*size_wf
437      shift_cprj_band1=(jband-1)*size_cprj
438      shift_jband_lambda=(jband-1)*2*nband_k
439 
440      if (abs(occ_k(jband))>tol8) then
441 
442 !      Extract first order wavefunction and eigenvalues for jband
443        eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir1_lambda:2*nband_k+shift_jband_lambda+shift_dir1_lambda)
444        cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
445        if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1)
446 
447 !      Compute H^(0) | du/dpert1 > (in h_cwave) and S^(0) | du/dpert1 > (in s_cwave)
448        call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_dudk,cprj_dudk,h_cwave,s_cwave,&
449 &       eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,0,0,ikpt,isppol,mkmem,&
450 &       mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf)
451 
452        if (gs_hamkq%usepaw==0) s_cwave(:,:)=cwave_dudk(:,:) ! Store | du/dpert1 > in s_cwave
453 
454 !      Copy infos in dsusdu
455        dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)&
456        +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
457 
458        if (debug_mode/=0) then
459          write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 1 kdir1 = ',kdir1,' jband = ',jband
460          call wrtout(std_out,msg)
461        end if
462 
463 !      For every occupied iband, we compute :
464 !      < du/dpert2(iband) | H^(0) | du/dpert1(jband) > and add it to lambda_mn
465 !      < du/dpert2(iband) | S^(0) | du/dpert1(jband) > and add it to amn
466        do iband=1,rf2%nband_k  ! = band m
467          if (abs(occ_k(iband))>tol8) then
468            shift_band2=(iband-1)*size_wf
469            cwave_dudk => dudk(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2)
470            call rf2_accumulate_bands(rf2,1,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
471            jband,debug_mode,cwave_dudk,h_cwave,s_cwave)
472          end if
473        end do
474 
475 !      Extract GS wavefunction for jband
476        cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
477        if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1)
478 
479        if (ipert1==natom+2) then
480 !        Extract ddk and multiply by i :
481          if(idir<=3) then ! in this case : idir1=idir2
482            gvnl1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1)
483            gvnl1(2,:) =  dudk(1,1+shift_band1:size_wf+shift_band1)
484          else
485            gvnl1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1)
486            gvnl1(2,:) =  dudk_dir2(1,1+shift_band1:size_wf+shift_band1)
487          end if
488        end if
489 
490 !      Compute dH/dpert1 | u^(0) > (in h_cwave) and dS/dpert1 | u^(0) > (in s_cwave)
491        call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_j,cprj_j,h_cwave,s_cwave,&
492 &       eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,idir1,ipert1,ikpt,isppol,&
493 &       mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf)
494 
495 !      Copy infos in dsusdu
496        if (gs_hamkq%usepaw==1) then
497          dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)&
498          +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
499        end if
500 
501        if (debug_mode/=0) then
502          write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 2 kdir1 = ',kdir1,' jband = ',jband
503          call wrtout(std_out,msg)
504        end if
505 
506 !      For every occupied iband, we compute :
507 !      < du/dpert2(iband) | dH/dpert1 | u^(0)(jband) > and add it to lambda_mn
508 !      < du/dpert2(iband) | dS/dpert1 | u^(0)(jband) > and add it to amn
509        do iband=1,rf2%nband_k  ! = band m
510          if (abs(occ_k(iband))>tol8) then
511            shift_band2=(iband-1)*size_wf
512            cwave_dudk => dudk(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2)
513            call rf2_accumulate_bands(rf2,2,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
514            jband,debug_mode,cwave_dudk,h_cwave,s_cwave)
515          end if
516        end do
517 
518      end if ! empty band test
519    end do ! jband
520  end do ! idir1
521 
522 ! Allgather dsusdu
523  nproc_band = xmpi_comm_size(mpi_enreg%comm_band)
524  if (nproc_band>1) then
525 
526    my_nband = nband_k/nproc_band;if (mod(nband_k,nproc_band)/=0) my_nband=my_nband+1
527    ABI_ALLOCATE(dsusdu_loc,(2,size_wf*my_nband*rf2%ndir))
528    ABI_ALLOCATE(dsusdu_gather,(2,size_wf*my_nband*rf2%ndir*nproc_band))
529    dsusdu_loc(:,:) = zero
530    dsusdu_gather(:,:) = zero
531 
532    do kdir1=1,rf2%ndir
533      indb = 1
534      shift_dir1=(kdir1-1)*size_wf*nband_k
535      shift_dir1_loc=(kdir1-1)*size_wf*my_nband
536      do jband=1,nband_k
537 !      Skip bands not treated by current proc
538        if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
539 
540        shift_band1=(jband-1)*size_wf
541        dsusdu_loc(:,indb+shift_dir1_loc:indb-1+size_wf+shift_dir1_loc) = &
542        dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
543        indb = indb + size_wf
544      end do
545    end do
546 
547    call xmpi_allgather(dsusdu_loc,2*size_wf*my_nband*rf2%ndir,dsusdu_gather,mpi_enreg%comm_band,ierr)
548 
549    do kdir1=1,rf2%ndir
550      shift_dir1=(kdir1-1)*size_wf*nband_k
551      shift_dir1_loc=(kdir1-1)*size_wf*my_nband
552      do iproc=1,nproc_band
553        shift_proc = (iproc-1)*size_wf*my_nband*rf2%ndir
554        indb = 1
555        do jband=1,my_nband
556          iband = jband+(iproc-1)*my_nband
557          if(iband<=nband_k) then
558            shift_band1=(iband-1)*size_wf
559            dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) = &
560            dsusdu_gather(:,indb+shift_dir1_loc+shift_proc:indb-1+size_wf+shift_dir1_loc+shift_proc)
561          end if
562          indb = indb + size_wf
563        end do
564      end do
565    end do
566    ABI_DEALLOCATE(dsusdu_loc)
567    ABI_DEALLOCATE(dsusdu_gather)
568 
569  end if
570 
571 ! **************************************************************************************************
572 ! COMPUTATION OF "RHS_Stern", THE LAST PART OF "A_mn" AND A PART OF "Lambda_mn"
573 ! **************************************************************************************************
574 
575  ABI_STAT_ALLOCATE(rf2%RHS_Stern,(2,nband_k*size_wf), ierr)
576  ABI_CHECK(ierr==0, "out of memory in m_rf2 : RHS_Stern")
577  rf2%RHS_Stern(:,:)=zero
578 
579 !Computation of terms containing H^(2)
580  if (ipert/=natom+11 .or. gs_hamkq%usepaw==1) then ! Otherwise H^(2) = 0
581 
582 !Load projected WF according to ipert and idir
583    cprj_j => cprj_empty
584    if (has_cprj_jband) then
585      call pawcprj_free(cprj_jband)
586      ncpgr_loc= 3;if(ipert==natom+1.or.ipert==natom+2) ncpgr_loc=1
587      icpgr_loc=-1;if(ipert==natom+1.or.ipert==natom+2) icpgr_loc=idir
588      call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj)
589      call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,&
590 &     isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,&
591 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc)
592    end if
593 
594    if (ipert==natom+10) then
595      rf_hamk_idir => rf_hamkq !     all info are in rf_hamkq
596    else if (ipert==natom+11) then
597      rf_hamk_idir => rf_hamk_dir2 ! all info are in rf_hamk_dir2
598    end if
599 
600    do jband=1,nband_k
601 
602 !    Skip bands not treated by current proc
603      if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
604 
605      if (abs(occ_k(jband))>tol8) then
606        shift_band1=(jband-1)*size_wf
607        shift_cprj_band1=(jband-1)*size_cprj
608        shift_jband_lambda=(jband-1)*2*nband_k
609 
610 !      Extract GS wavefunction
611        cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
612        if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1)
613 
614 !      Not used here, but a null pointer is not allowed... (called by rf2_apply_hamiltonian but not used)
615        eig1_k_jband => eig1_k_stored(1+shift_jband_lambda:2*nband_k+shift_jband_lambda)
616 
617        if (ipert==natom+11) then
618 !        Extract ddk and multiply by i :
619          if(idir<=3) then ! in this case : idir1=idir2
620            gvnl1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1)
621            gvnl1(2,:) =  dudk(1,1+shift_band1:size_wf+shift_band1)
622          else
623            gvnl1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1)
624            gvnl1(2,:) =  dudk_dir2(1,1+shift_band1:size_wf+shift_band1)
625          end if
626        end if
627 
628 !      Compute  : d^2H/(dpert1 dpert2)|u^(0)>  (in h_cwave)
629 !      and      : d^2S/(dpert1 dpert2)|u^(0)>  (in s_cwave)
630        call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_j,cprj_j,h_cwave,s_cwave,&
631 &       eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,idir,ipert,ikpt,isppol,&
632 &       mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf,&
633 &       ffnl1=ffnl1,ffnl1_test=ffnl1_test)
634 
635        if (debug_mode/=0) then
636          write(msg,'(a,i2)') 'RF2 TEST before accumulate_bands choice = 3 jband = ',jband
637          call wrtout(std_out,msg)
638        end if
639 
640 !      For every occupied iband, we compute :
641 !      < u^(0)(iband) | d^2H/(dpert1 dpert2) | u^(0)(jband) > and add it to lambda_mn
642 !      < u^(0)(iband) | d^2S/(dpert1 dpert2) | u^(0)(jband) > and add it to amn
643        do iband=1,rf2%nband_k  ! = band m
644          if (abs(occ_k(iband))>tol8) then
645            shift_band2=(iband-1)*size_wf
646            cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
647            if(ipert == natom+10) then
648              ipert1 = natom+1
649              ipert2 = natom+1
650            else
651              ipert1 = natom+1
652              ipert2 = natom+2
653            end if
654            call rf2_getidirs(idir,idir1,idir2)
655            call rf2_accumulate_bands(rf2,3,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
656            jband,debug_mode,cwave_i,h_cwave,s_cwave)
657          end if
658        end do
659 
660 !      Add d^2H/(dk_dir1 dk_dir2)|u^(0)> to RHS_Stern :
661        if (gs_hamkq%usepaw==1) h_cwave(:,:)=h_cwave(:,:)-eig0_k(jband)*s_cwave(:,:) ! if PAW : we add H^(2)-eps^(0) S^(2)
662        rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
663        call cg_zaxpy(size_wf,(/one,zero/),h_cwave,rhs_j)
664 
665      end if ! empty band test
666    end do ! jband
667  end if ! H^(2) exists
668 
669 !Computation of terms containing H^(1)
670  do kdir1=1,rf2%ndir
671 !  First iteration (kdir1=1) :
672 !  pert1 = rf2%iperts(1) along rf2%idirs(1)
673 !  pert2 = rf2%iperts(2) along rf2%idirs(2)
674 !  Second iteration (kdir1=2) :
675 !  pert1 = rf2%iperts(2) along rf2%idirs(2)
676 !  pert2 = rf2%iperts(1) along rf2%idirs(1)
677    shift_dir1=(kdir1-1)*nband_k*size_wf
678    shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj
679    shift_dir1_lambda=(kdir1-1)*2*nband_k**2
680    idir1=rf2%idirs(kdir1)
681    ipert1=rf2%iperts(kdir1)
682    if(ipert==natom+10 .and. idir<=3) then
683      idir2=idir1
684      ipert2=ipert1
685      shift_dir2=0
686      shift_dir2_lambda=0
687      rf_hamk_idir => rf_hamkq
688    else
689      idir2=rf2%idirs(2-kdir1+1)
690      ipert2=rf2%iperts(2-kdir1+1)
691      shift_dir2=(2-kdir1)*nband_k*size_wf
692      shift_dir2_lambda=(2-kdir1)*2*nband_k**2
693      if (kdir1==1) rf_hamk_idir => rf_hamk_dir2 ! dir2
694      if (kdir1==2) rf_hamk_idir => rf_hamkq ! dir1
695    end if
696 
697 !  Load projected WF according to ipert2 and idir2
698    cprj_j => cprj_empty ;  ; cprj_dudk => cprj_empty
699    if (has_cprj_jband) then
700      call pawcprj_free(cprj_jband)
701      ncpgr_loc= 3;if(ipert2==natom+1.or.ipert2==natom+2) ncpgr_loc=1
702      icpgr_loc=-1;if(ipert2==natom+1.or.ipert2==natom+2) icpgr_loc=idir2
703      call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj)
704      call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,&
705 &     isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,&
706 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc)
707    end if
708 
709    do jband=1,nband_k
710 
711 !    Skip bands not treated by current proc
712      if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
713 
714      if (abs(occ_k(jband))>tol8) then
715        shift_band1=(jband-1)*size_wf
716        shift_cprj_band1=(jband-1)*size_cprj
717        shift_jband_lambda=(jband-1)*2*nband_k
718 
719 !      Extract first order wavefunction | du/dpert1 > and eigenvalues
720        eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir2_lambda:2*nband_k+shift_jband_lambda+shift_dir2_lambda)
721        cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
722        if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1)
723 
724        if (ipert2==natom+2) then
725 !        Extract dkdk and multiply by i :
726          gvnl1(1,:) = -dudkdk(2,1+shift_band1:size_wf+shift_band1)
727          gvnl1(2,:) =  dudkdk(1,1+shift_band1:size_wf+shift_band1)
728        end if
729 
730 !      Compute dH/dpert2 | du/dpert1 > (in h_cwave) and dS/dpert2 | du/dpert1 > (in s_cwave)
731        call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_dudk,cprj_dudk,h_cwave,s_cwave,&
732 &       eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,idir2,ipert2,ikpt,isppol,&
733 &       mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf)
734 
735        if (debug_mode/=0) then
736          write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 4 kdir1 = ',kdir1,' jband = ',jband
737          call wrtout(std_out,msg)
738        end if
739 
740 !      For every occupied iband, we compute :
741 !      < u^(0)(iband) | dH/dpert2 | du/dpert1(jband) > and add it to lambda_mn
742 !      < u^(0)(iband) | dS/dpert2 | du/dpert1(jband) > and add it to amn
743        do iband=1,rf2%nband_k  ! = band m
744          if (abs(occ_k(iband))>tol8) then
745            shift_band2=(iband-1)*size_wf
746            cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
747            call rf2_accumulate_bands(rf2,4,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
748            jband,debug_mode,cwave_i,h_cwave,s_cwave)
749          end if
750        end do
751 
752 !      Add dH/dpert2 | du/dpert1 > to RHS_Stern :
753        if (gs_hamkq%usepaw==1) h_cwave(:,:)=h_cwave(:,:)-eig0_k(jband)*s_cwave(:,:) ! if PAW : we add H^(1)-eps^(0) S^(1)
754        rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
755        call cg_zaxpy(size_wf,(/factor*one,zero/),h_cwave,rhs_j)
756 
757 !      Compute : -factor * sum_iband ( dLambda/dpert1_{iband,jband} * dsusdu_{iband} )
758        do iband=1,nband_k
759          if (abs(occ_k(iband))>tol8) then ! if empty band, nothing to do
760 
761 !          Extract lambda_ij(iband,jband) for dir1
762            lambda_ij(1)=eig1_k_stored(2*iband-1+shift_jband_lambda+shift_dir1_lambda)
763            lambda_ij(2)=eig1_k_stored(2*iband  +shift_jband_lambda+shift_dir1_lambda)
764 
765 !          Extract dsusdu for iband and pert2 (in cwave_i)
766            shift_band2=(iband-1)*size_wf
767            cwave_i => dsusdu(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2)
768 
769 !          Compute Lambda_{iband,jband} * dsusdu_{iband} and add it to RHS_Stern
770            rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
771            call cg_zaxpy(size_wf,(/-factor*lambda_ij(1),-factor*lambda_ij(2)/),cwave_i,rhs_j) !do not forget the minus sign!
772 
773          end if ! empty iband test
774        end do ! iband
775 
776      end if ! empty jband test
777    end do ! jband
778 
779  end do ! kdir1
780 
781  ABI_DEALLOCATE(gvnl1)
782  ABI_DEALLOCATE(h_cwave)
783  ABI_DEALLOCATE(s_cwave)
784  ABI_DEALLOCATE(cg_jband)
785  ABI_DEALLOCATE(dudk)
786  ABI_DEALLOCATE(dudkdk)
787  ABI_DEALLOCATE(dudk_dir2)
788  ABI_DEALLOCATE(dsusdu)
789  ABI_DEALLOCATE(eig1_k_stored)
790  if (has_cprj_jband) call pawcprj_free(cprj_jband)
791  ABI_DATATYPE_DEALLOCATE(cprj_jband)
792  if (has_dudkprj) call pawcprj_free(dudkprj)
793  ABI_DATATYPE_DEALLOCATE(dudkprj)
794 
795 ! Compute the part of 2nd order wavefunction that belongs to the space of empty bands
796  do jband=1,nband_k
797 
798 !  Skip bands not treated by current proc
799    if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
800 
801    shift_band1=(jband-1)*size_wf
802    if (abs(occ_k(jband))>tol8) then
803      invocc = one/occ_k(jband)
804      rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
805      do iband=1,nband_k
806        if (iband /= jband) then
807          if (debug_mode/=0) then
808            if (abs(occ_k(iband) - occ_k(jband)) > tol12 .and. occ_k(iband) > tol8) then
809              write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband
810              call wrtout(std_out,msg)
811              call wrtout(ab_out,msg)
812              write(msg,'(a)') 'ERROR : occ_k(iband) /= occ_k(jband) (and both are >0)'
813              call wrtout(std_out,msg)
814              call wrtout(ab_out,msg)
815            end if
816            if (abs(eig0_k(iband) - eig0_k(jband)) < tol8 ) then
817              write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband
818              call wrtout(std_out,msg)
819              write(msg,'(a,es22.13e3)') 'WARNING : DEGENERATE BANDS  Eig(jband) = Eig(jband) = ',eig0_k(jband)
820              call wrtout(std_out,msg)
821            end if
822            if ( (eig0_k(iband) - eig0_k(jband) < -tol12) .and. (jband < iband) ) then
823              write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband
824              call wrtout(std_out,msg)
825              call wrtout(ab_out,msg)
826              write(msg,'(a)') 'ERROR : Eig(jband) < Eig(iband) with jband < iband'
827              call wrtout(std_out,msg)
828              call wrtout(ab_out,msg)
829              write(msg,'(a,es22.13e3)') 'Eig(jband) = ',eig0_k(jband)
830              call wrtout(std_out,msg)
831              call wrtout(ab_out,msg)
832              write(msg,'(a,es22.13e3)') 'Eig(iband) = ',eig0_k(iband)
833              call wrtout(std_out,msg)
834              call wrtout(ab_out,msg)
835            end if
836          end if ! end tests
837          if ( abs(occ_k(iband))<tol8 ) then ! for empty bands only
838            shift_band2=(iband-1)*size_wf
839            cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
840            call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_i,rhs_j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
841 !          Store it in a_mn
842 !          /!\ There is a factor "-2" to simplify the use of amn in the following.
843 !          /!\ Occupied and empty bands will be treated in a same way.
844            rf2%amn(:,iband+(jband-1)*nband_k)=-two*rocceig(iband,jband)*invocc*(/dotr,doti/)&
845            +rf2%amn(:,iband+(jband-1)*nband_k)
846          end if ! empty band test
847        end if ! iband \= jband
848      end do ! iband
849    end if  ! empty band test
850  end do ! jband
851 
852 ! **************************************************************************************************
853 !  COMPUTATION OF "dcwavef" AND "Lambda_mn" FROM "A_mn"
854 ! **************************************************************************************************
855 
856  ABI_STAT_ALLOCATE(rf2%dcwavef,(2,nband_k*size_wf), ierr)
857  ABI_CHECK(ierr==0, "out of memory in m_rf2 : dcwavef")
858  rf2%dcwavef=zero
859 
860  do jband=1,nband_k
861 
862 !  Skip bands not treated by current proc
863    if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
864 
865    shift_band1=(jband-1)*size_wf
866    if (abs(occ_k(jband))>tol8) then
867      do iband=1,nband_k
868        shift_band2=(iband-1)*size_wf
869 
870 !      Extract GS wavefunction for iband
871        cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
872 
873        call cg_zaxpy(size_wf,-half*rf2%amn(:,iband+(jband-1)*nband_k), &
874 &       cwave_i,rf2%dcwavef(:,1+shift_band1))
875 
876        if (abs(occ_k(iband))>tol8 .and. abs(occ_k(jband))>tol8) then
877          rf2%lambda_mn(:,iband+(jband-1)*nband_k) = rf2%lambda_mn(:,iband+(jband-1)*nband_k) &
878          -half*(eig0_k(iband)+eig0_k(jband))*rf2%amn(:,iband+(jband-1)*nband_k)
879 
880          eig1_k(2*iband-1+(jband-1)*2*nband_k) = rf2%lambda_mn(1,iband+(jband-1)*nband_k)
881          eig1_k(2*iband  +(jband-1)*2*nband_k) = rf2%lambda_mn(2,iband+(jband-1)*nband_k)
882 
883        end if ! empty band test
884      end do ! iband
885    end if ! empty band test
886  end do ! jband
887 
888 ! For the following, "rf2%lambda_mn" and "rf2%RHS_Stern" must be computed for every bands
889  call xmpi_barrier(mpi_enreg%comm_band)
890 
891 ! **************************************************************************************************
892 !  FINAL TEST
893 ! **************************************************************************************************
894 
895  tol_final = tol6
896  if (debug_mode/=0) then
897    do jband=1,nband_k
898 
899 !    Skip bands not treated by current proc
900      if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
901 
902      if (abs(occ_k(jband))>tol8) then
903 !       write(msg,'(3(a,i2))') 'RF2 TEST FINAL : ipert=',ipert-natom,' idir=',idir,' jband=',jband
904 !       call wrtout(std_out,msg)
905        shift_band1=(jband-1)*size_wf
906        rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
907        cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
908        call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,rhs_j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
909        dot2r = dotr - rf2%lambda_mn(1,jband+(jband-1)*nband_k)
910        dot2i = doti - rf2%lambda_mn(2,jband+(jband-1)*nband_k)
911        dot2r = sqrt(dot2r**2+dot2i**2)
912        if (dot2r > tol_final) then
913          write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST FINAL iband = ',jband,' : NOT PASSED dotr = ',dotr
914          call wrtout(std_out,msg)
915          call wrtout(ab_out,msg)
916          write(msg,'(2(a,es22.13E3))') ' < cwave_j | rhs_j > =',dotr,',',doti
917          call wrtout(std_out,msg)
918          call wrtout(ab_out,msg)
919          write(msg,'(2(a,es22.13E3))') '           lambda_jj =',&
920 &         rf2%lambda_mn(1,jband+(jband-1)*nband_k),',',rf2%lambda_mn(2,jband+(jband-1)*nband_k)
921          call wrtout(std_out,msg)
922          call wrtout(ab_out,msg)
923        else
924          if (dot2r<tol9) dot2r = zero ! in order to hide the numerical noise
925          write(msg,'(a,i2,a,es22.13E3,a,es7.1E2)') &
926          'RF2 TEST FINAL iband = ',jband,' : OK. |test| = ',dot2r,' < ',tol_final
927          call wrtout(std_out,msg)
928        end if
929      end if
930    end do
931  end if
932 
933 ! **************************************************************************************************
934 !  JOB FINISHED
935 ! **************************************************************************************************
936 
937 ! Deallocations of arrays
938  if (debug_mode==0) then
939    ABI_DEALLOCATE(rf2%amn)
940  end if
941 
942  call timab(514,2,tsec)
943 
944  DBG_EXIT("COLL")
945 
946 end subroutine rf2_init