!{\src2tex{textfont=tt}}
!!****f* ABINIT/loper3
!! NAME
!! loper3
!!
!! FUNCTION
!! Loop over perturbations
!!
!! COPYRIGHT
!! Copyright (C) 1999-2014 ABINIT group (XG, DRH, MB, XW, 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
!!  atindx(natom)=index table for atoms (see gstate.f)
!!  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
!!  codvsn=code version
!!  cpus=cpu time limit in seconds
!!  dimcprj(natom*usepaw)=array of dimensions of arrays cprj, cprjq (ordered by atom-type)
!!  dim_eigbrd=1 if eigbrd is to be computed
!!  dim_eig2nkq=1 if eig2nkq is to be computed
!!  doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy
!!  dtfil <type(datafiles_type)>=variables related to files
!!  dtset <type(dataset_type)>=all input variables for this dataset
!!  dyew(2,3,natom,3,natom)=Ewald part of the dynamical matrix
!!  dyfrlo(3,3,natom)=frozen wavefunctions local part of the dynamical matrix
!!  dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wavefunctions non-local part of the dynamical matrix
!!  dyfrx1(2,3,natom,3,natom)=frozen wf nonlin. xc core corr.(2) part of the dynamical matrix
!!  dyfrx2(3,3,natom)=frozen wf nonlin. xc core corr.(2) part of the dynamical matrix
!!  dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
!!  dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise
!!  eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eigbrd)=boradening factors for the electronic eigenvalues
!!  eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)=second derivatives of the electronic eigenvalues
!!  eltcore(6,6)=core contribution to the elastic tensor
!!  elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor
!!  eltfrhar(6,6)=Hartree contribution to the elastic tensor
!!  eltfrkin(6,6)=kinetic contribution to the elastic tensor
!!  eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor
!!  eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor
!!  eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor
!!  fermie=fermi energy (Hartree)
!!  gsqcut_eff=Fourier cutoff on G^2 for "large sphere" of radius double
!!    that of the basis sphere--appropriate for charge density rho(G),
!!    Hartree potential, and pseudopotentials, corresponding to ecut_eff
!!  iexit=index of "exit" on first line of file (0 if not found)
!!  indsym(4,nsym,natom)=indirect indexing array for atom labels
!!  kxc(nfftf,nkxc)=exchange and correlation kernel (see rhohxc.f)
!!  mkmem =Number of k points treated by this node (GS data)
!!  mkqmem=Number of k+q points treated by this node (GS data)
!!  mk1mem=Number of k points treated by this node (RF data)
!!  mpert=maximum number of ipert
!!  mpi_enreg=informations about MPI parallelization
!!  my_natom=number of atoms treated by current processor
!!  nattyp(ntypat)= # atoms of each type.
!!  nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
!!  nkpt=number of k points
!!  nkxc=second dimension of the kxc array
!!  nspden=number of spin-density components
!!  nsym=number of symmetry elements in space group
!!  occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
!!  paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS
!!  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
!!  pawang <type(pawang_type)>=paw angular mesh and related data
!!  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
!!  pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
!!  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
!!  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
!!  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
!!  pertsy(3,mpert)=set of perturbations that form a basis for all other perturbations
!!  prtbbb=if 1, bbb decomposition, also dimension d2bbb
!!  psps <type(pseudopotential_type)>=variables related to pseudopotentials
!!  rfpert(mpert)=array defining the type of perturbations that have to be computed
!!                1   ->   element has to be computed explicitely
!!               -1   ->   use symmetry operations to obtain the corresponding element
!!  rhog(2,nfftf)=array for Fourier transform of GS electron density
!!  rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
!!  symq(4,2,nsym)=1 if symmetry preserves present qpoint. From symq3
!!  symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
!!  timrev=1 if time-reversal preserves the q wavevector; 0 otherwise.
!!  usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
!!  vtrial(nfftf,nspden)=GS potential (Hartree)
!!  vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree)
!!  vxcavg=average of vxc potential
!!  xred(3,natom)=reduced dimensionless atomic coordinates
!!
!! OUTPUT
!!  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
!!  ddkfil(3)=unit numbers for the three possible ddk files
!!  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives
!!  d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs
!!  d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs
!!  d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs
!!  etotal=total energy (sum of 8 contributions) (hartree)
!!
!! PARENTS
!!      respfn
!!
!! CHILDREN
!!      appdig,atm2fft3,ebands_free,ctocprj,distrb2,dtset_copy,dtset_free
!!      ebands_init,eig2stern,eigen_meandege,exit_check,fourdp,getcgqphase
!!      getmpw,getnel,getph,hdr_free,hdr_init,hdr_update,initmpi_band,initylmg
!!      insy3,inwffil,ioarr,ioddb8_out,kpgio,localfilnam,localrdfile
!!      localredirect,localwrfile,metric,mkcor3,mkrdim,mkrho3,outbsd,outgkk
!!      outwf,pawang_destroy,pawang_init,pawang_nullify,pawcprj_alloc
!!      pawcprj_copy,pawcprj_destroy,pawrhoij_alloc,pawrhoij_bcast
!!      pawrhoij_copy,pawrhoij_destroy,pawrhoij_nullify,prteigrs,prtene3,psddb8
!!      rotate_rho,scfcv3,set_pert_comm,set_pert_paw,setsym,setsymrhoij,status
!!      symkpt,timab,transgrid,unset_pert_comm,unset_pert_paw,vloca3,vlocalstr
!!      wffclose,wffopen,wfk_read_eigenvalues,wrtout,xmpi_bcast,xmpi_sum
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

subroutine loper3(atindx,atindx1,blkflg,codvsn,cpus,dimcprj,dim_eigbrd,dim_eig2nkq,doccde,&
&  ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,&
&  dyfr_cplex,dyfr_nondiag,d2bbb,d2lo,d2nl,d2ovl,eigbrd,eig2nkq,&
&  eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,&
&  etotal,fermie,gsqcut_eff,iexit,indsym,kxc,&
&  mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,&
&  nfftf,nhat,nkpt,nkxc,nspden,nsym,occ,&
&  paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
&  pertsy,prtbbb,psps,rfpert,rhog,rhor,symq,symrec,timrev,&
&  usecprj,vtrial,vxc,vxcavg,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,&
&  eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0)

 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use defs_wvltypes
 use m_profiling
 use m_xmpi
 use m_errors
 use m_exit
 use m_wffile
 use m_io_redirect
 use m_paral_pert
 use m_abi_etsf
 use m_ncfile
#ifdef HAVE_TRIO_NETCDF
 use netcdf
#endif

 use m_io_tools,   only : get_unit
 use m_fstrings,   only : strcat
 use m_ebands,     only : ebands_init, ebands_ncwrite, ebands_free, ebands_3dprint,get_nelect_per_spin
 use m_eig2d,      only : eigr2d_init,eigr2d_t, eigr2d_ncwrite,eigr2d_free
 use m_crystal,    only : crystal_init, crystal_free, crystal_t,isalchemical
 use m_crystal_io, only : crystal_ncwrite
 use m_ddb,        only : psddb8
 use m_header,     only : hdr_init, hdr_free, hdr_update
 use m_wfk,        only : wfk_read_eigenvalues
 use m_dtset,      only : dtset_copy, dtset_free
 use m_pawang,     only : pawang_type, pawang_nullify, pawang_init, pawang_destroy
 use m_pawrad,     only : pawrad_type
 use m_pawtab,     only : pawtab_type
 use m_paw_an,     only : paw_an_type
 use m_paw_ij,     only : paw_ij_type
 use m_pawfgrtab,  only : pawfgrtab_type
 use m_pawrhoij,   only : pawrhoij_type, pawrhoij_alloc, pawrhoij_destroy, pawrhoij_bcast, pawrhoij_copy, &
&                         pawrhoij_nullify,pawrhoij_redistribute
 use m_pawcprj,    only : pawcprj_type, pawcprj_alloc, pawcprj_destroy, pawcprj_copy
 use m_pawfgr,     only : pawfgr_type

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'loper3'
 use interfaces_14_hidewrite
 use interfaces_18_timing
 use interfaces_32_util
 use interfaces_41_geometry
 use interfaces_51_manage_mpi
 use interfaces_53_ffts
 use interfaces_56_recipspace
 use interfaces_62_iowfdenpot
 use interfaces_62_occeig
 use interfaces_65_nonlocal
 use interfaces_65_psp
 use interfaces_66_paw
 use interfaces_67_common
 use interfaces_72_response
 use interfaces_79_seqpar_mpi
!End of the abilint section

 implicit none

!Arguments ------------------------------------
 integer, intent(in) :: dim_eigbrd,dim_eig2nkq,dyfr_cplex,dyfr_nondiag,mk1mem,mkmem,mkqmem,mpert
 integer, intent(in) :: nfftf,nkpt,nkxc,nspden,nsym,prtbbb,timrev,usecprj
 integer, intent(out) :: iexit
 integer, intent(inout) :: my_natom
 real(dp), intent(in) :: cpus,gsqcut_eff,vxcavg
 real(dp), intent(inout) :: fermie
 real(dp), intent(inout) :: etotal !vz_i
 character(len=6), intent(in) :: codvsn
 type(MPI_type), intent(inout) :: mpi_enreg
 type(datafiles_type), intent(in) :: dtfil
 type(dataset_type), intent(in), target :: dtset
 type(pawang_type),intent(in) :: pawang
 type(pawfgr_type),intent(in) :: pawfgr
 type(pseudopotential_type), intent(inout) :: psps
 integer, intent(in) :: atindx(dtset%natom),atindx1(dtset%natom)
 integer, intent(in) :: dimcprj(dtset%natom*psps%usepaw),indsym(4,nsym,dtset%natom)
 integer, intent(in) :: nattyp(dtset%ntypat),pertsy(3,mpert)
 integer, intent(in) :: rfpert(mpert),symq(4,2,nsym),symrec(3,3,nsym)
 integer, intent(out) :: ddkfil(3)
 integer, intent(inout) :: blkflg(3,mpert,3,mpert) !vz_i
 integer, intent(out) :: clflg(3,mpert)
 real(dp), intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol)
 real(dp), intent(in) :: dyew(2,3,dtset%natom,3,dtset%natom)
 real(dp), intent(in) :: dyfrlo(3,3,dtset%natom)
 real(dp), intent(in) :: dyfrnl(dyfr_cplex,3,3,dtset%natom,1+(dtset%natom-1)*dyfr_nondiag)
 real(dp), intent(in) :: dyfrx1(2,3,dtset%natom,3,dtset%natom)
 real(dp), intent(in) :: dyfrx2(3,3,dtset%natom),eltcore(6,6)
 real(dp), intent(in) :: elteew(6+3*dtset%natom,6),eltfrhar(6,6)
 real(dp), intent(in) :: eltfrkin(6,6),eltfrloc(6+3*dtset%natom,6)
 real(dp), intent(in) :: eltfrnl(6+3*dtset%natom,6)
 real(dp), intent(in) :: eltfrxc(6+3*dtset%natom,6),kxc(nfftf,nkxc)
 real(dp), intent(in) :: nhat(nfftf,nspden)
 real(dp), intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol)
 real(dp), intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),vxc(nfftf,nspden)
 real(dp), intent(in) :: vtrial(nfftf,nspden)
 real(dp), intent(inout) :: xred(3,dtset%natom)
 real(dp), intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)!vz_i
 real(dp), intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i
 real(dp), intent(inout) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) !vz_i
 real(dp), intent(out) :: eigbrd(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eigbrd)
 real(dp), intent(out) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eig2nkq)
 type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:)
 type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:)
 type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:)
 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:)
 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 real(dp),pointer,optional :: eigen1_pert(:,:,:)
 real(dp),intent(out),optional :: occ_rbz_pert(:),eigen0_pert(:),eigenq_pert(:)
 real(dp),pointer,optional :: eigenq_fine(:,:,:)
 integer, intent(out),optional :: nkpt_rbz
 type(hdr_type),intent(out),optional :: hdr0,hdr_fine

!Local variables-------------------------------
!scalars
!Define file format for different type of files. Presently,
!only one file format is supported for each type of files, but this might
!change soon ...
!1   for wavefunction file, old format (version prior to 2.0)
!2   for wavefunction file, new format (version 2.0 and after)    (fform)
!(51 or 52   for density rho(r)       (fformr)
!101 or 102 for potential V(r) file. (fformv)
 integer,parameter :: fform=2,fformv=102,level=11,response=1
 integer :: fformr=52
 integer :: accessfil,ask_accurate,band_index,bantot,bantot_rbz,bdeigrf,bdtot1_index
 integer :: bdtot_index,choice,cplex,cplex_rhoij,ddkcase,dim_eig2rf,formeig
 integer :: fullinit,gscase,iband,icase,icase_eq,idir,idir0,idir_eq,ierr,ii,ikpt,ikpt1
 integer :: initialized,iorder_cprj,ipert,ipert_cnt,ipert_eq,ipert_me,ireadwf0
 integer :: iscf_mod,iscf_mod_save,isppol,istr,isym,master,mcg,mcgq,mcg1,mcprj,mcprjq
 integer :: me,mgfftf,mkmem_rbz,mk1mem_rbz,mkqmem_rbz,mpw,mpw1,mxfh
 integer :: n3xccc,nband_k,nblok,ncpgr,ndir,nkpt_eff,nkpt_max,nline_save,nmatel,npert_io,npert_me,nspden_rhoij
 integer :: nstep_save,nsym1,ntypat,nxfh,old_comm_atom,openexit,option,optorth,optthm,pertcase,rdwr
 integer :: rdwrpaw,spaceComm,smdelta,t_iostat,timrev_pert,unitout,useylmgr,useylmgr1,vrsddb
 real(dp) :: dosdeltae,eberry,ecore,ecut_eff,edocc,eei,eeig0,eew,efrhar,efrkin,efrloc
 real(dp) :: efrnl,efrx1,efrx2,ehart,ehart01,ehart1,eii,ek,ek0,ek1,ek2,eloc0
 real(dp) :: elpsp1,enl,enl0,enl1,entropy,enxc,eovl1,epaw1,exc1,fsum,gsqcut,maxocc,nelectkq
 real(dp) :: residm,tolwfr,tolwfr_save,toldfe_save,toldff_save,tolrff_save,tolvrs_save
 real(dp) :: ucvol
 logical,parameter :: paral_pert_inplace=.true.
 logical :: first_entry,found_eq_gkk,t_exist,paral_atom,remove_inv,write_1wfk
 character(len=fnlen) :: dscrpt,fiden1i,fiwf1i,fiwf1o,fiwfddk,gkkfilnam,fname
 character(len=500) :: message
 type(crystal_t) :: Crystal
 type(dataset_type), pointer :: dtset_tmp
 type(ebands_t) :: bs_rbz,Bands
 type(eigr2d_t)  :: eigr2d
 type(hdr_type) :: hdr
 type(ncfile_t) :: ncf
 type(pawang_type) :: pawang1
 type(wffile_type) :: wff1,wffddk,wffgs,wffkq,wffnow,wfftgs,wfftkq
 type(wvl_data) :: wvl
!arrays
 integer :: eq_symop(3,3),ngfftf(18)
 integer,allocatable :: blkflg_save(:,:,:,:),dyn(:),indkpt1(:),indkpt1_tmp(:),indsy1(:,:,:)
 integer,allocatable :: irrzon1(:,:,:),istwfk_rbz(:),istwfk_pert(:,:,:)
 integer,allocatable :: kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:),npwarr(:),npwtot(:)
 integer, allocatable :: npwtot1(:),npwar1_pert(:,:),npwarr_pert(:,:),npwtot_pert(:,:),pert_calc(:),pert_tmp(:)
 integer,allocatable :: symaf1(:),symaf1_tmp(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:)
 integer, pointer :: old_atmtab(:)
 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2)
 real(dp),allocatable :: buffer1(:,:,:,:,:),cg(:,:),cg1(:,:),cg1_active(:,:)
 real(dp),allocatable :: cg1_pert(:,:,:,:),cgq(:,:),gh0c1_pert(:,:,:,:)
 real(dp),allocatable :: doccde_rbz(:),docckqde(:)
 real(dp),allocatable :: gh1c_pert(:,:,:,:),eigen0(:),eigen0_copy(:),eigen1(:),eigen1_mean(:)
 real(dp),allocatable :: eigenq(:),gh1c_set(:,:),gh0c1_set(:,:),kpq(:,:)
 real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_pert(:),occ_rbz(:),occkq(:)
 real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons1(:,:,:),resid(:),rhog1(:,:)
 real(dp),allocatable :: rhor1_save(:,:,:)
 real(dp),allocatable :: rhor1(:,:),rho1wfg(:,:),rho1wfr(:,:),tnons1(:,:),tnons1_tmp(:,:)
 real(dp),allocatable :: vpsp1(:),work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:)
 real(dp),allocatable :: xfhist(:,:,:,:)
 real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:)
 real(dp),allocatable :: phasecg(:,:)
 type(pawcprj_type),allocatable :: cprj(:,:),cprjq(:,:)
 type(paw_ij_type),pointer :: paw_ij_pert(:)
 type(paw_an_type),pointer :: paw_an_pert(:)
 type(pawfgrtab_type),pointer :: pawfgrtab_pert(:)
 type(pawrhoij_type),allocatable :: pawrhoij1(:),pawrhoij_read(:)
 type(pawrhoij_type),pointer :: pawrhoij_pert(:)

! ***********************************************************************

 DBG_ENTER("COLL")

 _IBM6("In loper3")

 call timab(141,1,tsec)

!Structured debugging if prtvol==-level
 if(dtset%prtvol==-level)then
   write(message,'(80a,a,a)')  ('=',ii=1,80),ch10,' loper3 : enter , debug mode '
   call wrtout(std_out,message,'COLL')
 end if

!Obtain dimensional translations in reciprocal space gprimd,
!metrics and unit cell volume, from rprimd. Also output rprimd, gprimd and ucvol
 call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd)
 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol)

!Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90
 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
   mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:)
 else
   mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:)
 end if
 ecut_eff=dtset%ecut*(dtset%dilatmx)**2

!Various initializations/allocations
 iscf_mod=dtset%iscf
 ntypat=psps%ntypat
 nkpt_max=50;if (xmpi_paral==1) nkpt_max=-1
 paral_atom=(dtset%natom/=my_natom)
 cplex=2-timrev !cplex=2 ! DEBUG: impose cplex=2
 first_entry=.true.
 initialized=0
 call pawang_nullify(pawang1)
 ecore=zero ; ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero ; eii=zero
 clflg(:,:)=0 ! Array on calculated perturbations for eig2rf

!Save values of SCF cycle parameters
 iscf_mod_save = iscf_mod
 nstep_save = dtset%nstep
 nline_save = dtset%nline
 tolwfr_save = dtset%tolwfr
 toldfe_save = dtset%toldfe
 toldff_save = dtset%toldff
 tolrff_save = dtset%tolrff
 tolvrs_save = dtset%tolvrs

!This dtset will be used in scfcv3 to force non scf calculations for equivalent perturbations
 nullify(dtset_tmp)
 if (dtset%prepgkk/=0) then ! .and. dtset%use_nonscf_gkk==1) then !Later uncomment this - in scf case rhor1_save is used below only for testing
   ABI_DATATYPE_ALLOCATE(dtset_tmp,)
   call dtset_copy(dtset_tmp, dtset)
 else
   dtset_tmp => dtset
 end if

!If dtset%accesswff == 2 set all array outputs to netcdf format
 accessfil = 0
 if (dtset%accesswff == IO_MODE_NETCDF) accessfil = 1
 if (dtset%accesswff == IO_MODE_ETSF)   accessfil = 3

!Generate the 1-dimensional phases
 ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
 ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
 call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred)
 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
   call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
 else
   ph1df(:,:)=ph1d(:,:)
 end if

!Determine existence of pertubations and of pertubation symmetries
!Create array with pertubations which have to be calculated
 ABI_ALLOCATE(pert_tmp,(3*mpert))
 ipert_cnt=0
 do ipert=1,mpert
   do idir=1,3
     if( rfpert(ipert)==1 .and. dtset%rfdir(idir) == 1 )then
       if ((pertsy(idir,ipert)==1).or.&
&       ((dtset%prepanl == 1).and.(ipert == dtset%natom+2.or.ipert==dtset%natom+5)).or.&
&       ((dtset%prepgkk == 1).and.(ipert <= dtset%natom))  ) then
         ipert_cnt = ipert_cnt+1;
         pert_tmp(ipert_cnt) = idir+(ipert-1)*3;
       else
         write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
&         ' The perturbation idir=',idir,'  ipert=',ipert,' is',ch10,&
&         ' symmetric of a previously calculated perturbation.',ch10,&
&         ' So, its SCF calculation is not needed.',ch10
         call wrtout(std_out,message,'COLL')
         call wrtout(ab_out,message,'COLL')
       end if ! Test of existence of symmetry of perturbation
     end if ! Test of existence of perturbation
   end do
 end do
 ABI_ALLOCATE(pert_calc,(ipert_cnt))
 do icase=1,ipert_cnt
   pert_calc(icase)=pert_tmp(icase)
 end do
 ABI_DEALLOCATE(pert_tmp)

 if (dtset%prepgkk/=0) then ! .and. dtset%use_nonscf_gkk==1) then !Later uncomment this - in scf case rhor1_save is used below only for testing
   ABI_ALLOCATE(rhor1_save,(cplex*nfftf,nspden,ipert_cnt))
   rhor1_save=zero
   ABI_ALLOCATE(blkflg_save,(3,mpert,3,mpert))
 end if

! Initialize quantities for netcdf print 
 ABI_ALLOCATE(eigen0_copy,(dtset%mband*nkpt*dtset%nsppol))
 eigen0_copy(:)=zero

!%%%% Parallelization over perturbations %%%%%
!*Define file output/log file names
 npert_io=ipert_cnt;if (dtset%nppert<=1) npert_io=0
 call localfilnam(mpi_enreg%comm_pert,mpi_enreg%comm_cell_pert,mpi_enreg%comm_world,dtfil%filnam_ds,'_PRT',npert_io)
! Compute the number of perturbation done by the current cpu
 if(mpi_enreg%paral_pert==1) then
   npert_me = 0
   ipert_me = 0
   do icase=1,ipert_cnt
     if (mpi_enreg%distrb_pert(icase)==mpi_enreg%me_pert) then
       npert_me = npert_me +1
     endif
   enddo
 endif

!*Redefine communicators
 call set_pert_comm(mpi_enreg,dtset%nppert)

!*Redistribute PAW on-site data
 if (paral_pert_inplace) then
   call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
&   paw_an,paw_ij,pawfgrtab,pawrhoij)
   pawfgrtab_pert=>pawfgrtab ; pawrhoij_pert=>pawrhoij
   paw_an_pert   =>paw_an    ; paw_ij_pert  =>paw_ij
 else
   call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
&   paw_an,paw_ij,pawfgrtab,pawrhoij,&
&   paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,&
&   pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert)
 end if

!Loop on perturbations
!==========================================================================
 do icase=1,ipert_cnt
   _IBM6("In loop on perts")

!  %%%% Parallelization over perturbations %%%%%
!  Select the perturbations treated by curent processor
   if(mpi_enreg%paral_pert==1) then
     if (mpi_enreg%distrb_pert(icase)/=mpi_enreg%me_pert) cycle
   end if
!  *Redefine output/log files
   call localwrfile(mpi_enreg%comm_cell,icase,npert_io,mpi_enreg%paral_pert,0)

!  Retrieve type and direction of the perturbation
   if (pert_calc(icase) <= dtset%natom*3) then
     idir = mod(pert_calc(icase),3)
     if (idir==0) idir=3
     ipert=( (pert_calc(icase)-idir) / 3 + 1)
   else
     ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1
     idir = mod(pert_calc(icase),3)
     if (idir==0) idir=3
   end if
   pertcase=idir+(ipert-1)*3

!  Init MPI communicator
   spaceComm=mpi_enreg%comm_cell
   me=mpi_enreg%me_kpt

!  ===== Describe the perturbation in output/log file
   _IBM6("IBM6 before print perturbation")

   write(message, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,&
&   ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:)
   call wrtout(std_out,message,'COLL')
   call wrtout(ab_out,message,'COLL')
   if(ipert>=1 .and. ipert<=dtset%natom)then
     write(message, '(a,i4,a,i4)' )' Perturbation : displacement of atom',ipert,'   along direction',idir
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
     if(iscf_mod == -3)then
       write(message, '(a,a,a,a,a,a,a,a)' )ch10,&
&       ' loper3 : COMMENT -',ch10,&
&       '  The first-order density is imposed to be zero (iscf=-3).',ch10,&
&       '  Although this is strange in the case of phonons,',ch10,&
&       '  you are allowed to do so.'
       call wrtout(std_out,message,'COLL')
       call wrtout(ab_out,message,'COLL')
     end if
   else if(ipert==dtset%natom+1)then
     write(message,'(a,i4)')' Perturbation : derivative vs k along direction',idir
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
     if( iscf_mod /= -3 )then
       write(message, '(4a)' )ch10,&
&       ' loper3 : COMMENT -',ch10,&
&       '  In a d/dk calculation, iscf is set to -3 automatically.'
       call wrtout(std_out,message,'COLL')
       call wrtout(ab_out,message,'COLL')
       iscf_mod=-3
     end if
     if( abs(dtset%sciss) > 1.0d-8 )then
       write(message, '(a,a,a,a,f14.8,a,a)' )ch10,&
&       ' loper3 : WARNING -',ch10,&
&       '  Value of sciss=',dtset%sciss,ch10,&
&       '  Scissor with d/dk calculation : you are using a "naive" approach !'
       call wrtout(std_out,message,'COLL')
       call wrtout(ab_out,message,'COLL')
     end if
   else if(ipert==dtset%natom+2)then
     write(message, '(a,i4)' )' Perturbation : homogeneous electric field along direction',idir
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
     if( iscf_mod == -3 )then
       write(message, '(a,a,a,a,a,a)' )ch10,&
&       ' loper3 : COMMENT -',ch10,&
&       '  The first-order density is imposed to be zero (iscf=-3).',ch10,&
&       '  This corresponds to a calculation without local fields.'
       call wrtout(std_out,message,'COLL')
       call wrtout(ab_out,message,'COLL')
     end if
   else if(ipert==dtset%natom+5)then
     write(message, '(a,i4)' )&
&     ' Perturbation : homogeneous magnetic field along direction, presently set to electric field for testing',idir
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
     if( iscf_mod == -3 )then
       write(message, '(a,a,a,a,a,a)' )ch10,&
&       ' loper3 : COMMENT -',ch10,&
&       '  The first-order density is imposed to be zero (iscf=-3).',ch10,&
&       '  This corresponds to a calculation without local fields.'
       call wrtout(std_out,message,'COLL')
       call wrtout(ab_out,message,'COLL')
     end if
   else if(ipert>dtset%natom+7 .or. ipert<=0 )then
     write(message, '(a,i4,a,a,a)' ) &
&     '  ipert=',ipert,' is outside the [1,dtset%natom+7] interval.',ch10,&
&     '  This perturbation is not (yet) allowed.'
     MSG_BUG(message)
   end if
!  Initialize the diverse parts of energy :
   eew=zero ; efrloc=zero ; efrnl=zero ; efrx1=zero ; efrx2=zero
   efrhar=zero ; efrkin=zero
   if(ipert<=dtset%natom)then
     eew=dyew(1,idir,ipert,idir,ipert)
     efrloc=dyfrlo(idir,idir,ipert)
     if (dyfr_nondiag==0) efrnl=dyfrnl(1,idir,idir,ipert,1)
     if (dyfr_nondiag/=0) efrnl=dyfrnl(1,idir,idir,ipert,ipert)
     efrx1=dyfrx1(1,idir,ipert,idir,ipert)
     efrx2=dyfrx2(idir,idir,ipert)
   else if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
!    istr = 1,2,...,6 and indicates the cartesian strain component
     if(ipert==dtset%natom+3) then
       istr=idir
     else
       istr=idir+3
     end if
     eii=eltcore(istr,istr)
     eew=elteew(istr,istr)
     efrhar=eltfrhar(istr,istr)
     efrkin=eltfrkin(istr,istr)
     efrloc=eltfrloc(istr,istr)
     efrnl=eltfrnl(istr,istr)
     efrx1=eltfrxc(istr,istr)
   end if

!  Determine the subset of symmetry operations (nsym1 operations)
!  that leaves the perturbation invariant, and initialize corresponding arrays
!  symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW)..
   ABI_ALLOCATE(symaf1_tmp,(nsym))
   ABI_ALLOCATE(symrl1_tmp,(3,3,nsym))
   ABI_ALLOCATE(tnons1_tmp,(3,nsym))
   if (dtset%prepanl/=1.and.&
&   dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
&   dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17) then
     call insy3(gprimd,idir,indsym,ab_out,ipert,dtset%natom,nsym,nsym1,2,&
&     dtset%symafm,symaf1_tmp,symq,symrec,dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp)
   else
     nsym1 = 1
     symaf1_tmp(1) = 1
     symrl1_tmp(:,:,1) = dtset%symrel(:,:,1)
     tnons1_tmp(:,1) = 0_dp
   end if
   ABI_ALLOCATE(indsy1,(4,nsym1,dtset%natom))
   ABI_ALLOCATE(symrc1,(3,3,nsym1))
   ABI_ALLOCATE(symaf1,(nsym1))
   ABI_ALLOCATE(symrl1,(3,3,nsym1))
   ABI_ALLOCATE(tnons1,(3,nsym1))
   symaf1(1:nsym1)=symaf1_tmp(1:nsym1)
   symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1)
   tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1)
   ABI_DEALLOCATE(symaf1_tmp)
   ABI_DEALLOCATE(symrl1_tmp)
   ABI_DEALLOCATE(tnons1_tmp)

!  Set up corresponding symmetry data
   ABI_ALLOCATE(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4)))
   ABI_ALLOCATE(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4)))
   call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,&
&   nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)
   if (psps%usepaw==1) then
!    Allocate/initialize only zarot in pawang1 datastructure
     call pawang_init(pawang1,0,pawang%l_max-1,0,nsym1,0,1,0,0,0)
     call setsymrhoij(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot)
   end if

!  Initialize k+q array
   ABI_ALLOCATE(kpq,(3,nkpt))
   if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
     kpq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran
   else
     do ikpt=1,nkpt
       kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt)
     end do
   end if

!  Determine the subset of k-points needed in the "reduced Brillouin zone",
!  and initialize other quantities
   ABI_ALLOCATE(indkpt1_tmp,(nkpt))
   ABI_ALLOCATE(wtk_folded,(nkpt))
   indkpt1_tmp(:)=0 ; optthm=0
   timrev_pert=timrev
   if(dtset%ieig2rf>0) then
     timrev_pert=0
     call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,&
&     1,symrc1,timrev_pert,dtset%wtk,wtk_folded)
   else
!    For the time being, the time reversal symmetry is not used
!    for ddk, elfd, mgfd perturbations.
     timrev_pert=timrev
     if(ipert==dtset%natom+1.or.ipert==dtset%natom+2.or.ipert==dtset%natom+5.or.&
&     dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.  &
&     dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) timrev_pert=0
     call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,&
     nsym1,symrc1,timrev_pert,dtset%wtk,wtk_folded)
   end if
!  Be careful: when parallelization over perturbation is activated, mkmem/mk1mem is not modified
   mkmem_rbz=mkmem   ; if (mpi_enreg%nproc_pert>1) mkmem_rbz =nkpt_rbz
   mkqmem_rbz=mkqmem ; if (mpi_enreg%nproc_pert>1) mkqmem_rbz=nkpt_rbz
   mk1mem_rbz=mk1mem ; if (mpi_enreg%nproc_pert>1) mk1mem_rbz=nkpt_rbz

   ABI_ALLOCATE(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
   ABI_ALLOCATE(indkpt1,(nkpt_rbz))
   ABI_ALLOCATE(istwfk_rbz,(nkpt_rbz))
   ABI_ALLOCATE(kpq_rbz,(3,nkpt_rbz))
   ABI_ALLOCATE(kpt_rbz,(3,nkpt_rbz))
   ABI_ALLOCATE(nband_rbz,(nkpt_rbz*dtset%nsppol))
   ABI_ALLOCATE(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
   ABI_ALLOCATE(wtk_rbz,(nkpt_rbz))
   indkpt1(:)=indkpt1_tmp(1:nkpt_rbz)
   do ikpt=1,nkpt_rbz
     istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt))
     kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt))
     kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt))
     wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt))
   end do
   ABI_DEALLOCATE(indkpt1_tmp)
   ABI_DEALLOCATE(wtk_folded)

!  Transfer occ to occ_rbz and doccde to doccde_rbz :
!  this is a more delicate issue
!  NOTE : this takes into account that indkpt1 is ordered
!  MG: What about using occ(band,kpt,spin) ???
   bdtot_index=0;bdtot1_index=0
   do isppol=1,dtset%nsppol
     ikpt1=1
     do ikpt=1,nkpt
       nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
!      Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1)
       if(ikpt1/=nkpt_rbz+1)then
         if(ikpt==indkpt1(ikpt1))then
           nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k
           occ_rbz(1+bdtot1_index:nband_k+bdtot1_index)    = occ(1+bdtot_index:nband_k+bdtot_index)
           doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index) = doccde(1+bdtot_index:nband_k+bdtot_index)
           ikpt1=ikpt1+1
           bdtot1_index=bdtot1_index+nband_k
         end if
       end if
       bdtot_index=bdtot_index+nband_k
     end do
   end do

   _IBM6("IBM6 in loper3 before getmpw")

!  Compute maximum number of planewaves at k
   call timab(142,1,tsec)
   call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz)
   call timab(142,2,tsec)

!  Allocate some k-dependent arrays at k
   ABI_ALLOCATE(kg,(3,mpw*mkmem_rbz))
   ABI_ALLOCATE(npwarr,(nkpt_rbz))
   ABI_ALLOCATE(npwtot,(nkpt_rbz))

!  Determine distribution of k-points/bands over MPI processes
   if(xmpi_paral==1) then
     ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol))
     if (allocated(mpi_enreg%my_kpttab)) then
       ABI_DEALLOCATE(mpi_enreg%my_kpttab)
     end if
     ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz))
     call distrb2(dtset%mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg)
   end if
   call initmpi_band(mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol)

   _IBM6("IBM6 before kpgio")

!  Set up the basis sphere of planewaves at k
   call timab(143,1,tsec)
   call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,dtfil%fnametmp_kg,&
&   kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw,npwarr,npwtot,dtset%nsppol,dtfil%unkg)
   call timab(143,2,tsec)

!  Set up the spherical harmonics (Ylm) at k
   useylmgr=0
   if (psps%useylm==1.and. &
&   (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. &
&   (psps%usepaw==1.and.(ipert==dtset%natom+2.or.ipert==dtset%natom+5)))) &
&   useylmgr=1
   ABI_ALLOCATE(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
   ABI_ALLOCATE(ylmgr,(mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
   if (psps%useylm==1) then
     option=0
     if (useylmgr==1) option=1
     call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,&
&     npwarr,dtset%nsppol,option,rprimd,dtfil%unkg,dtfil%unylm,ylm,ylmgr)
   end if

   _IBM6("Before ieig2rf > 0")

!  Set up occupations for this perturbation
   if (dtset%ieig2rf>0) then
     if (.not.allocated(istwfk_pert)) then
       ABI_ALLOCATE(istwfk_pert,(nkpt,3,mpert))
       ABI_ALLOCATE(occ_pert,(dtset%mband*nkpt*dtset%nsppol))
       istwfk_pert(:,:,:)=0 ; occ_pert(:)=zero
     end if
     istwfk_pert(:,idir,ipert)=istwfk_rbz(:)
     occ_pert(:)=occ_rbz(:)
   end if

!  Print a separator in output file
   write(message,'(3a)')ch10,'--------------------------------------------------------------------------------',ch10
   call wrtout(ab_out,message,'COLL')

!  Initialize band structure datatype at k
   bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol))
   ABI_ALLOCATE(eigen0,(bantot_rbz))
   eigen0(:)=zero
   call ebands_init(bantot_rbz,bs_rbz,dtset%nelect,doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,&
&   nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,&
&   dtset%occopt,occ_rbz,wtk_rbz)
   ABI_DEALLOCATE(eigen0)

!  Initialize header, update it with evolving variables
   gscase=0 ! A GS WF file is read
   call hdr_init(bs_rbz,codvsn,dtset,hdr0,pawtab,gscase,psps,wvl%descr,&
&   mpi_comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)

   call hdr_update(bantot_rbz,etotal,fermie,hdr0,dtset%natom,&
&   residm,rprimd,occ_rbz,pawrhoij_pert,psps%usepaw,xred, &
&   mpi_comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)

!  Clean band structure datatype (should use it more in the future !)
   call ebands_free(bs_rbz)

   _IBM6("before inwffil")

!  Initialize GS wavefunctions at k
   ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0
   mcg=mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol
   ABI_ALLOCATE(cg,(2,mcg))
   ABI_CHECK_ALLOC("out-of-memory in cg")
   ABI_ALLOCATE(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol))
   call timab(144,1,tsec)
   call status(0,dtfil%filstat,iexit,level,'call inwffil-k')
   call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,&
&   formeig,gmet,hdr0,ireadwf0,istwfk_rbz,kg,&
&   kpt_rbz,dtset%localrdwf,dtset%mband,mcg,&
&   mkmem_rbz,mpi_enreg,mpw,nband_rbz,dtset%ngfft,nkpt_rbz,npwarr,&
&   dtset%nsppol,nsym,occ_rbz,optorth,rprimd,dtset%symafm,&
&   dtset%symrel,dtset%tnons,dtfil%unkg,wffgs,wfftgs,&
&   dtfil%unwffgs,dtfil%unwftgs,dtfil%fnamewffk,dtfil%fnametmp_wfgs,wvl)
   call timab(144,2,tsec)
!  Close wffgs%unwff, if it was ever opened (in inwffil)
   if (ireadwf0==1) then
     call WffClose(wffgs,ierr)
   end if

!  PAW: compute on-site projections of GS wavefunctions (cprj) (and derivatives) at k
   ncpgr=0
   ABI_DATATYPE_ALLOCATE(cprj,(0,0))
   if (psps%usepaw==1) then
     ncpgr=3 ! Valid for ipert<=natom (phonons), ipert=natom+2 (elec. field) or ipert=natom+5 (magn. field)
     if (ipert==dtset%natom+1) ncpgr=1
     if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) ncpgr=1
     if (usecprj==1) then
       mcprj=dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol
       ABI_DATATYPE_DEALLOCATE(cprj)
       ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj))
       call pawcprj_alloc(cprj,ncpgr,dimcprj)
       if (ipert<=dtset%natom) then
         choice=2; iorder_cprj=0; idir0=0
       else if (ipert==dtset%natom+1) then
         choice=5; iorder_cprj=0; idir0=idir
       else if (ipert==dtset%natom+2) then
         choice=5; iorder_cprj=0; idir0=0
       else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
         choice=3; iorder_cprj=0; idir0=istr
       else
         choice=1; iorder_cprj=0; idir0=idir
       end if
       call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,-1,idir0,iorder_cprj,istwfk_rbz,&
&       kg,kpt_rbz,dtset%mband,mcg,mcprj,dtset%mgfft,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,&
&       dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,&
&       npwarr,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,&
&       rmet,dtset%typat,ucvol,dtfil%unpaw,dtfil%unkg,dtfil%unylm,useylmgr,wfftgs,xred,ylm,ylmgr)
     end if
   end if

!  Compute maximum number of planewaves at k+q
!  Will be useful for both GS wfs at k+q and RF wavefunctions
   call timab(143,1,tsec)
   call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpq_rbz,mpi_enreg,mpw1,nkpt_rbz)
   call timab(143,2,tsec)

!  Allocate some arrays at k+q
   ABI_ALLOCATE(kg1,(3,mpw1*mk1mem_rbz))
   ABI_ALLOCATE(npwar1,(nkpt_rbz))
   ABI_ALLOCATE(npwtot1,(nkpt_rbz))

!  Set up the basis sphere of planewaves at k+q
!  Will be useful for both GS wfs at k+q and RF wavefunctions
   call timab(142,1,tsec)
   call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1,dtfil%fnametmp_kg1,&
&   kpq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1,&
&   npwar1,npwtot1,dtset%nsppol,dtfil%unkg1)
   call timab(142,2,tsec)

!  Set up the spherical harmonics (Ylm) at k+q
   useylmgr1=0
   if (psps%useylm==1.and. &
&   (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. &
&   (psps%usepaw==1.and.(ipert==dtset%natom+2.or.ipert==dtset%natom+5)))) &
   useylmgr1=1
   ABI_ALLOCATE(ylm1,(mpw1*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
   ABI_ALLOCATE(ylmgr1,(mpw1*mk1mem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
   if (psps%useylm==1) then
     option=0
     if (useylmgr1==1) option=1
     call initylmg(gprimd,kg1,kpq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1,nband_rbz,nkpt_rbz,&
&     npwar1,dtset%nsppol,option,rprimd,dtfil%unkg1,dtfil%unylm1,ylm1,ylmgr1)
   end if

!  Print a separator in output file
   write(message, '(a,a)' )'--------------------------------------------------------------------------------',ch10
   call wrtout(ab_out,message,'COLL')

!  Initialize band structure datatype at k+q
   ABI_ALLOCATE(eigenq,(bantot_rbz))
   eigenq(:)=zero
   call ebands_init(bantot_rbz,bs_rbz,dtset%nelect,doccde_rbz,eigenq,istwfk_rbz,kpq_rbz,&
&   nband_rbz,nkpt_rbz,npwar1,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz)
   ABI_DEALLOCATE(eigenq)
!  Initialize header
   call hdr_init(bs_rbz,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, &
&   mpi_comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )

!  Clean band structure datatype (should use it more in the future !)
   call ebands_free(bs_rbz)

!  Initialize wavefunctions at k+q
!  MG: Here it is possible to avoid the extra reading if the same k mesh can be used.
   ireadwf0=1 ; formeig=0 ; ask_accurate=1 ; optorth=0
   mcgq=mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol
   ABI_ALLOCATE(cgq,(2,mcgq))
   ABI_CHECK_ALLOC("out-of-memory in cgq")
   ABI_ALLOCATE(eigenq,(dtset%mband*nkpt_rbz*dtset%nsppol))
   call timab(144,1,tsec)
   call status(0,dtfil%filstat,iexit,level,'call inwffilkq')
   call inwffil(ask_accurate,cgq,dtset,dtset%ecut,ecut_eff,eigenq,dtset%exchn2n3d,&
&   formeig,gmet,hdr,&
&   ireadwf0,istwfk_rbz,kg1,kpq_rbz,dtset%localrdwf,dtset%mband,mcgq,&
&   mkqmem_rbz,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,&
&   dtset%nsppol,nsym,occ_rbz,optorth,&
&   rprimd,dtset%symafm,dtset%symrel,dtset%tnons,&
&   dtfil%unkg1,wffkq,wfftkq,dtfil%unwffkq,dtfil%unwftkq,dtfil%fnamewffq,dtfil%fnametmp_wfkq,wvl)
   call timab(144,2,tsec)
!  Close dtfil%unwffkq, if it was ever opened (in inwffil)
   if (ireadwf0==1) then
     call WffClose(wffkq,ierr)
   end if

!  PAW: compute on-site projections of GS wavefunctions (cprjq) (and derivatives) at k+q
   ABI_DATATYPE_ALLOCATE(cprjq,(0,0))
   if (psps%usepaw==1) then
     if (usecprj==1) then
       call status(0,dtfil%filstat,iexit,level,'call ctocprjq')
       mcprjq=dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol
       ABI_DATATYPE_DEALLOCATE(cprjq)
       ABI_DATATYPE_ALLOCATE(cprjq,(dtset%natom,mcprjq))
       call pawcprj_alloc(cprjq,0,dimcprj)
       if (ipert<=dtset%natom.and.(sum(dtset%qptn(1:3)**2)>=1.d-14)) then ! phonons at non-zero q
         choice=1 ; iorder_cprj=0 ; idir0=0
         call ctocprj(atindx,cgq,choice,cprjq,gmet,gprimd,-1,idir0,0,istwfk_rbz,&
&         kg1,kpq_rbz,dtset%mband,mcgq,mcprjq,dtset%mgfft,mkqmem_rbz,mpi_enreg,psps%mpsang,mpw1,&
&         dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,&
&         npwar1,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,&
&         psps,rmet,dtset%typat,ucvol,dtfil%unpawq,dtfil%unkg1,dtfil%unylm1,useylmgr1,&
&         wfftkq,xred,ylm1,ylmgr1)
       else
         call pawcprj_copy(cprj,cprjq)
       end if
     end if
   end if

!  ===== Report on eigenq values
   if (dtset%ieig2rf>0.and.icase==ipert_cnt) then
     eigen0_pert(:) = eigen0(:)
     eigenq_pert(:) = eigenq(:)
     occ_rbz_pert(:) = occ_rbz(:)
   end if
   call wrtout(std_out,ch10//' loper3 : eigenq array',"COLL")
   nkpt_eff=nkpt
   if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. nkpt>nkpt_max ) nkpt_eff=nkpt_max
   band_index=0
   do isppol=1,dtset%nsppol
     do ikpt=1,nkpt_rbz
       nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
       if(ikpt<=nkpt_eff)then
         write(message, '(a,i2,a,i5)' )'  isppol=',isppol,', k point number',ikpt
         call wrtout(std_out,message,'COLL')
         do iband=1,nband_k,4
           write(message, '(a,4es16.6)')'  ',eigenq(iband+band_index:min(iband+3,nband_k)+band_index)
           call wrtout(std_out,message,'COLL')
         end do
       else if(ikpt==nkpt_eff+1)then
         write(message,'(a,a)' )'  respfn : prtvol=0, 1 or 2, stop printing eigenq.',ch10
         call wrtout(std_out,message,'COLL')
       end if
       band_index=band_index+nband_k
     end do
   end do

!  Generate occupation numbers for the reduced BZ at k+q
   ABI_ALLOCATE(docckqde,(dtset%mband*nkpt_rbz*dtset%nsppol))
   ABI_ALLOCATE(occkq,(dtset%mband*nkpt_rbz*dtset%nsppol))
   if(0<=dtset%occopt .and. dtset%occopt<=2)then
!    Same occupation numbers at k and k+q (usually, insulating)
     occkq(:)=occ_rbz(:)
     docckqde(:)=zero  ! docckqde is irrelevant in this case
   else
!    Metallic occupation numbers
     option=1
     dosdeltae=zero ! the DOS is not computed with option=1
     maxocc=two/(dtset%nspinor*dtset%nsppol)
     call getnel(docckqde,dosdeltae,eigenq,entropy,fermie,maxocc,dtset%mband,&
&     nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occkq,dtset%occopt,option,&
&     dtset%tphysel,dtset%tsmear,tmp_unit,wtk_rbz)
!    Compare nelect at k and nelelect at k+q
     write(message, '(a,a,a,es16.6,a,es16.6,a)')&
&     ' loper3 : total number of electrons, from k and k+q',ch10,&
&     '  fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.'
     call wrtout(ab_out,message,'COLL')
     call wrtout(std_out,message,'COLL')
   end if

!  Debug message
   if(dtset%prtvol==-level)then
     call wrtout(std_out,'loper3 : initialisation of q part done. ','COLL')
   end if

!  Initialisation of first-order wavefunctions
   write(message,'(3a,i4)')&
&   ' Initialisation of the first-order wave-functions :',ch10,&
&   '  ireadwf=',dtfil%ireadwf
   call wrtout(std_out,message,'COLL')
   call wrtout(ab_out,message,'COLL')
   call appdig(pertcase,dtfil%fnamewff1,fiwf1i)
   call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o)

!  Allocate 1st-order PAW occupancies (rhoij1)
   if (psps%usepaw==1) then
     ABI_DATATYPE_ALLOCATE(pawrhoij1,(my_natom))
     call pawrhoij_nullify(pawrhoij1)
     cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=dtset%nspden
     call pawrhoij_alloc(pawrhoij1,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
&     dtset%typat,pawtab=pawtab,mpi_comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
     if (cplex_rhoij/=hdr%pawrhoij(1)%cplex) then
!      Eventually reallocate hdr%pawrhoij
       call pawrhoij_destroy(hdr%pawrhoij)
       call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
&       dtset%typat,pawtab=pawtab,mpi_comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
     end if
   else
     ABI_DATATYPE_ALLOCATE(pawrhoij1,(0))
   end if

!  Initialize 1st-order wavefunctions
   formeig=1; ask_accurate=0; optorth=0
! NB: 4 Sept 2013: this was introducing a bug - for ieig2rf ==0 the dim was being set to 1 and passed to vtowfk3
   dim_eig2rf=0
   if (dtset%ieig2rf > 0 .and. dtset%ieig2rf/=2) then
     dim_eig2rf=1
   end if
   mcg1=mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol
   ABI_ALLOCATE(cg1,(2,mcg1))
   ABI_CHECK_ALLOC("out of memory in cg1")
   ABI_ALLOCATE(cg1_active,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
   ABI_ALLOCATE(gh1c_set,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
   ABI_ALLOCATE(gh0c1_set,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
!  XG090606 This is needed in the present 5.8.2 , for portability for the pathscale machine.
!  However, it is due to a bug to be corrected by Paul Boulanger. When the bug will be corrected,
!  this line should be removed.
   if(mk1mem_rbz/=0 .and. dtset%ieig2rf/=0)then
     cg1_active=zero
     gh1c_set=zero
     gh0c1_set=zero
   end if
   ABI_ALLOCATE(eigen1,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol))
   ABI_ALLOCATE(resid,(dtset%mband*nkpt_rbz*dtset%nsppol))
   call timab(144,1,tsec)
   call status(pertcase,dtfil%filstat,iexit,level,'call inwffil  ')
   call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
&   formeig,gmet,hdr,&
&   dtfil%ireadwf,istwfk_rbz,kg1,kpq_rbz,dtset%localrdwf,&
&   dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,&
&   dtset%nsppol,nsym1,occ_rbz,optorth,rprimd,&
&   symaf1,symrl1,tnons1,dtfil%unkg1,wff1,wffnow,dtfil%unwff1,dtfil%unwft1,&
&   fiwf1i,dtfil%fnametmp_1wf1,wvl)
   call timab(144,2,tsec)
!  Close wff1 (filename fiwf1i), if it was ever opened (in inwffil)
   if (dtfil%ireadwf==1) then
     call WffClose(wff1,ierr)
   end if
!  Eventually reytrieve 1st-order PAW occupancies from file header
   if (psps%usepaw==1.and.dtfil%ireadwf/=0) then
     call pawrhoij_copy(hdr%pawrhoij,pawrhoij1,mpi_comm_atom=mpi_enreg%comm_atom , &
&     mpi_atmtab=mpi_enreg%my_atmtab)
   end if

!  In case of electric field, open the ddk wf file
   if (ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and. &
&   (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and. &
&   dtset%berryopt/= 7.and.dtset%berryopt/=14.and. &
&   dtset%berryopt/=16.and.dtset%berryopt/=17)) then
     ddkcase=idir+dtset%natom*3
     call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk)
     write(message,'(2a)')'-loper3 : read the ddk wavefunctions from file: ',fiwfddk
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
     master=0
     call WffOpen(dtset%accesswff,spaceComm,fiwfddk,ierr,wffddk,master,me,dtfil%unddk)
   end if

!  In case of magnetic field, open the ddk wf file
   if (ipert==dtset%natom+5.and.sum( (dtset%qptn(1:3))**2)<1.0d-7.and. &
&   (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and. &
&   dtset%berryopt/= 7.and.dtset%berryopt/=14.and. &
&   dtset%berryopt/=16.and.dtset%berryopt/=17)) then
     ddkcase=idir+dtset%natom*3
     call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk)
     write(message,'(2a)')'-loper3 : read the ddk wavefunctions from file: ',fiwfddk
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
     master=0
     call WffOpen(dtset%accesswff,spaceComm,fiwfddk,ierr,wffddk,master,me,dtfil%unddk)
   end if

!  Get first-order local potentials and 1st-order core correction density change
!  (do NOT include xccc3d1 in vpsp1 : this will be done in scfcv3 because vpsp1
!  might become spin-polarized)

!  The correct value of gsqcut should be derived from ecutf and not ecutsm, but the present trick works
   gsqcut=gsqcut_eff

   n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
   ABI_ALLOCATE(xccc3d1,(cplex*n3xccc))
   ABI_ALLOCATE(vpsp1,(cplex*nfftf))

!  PAW: compute Vloc(1) and core(1) together in reciprocal space
!  --------------------------------------------------------------
   if (psps%usepaw==1) then
     ndir=1
     call atm2fft3(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
&     mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
&     ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,&
&     atmrho1=xccc3d1,atmvloc1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,pawtab=pawtab,vspl=psps%vlspl)
   else

!    Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
!    ------------------------------------------------------------------------------

     if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
!      Section for strain perturbation
       call vlocalstr(gmet,gprimd,gsqcut,istr,mgfftf,mpi_enreg,&
&       psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,dtset%paral_kgb,ph1df,psps%qgrid_vl,&
&       ucvol,psps%vlspl,vpsp1)
     else
       call vloca3(atindx,cplex,gmet,gsqcut,idir,ipert,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
&       nattyp,nfftf,ngfftf,ntypat,ngfftf(1),ngfftf(2),ngfftf(3),dtset%paral_kgb,ph1df,psps%qgrid_vl,&
&       dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
     end if

     if(psps%n1xccc/=0)then
       call mkcor3(cplex,idir,ipert,dtset%natom,ntypat,ngfftf(1),psps%n1xccc,&
&       ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
     end if ! psps%n1xccc/=0
   end if ! usepaw

   eigen1(:)=zero; resid(:)=zero

!  Get starting charge density and Hartree + xc potential
   ABI_ALLOCATE(rhor1,(cplex*nfftf,nspden))
   ABI_ALLOCATE(rhog1,(2,nfftf))

!  can we get this set of gkk matrices from previously calculated rhog1 through a non-scf calculation?
   found_eq_gkk=.false.
   if (dtset%prepgkk == 1 .and. ipert <= dtset%natom) then
!    NOTE: this does not take into account combinations e.g. x+y -> z
!    if rhor1 add linearly this could be done...
     do icase_eq = 1, icase-1
       idir_eq = mod(icase_eq,3)
       if (idir_eq==0) idir_eq=3
       ipert_eq = ( (icase_eq-idir_eq) / 3 + 1)

! find sym which links old perturbation to present one
       do isym=1, nsym
         ! check that isym preserves qpt to begin with
         if (symq(4,1,isym) /= 1 .or. &
&         symq(1,1,isym) /= 0 .or. &
&         symq(2,1,isym) /= 0 .or. &
&         symq(3,1,isym) /= 0      ) cycle

         eq_symop = dtset%symrel(:,:,isym)
         if (indsym(4,isym,ipert) == ipert_eq .and. &
&         abs(eq_symop(idir,idir_eq)) == 1 .and. &
&         sum(abs(eq_symop(:,idir_eq))) == 1) then
           found_eq_gkk = .true.
           exit
         end if
       end do ! isym
       if (found_eq_gkk) exit
     end do ! icase_eq
   end if ! check for prepgkk with symmetric pert

   if (found_eq_gkk) then
     write (message, '(a,l6,i6,2a,3i6,2a,3i6,a)')  &
&     ' found_eq_gkk,isym = ', found_eq_gkk, isym, ch10, &
&     ' idir,  ipert,  icase   =  ', idir, ipert, icase, ch10, &
&     ' idireq iperteq icaseeq =  ', idir_eq, ipert_eq, icase_eq, ch10
     call wrtout(std_out,message,'COLL')
!
!    Make density for present perturbation, which is symmetric of 1 or more previous perturbations:
!    rotate 1DEN arrays with symrel(isym) to produce rhog1_eq rhor1_eq
!
     if (dtset%use_nonscf_gkk == 1) then
       call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, &
&       dtset%paral_kgb, rhor1_save(:,:,icase_eq), rhog1, rhor1, eq_symop, dtset%tnons(:,isym))

       rhor1 = rhor1 * eq_symop(idir,idir_eq)

! TODO: rotate rhoij in PAW case

       blkflg_save = blkflg
       dtset_tmp%iscf = -2
       iscf_mod = -2
       dtset_tmp%nstep = 1
       dtset_tmp%nline = 1
       if (abs(dtset_tmp%tolwfr) < 1.e-24) dtset_tmp%tolwfr = 1.e-24
       dtset_tmp%toldfe = zero
       dtset_tmp%toldff = zero
       dtset_tmp%tolrff = zero
       dtset_tmp%tolvrs = zero
       write (message, '(a,i6,a)') ' NOTE: doing GKK calculation for icase ', icase, ' with non-SCF calculation'
       call wrtout(std_out,message,'COLL')
       !call wrtout(ab_out,message,'COLL') ! decomment and update output files

     else ! do not use non-scf shortcut, but save rotated 1DEN for comparison
! saves the rotated rho, for later comparison with the full SCF rhor1: comment lines below for iscf = -2
       call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, &
&       dtset%paral_kgb, rhor1_save(:,:,icase_eq), rhog1, rhor1_save(:,:,icase), eq_symop, &
&       dtset%tnons(:,isym))
       rhor1_save(:,:,icase) = rhor1_save(:,:,icase) * eq_symop(idir,idir_eq)

     end if ! force non scf calculation of other gkk, or not
   end if ! found an equiv perturbation for the gkk

   if ( (dtfil%ireadwf==0 .and. iscf_mod/=-4 .and. dtset%get1den==0 .and. dtset%ird1den==0) .or. (iscf_mod== -3 ) ) then

     rhor1(:,:)=zero ; rhog1(:,:)=zero
!    PAW: rhoij have been set to zero in call to pawrhoij_alloc above

   else  ! rhor1 not being forced to 0.0
     if(iscf_mod>0)then
!      cplex=2 gets the complex density, =1 only real part
       if (psps%usepaw==1) then
!        Be careful: in PAW, rho does not include the 1st-order compensation
!        density (to be added in scfcv3.F90) !
         ABI_ALLOCATE(rho1wfg,(2,dtset%nfft))
         ABI_ALLOCATE(rho1wfr,(dtset%nfft,nspden))
         call mkrho3(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,&
&         kg,kg1,dtset%mband,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,&
&         dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
&         occ_rbz,dtset%paral_kgb,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,ucvol,&
&         wtk_rbz)
         call transgrid(cplex,mpi_enreg,nspden,+1,1,1,dtset%paral_kgb,pawfgr,rho1wfg,rhog1,rho1wfr,rhor1)
         ABI_DEALLOCATE(rho1wfg)
         ABI_DEALLOCATE(rho1wfr)
       else
         call mkrho3(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,&
&         kg,kg1,dtset%mband,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,&
&         dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
&         occ_rbz,dtset%paral_kgb,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,ucvol,&
&         wtk_rbz)
       end if

     else if (.not. found_eq_gkk) then ! negative iscf_mod and no symmetric rotation of rhor1
       rdwr=1;rdwrpaw=psps%usepaw;if(dtfil%ireadwf/=0) rdwrpaw=0
       if (me==0) then
!        Read rho1(r) from a disk file
         if (rdwrpaw/=0) then
           ABI_DATATYPE_ALLOCATE(pawrhoij_read,(dtset%natom))
!          Values of nspden/nsppol are dummy ones;
!          they will be overwritten later (pawrhoij_bcast below)
           cplex_rhoij=max(cplex,dtset%pawcpxocc)
           call pawrhoij_nullify(pawrhoij_read)
           call pawrhoij_alloc(pawrhoij_read,cplex_rhoij,hdr%nspden,hdr%nspinor,&
&           hdr%nsppol,dtset%typat,pawtab=pawtab)
         end if
         call appdig(pertcase,dtfil%fildens1in,fiden1i)
         call ioarr(accessfil,rhor1, dtset, etotal,fformr,fiden1i,hdr, mpi_enreg,&
&         cplex*nfftf,pawrhoij_read,rdwr,rdwrpaw,wvl%den)
       end if
!      MPI communication
       call xmpi_bcast(rhor1,0,spaceComm,ierr)
       if (rdwrpaw/=0) then
         call pawrhoij_bcast(pawrhoij_read,pawrhoij1,0,spaceComm,mpi_comm_atom=mpi_enreg%comm_atom)
         call pawrhoij_destroy(pawrhoij_read)
         ABI_DATATYPE_DEALLOCATE(pawrhoij_read)
       end if

!      Compute up+down rho1(G) by fft
       ABI_ALLOCATE(work,(cplex*nfftf))
       work(:)=rhor1(:,1)
       call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
       ABI_DEALLOCATE(work)

     end if ! rhor1 generated or read in from file

   end if ! rhor1 set to 0 or read in from file

!  Check whether exiting was required by the user.
!  If found then do not start minimization steps
   openexit=1 ; if(dtset%chkexit==0) openexit=0
   call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
!  If immediate exit, and wavefunctions were not read, must zero eigenvalues
   if (iexit/=0) eigen1(:)=zero

   if (iexit==0) then
     _IBM6("before scfcv3")

!    Main calculation: get 1st-order wavefunctions from Sternheimer equation (SCF cycle)
     call scfcv3(atindx,atindx1,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,&
&     dimcprj,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,&
&     d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
&     ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,&
&     enl0,enl1,eovl1,epaw1,etotal,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,indkpt1,&
&     indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
&     kg,kg1,kpt_rbz,kxc,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,&
&     mpert,mpi_enreg,mpw,mpw1,my_natom,&
&     nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,&
&     npwarr,npwar1,nspden,&
&     nsym1,n3xccc,occkq,occ_rbz,&
&     paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,&
&     pertcase,phnons1,ph1d,ph1df,prtbbb,psps,&
&     dtset%qptn,resid,residm,rhog,rhog1,&
&     rhor,rhor1,rprimd,symaf1,symrc1,symrl1,&
&     usecprj,useylmgr,useylmgr1,wffddk,vpsp1,vtrial,vxc,&
&     wtk_rbz,wvl%den,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1)

     _IBM6("after scfcv3")

!    2nd-order eigenvalues stuff
     if (dtset%ieig2rf>0) then
       if (first_entry) then
         nullify(eigen1_pert)
         first_entry = .false.
       end if
       if (.not.associated(eigen1_pert)) then
         ABI_ALLOCATE(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert))
         ABI_ALLOCATE(cg1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
         ABI_ALLOCATE(gh0c1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
         ABI_ALLOCATE(gh1c_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
         ABI_ALLOCATE(npwarr_pert,(nkpt_rbz,mpert))
         ABI_ALLOCATE(npwar1_pert,(nkpt_rbz,mpert))
         ABI_ALLOCATE(npwtot_pert,(nkpt_rbz,mpert))
         eigen1_pert(:,:,:) = zero
         cg1_pert(:,:,:,:) = zero
         gh0c1_pert(:,:,:,:) = zero
         gh1c_pert(:,:,:,:) = zero
         npwar1_pert (:,:) = 0
         npwarr_pert (:,:) = 0
       end if
       clflg(idir,ipert)=1
       eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:)
       if (dtset%ieig2rf==1.or.dtset%ieig2rf==3.or.dtset%ieig2rf==4) then
         cg1_pert(:,:,idir,ipert)=cg1_active(:,:)
         gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:)
         gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:)
       end if
       npwarr_pert(:,ipert)=npwarr(:)
       npwar1_pert(:,ipert)=npwar1(:)
       npwtot_pert(:,ipert)=npwtot(:)
     end if
     ABI_DEALLOCATE(gh1c_set)
     ABI_DEALLOCATE(gh0c1_set)
     ABI_DEALLOCATE(cg1_active)

   end if ! End of the check of hasty exit

   call timab(146,1,tsec)

!  Print out message at the end of the iterations
   write(message, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,&
&   ' ----iterations are completed or convergence reached----',ch10
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,  message,'COLL')

!  Update the content of the header (evolving variables)
   call hdr_update(bantot_rbz,etotal,fermie,hdr,dtset%natom,&
&   residm,rprimd,occ_rbz,pawrhoij1,psps%usepaw,xred,mpi_comm_atom=mpi_enreg%comm_atom , &
&   mpi_atmtab=mpi_enreg%my_atmtab )

!  Print _gkk file for this perturbation
   if (dtset%prtgkk == 1) then
     call appdig(3*(ipert-1)+idir,dtfil%fnameabo_gkk,gkkfilnam)
     nmatel = dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol
     ABI_ALLOCATE(phasecg, (2, nmatel))
     call getcgqphase(dtset, timrev, cg,  mcg,  cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg)
     phasecg(1,:) = one
     phasecg(2,:) = zero
! NB: phasecg not actually used in outgkk for the moment (2013/08/15)
     call outgkk(bantot_rbz, nmatel,gkkfilnam,eigen0,eigen1,hdr0,hdr,mpi_enreg,phasecg)
     ABI_DEALLOCATE(phasecg)
   end if

   if (dtset%prepgkk == 1 .and. found_eq_gkk) then
     if (dtset%use_nonscf_gkk == 1) then
!      Restore old values of SCF cycle parameters
       iscf_mod = iscf_mod_save
       dtset_tmp%iscf = iscf_mod_save
       dtset_tmp%nstep = nstep_save
       dtset_tmp%nline = nline_save
       dtset_tmp%tolwfr = tolwfr_save
       dtset_tmp%toldfe = toldfe_save
       dtset_tmp%toldff = toldff_save
       dtset_tmp%tolrff = tolrff_save
       dtset_tmp%tolvrs = tolvrs_save
       blkflg = blkflg_save ! this ensures we do not use the (unconverged) 2DTE from this non scf run
!      Save density for present perturbation, for future use in symmetric perturbations
       rhor1_save(:,:,icase) = rhor1
     else
       write (message, '(a,3E20.10)') 'norm diff = ', sum(abs(rhor1_save(:,:,icase) - rhor1)), &
&       sum(abs(rhor1)), sum(abs(rhor1_save(:,:,icase) - rhor1))/sum(abs(rhor1))
       call wrtout(std_out,  message,'COLL')
     end if
   end if

   write_1wfk = .True.
   !if (dtset%prtwf==-1) then
   ! Write wavefunctions file only if convergence was not achieved.
   ! MG TODO
   ! write_1wfk = (scfcv3_retcode /= 0) 
   !end if

   if (write_1wfk) then
     ! Output 1st-order wavefunctions in file
     mxfh=0; nxfh=0
     ABI_ALLOCATE(xfhist,(3,dtset%natom+4,2,mxfh))
     call status(pertcase,dtfil%filstat,iexit,level,'call outwf    ')
     call outwf(cg1,dtset,eigen1,fiwf1o,hdr,kg1,kpt_rbz,&
&     dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,mxfh,dtset%natom,nband_rbz,&
&     nkpt_rbz,npwar1,dtset%nsppol,dtset%nstep,&
&     nxfh,occ_rbz,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr,xfhist)
     ABI_DEALLOCATE(xfhist)
   end if

!  If the perturbation is d/dk, evaluate the f-sum rule.
   if (ipert==dtset%natom+1 )then
!    Note : the factor of two is related to the difference
!    between Taylor expansion and perturbation expansion
!    Note : this expression should be modified for ecutsm.
!    Indeed, the present one will NOT tend to 1.0_dp.
     ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect
     fsum=-ek1/ek2
     if(dtset%ecutsm<tol6)then
       write(message, '(a,es20.10,a,a,es20.10)' ) &
&       ' loper3 : ek2=',ek2,ch10,&
&       '          f-sum rule ratio=',fsum
     else
       write(message, '(a,es20.10,a,a,es20.10,a)' ) &
&       ' loper3 : ek2=',ek2,ch10,&
&       '          f-sum rule ratio=',fsum,' (note : ecutsm/=0)'
     end if
     call wrtout(std_out,message,'COLL')
     call wrtout(ab_out,message,'COLL')
!    Write the diagonal elements of the dH/dk operator, after averaging over degenerate states
     ABI_ALLOCATE(eigen1_mean,(dtset%mband*nkpt_rbz*dtset%nsppol))
     call eigen_meandege(eigen0,eigen1,eigen1_mean,dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,1)
     option=4
     call prteigrs(eigen1_mean,dtset%enunit,fermie,dtfil%fnametmp_1wf1_eig,ab_out,iscf_mod,kpt_rbz,dtset%kptopt,&
&     dtset%mband,nband_rbz,nkpt_rbz,dtset%nnsclo,dtset%nsppol,occ_rbz,dtset%occopt,&
&     option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk_rbz)
     ABI_DEALLOCATE(eigen1_mean)
   end if

!  Print the energies
   if (dtset%nline/=0 .or. dtset%nstep/=0)then
     call prtene3(dtset%berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
&     ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,eovl1,epaw1,exc1,ab_out,&
&     ipert,dtset%natom,psps%usepaw)
   end if

   if(mpi_enreg%paral_pert==1) then
     if (ipert_me < npert_me -1) then
       call hdr_free(hdr0)
     else   
       eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
     endif
     ipert_me = ipert_me +1
   else
     if (icase == ipert_cnt) then
       eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
     else
       call hdr_free(hdr0)
     end if
   endif

!Release the temporary arrays (for k, k+q and 1st-order)
   ABI_DEALLOCATE(cg)
   ABI_DEALLOCATE(cgq)
   ABI_DEALLOCATE(cg1)
   ABI_DEALLOCATE(docckqde)
   ABI_DEALLOCATE(doccde_rbz)
   ABI_DEALLOCATE(eigen0)
   ABI_DEALLOCATE(eigenq)
   ABI_DEALLOCATE(eigen1)
   ABI_DEALLOCATE(kpq)
   ABI_DEALLOCATE(indkpt1)
   ABI_DEALLOCATE(indsy1)
   ABI_DEALLOCATE(istwfk_rbz)
   ABI_DEALLOCATE(irrzon1)
   ABI_DEALLOCATE(kg)
   ABI_DEALLOCATE(kg1)
   ABI_DEALLOCATE(kpq_rbz)
   ABI_DEALLOCATE(kpt_rbz)
   ABI_DEALLOCATE(nband_rbz)
   ABI_DEALLOCATE(npwarr)
   ABI_DEALLOCATE(npwar1)
   ABI_DEALLOCATE(npwtot)
   ABI_DEALLOCATE(npwtot1)
   ABI_DEALLOCATE(occkq)
   ABI_DEALLOCATE(occ_rbz)
   ABI_DEALLOCATE(phnons1)
   ABI_DEALLOCATE(resid)
   ABI_DEALLOCATE(rhog1)
   ABI_DEALLOCATE(rhor1)
   ABI_DEALLOCATE(symaf1)
   ABI_DEALLOCATE(symrc1)
   ABI_DEALLOCATE(symrl1)
   ABI_DEALLOCATE(tnons1)
   ABI_DEALLOCATE(wtk_rbz)
   ABI_DEALLOCATE(xccc3d1)
   ABI_DEALLOCATE(vpsp1)
   ABI_DEALLOCATE(ylm)
   ABI_DEALLOCATE(ylm1)
   ABI_DEALLOCATE(ylmgr)
   ABI_DEALLOCATE(ylmgr1)
   if (psps%usepaw==1) then
     call pawang_destroy(pawang1)
     call pawrhoij_destroy(pawrhoij1)
     if (usecprj==1) then
       call pawcprj_destroy(cprj)
       call pawcprj_destroy(cprjq)
     end if
   end if
   ABI_DATATYPE_DEALLOCATE(pawrhoij1)
   ABI_DATATYPE_DEALLOCATE(cprjq)
   ABI_DATATYPE_DEALLOCATE(cprj)
   if(xmpi_paral==1)  then
     ABI_DEALLOCATE(mpi_enreg%proc_distrb)
     ABI_DEALLOCATE(mpi_enreg%my_kpttab)
   end if
   call hdr_free(hdr)

!  %%%% Parallelization over perturbations %%%%%
!  *Redefine output/log files
   call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,npert_io,&
&   mpi_enreg%paral_pert,0)

   call timab(146,2,tsec)
   if(iexit/=0) exit
 end do ! End loop on perturbations

!%%%% Parallelization over perturbations %%%%%
!*Restore default communicators
 call unset_pert_comm(mpi_enreg)
!*Gather output/log files
 ABI_ALLOCATE(dyn,(npert_io))
 if (npert_io>0) dyn=1
 call localrdfile(mpi_enreg%comm_pert,mpi_enreg%comm_world,.true.,npert_io,&
& mpi_enreg%paral_pert,0,dyn)
 ABI_DEALLOCATE(dyn)
!*Restore PAW on-site data
 if (paral_pert_inplace) then
   call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
&   paw_an,paw_ij,pawfgrtab,pawrhoij)
 else
   call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
&   paw_an,paw_ij,pawfgrtab,pawrhoij,&
&   paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,&
&   pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert)
 end if

!#################################################################################
!Calculate the second-order eigenvalues for a wavevector Q

 call timab(147,1,tsec)
 smdelta = dtset%smdelta
 bdeigrf = dtset%bdeigrf
 if(dtset%bdeigrf == -1) bdeigrf = dtset%mband

 if(dtset%ieig2rf > 0) then

   if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
&   (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
     call wrtout(std_out,'Reading the dense grid WF file',"COLL")
!    We get the Abinit header of the file hdr_fine as ouput
!    We get eigenq_fine(mband,hdr_fine%nkpt,hdr_fine%nsppol) as ouput
     call wfk_read_eigenvalues(dtfil%fnameabi_wfkfine,eigenq_fine,hdr_fine,mpi_enreg%comm_world)
     ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband")
   end if
   if(dtset%kptopt==3 .or. dtset%kptopt==0)then
     if (dtset%nsym > 1) then
       MSG_ERROR("Symmetries are not implemented for temperature dependence calculations")
     endif
     write(std_out,*) 'Entering: eig2stern'
     if(smdelta>0)then
       if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
&       (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
         call eig2stern(bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
&         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
&         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
&         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,eigbrd,dtset,eigenq_fine,hdr_fine,hdr0)
       else
         call eig2stern(bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
&         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
&         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
&         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,eigbrd)
       end if
     else
       if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
&       (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
         call eig2stern(bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
&         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
&         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
&         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,eigbrd,dtset,eigenq_fine,hdr_fine,hdr0)
       else
         call eig2stern(bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
&         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
&         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
&         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta)
       end if
     end if
     write(std_out,*) 'Leaving: eig2stern'

     if (dtset%ieig2rf==1.or.dtset%ieig2rf==2) then
       master=0
       if (me==master) then
!        print _EIGR2D file for this perturbation
         unitout = dtfil%unddb
         vrsddb=100401
         dscrpt=' Note : temporary (transfer) database '
!        tolwfr must be initialized here, but it is a dummy value
         tolwfr=1.0_dp
         call ioddb8_out (dscrpt,dtfil%fnameabo_eigr2d,dtset%natom,dtset%mband,&
&         dtset%nkpt,dtset%nsym,dtset%ntypat,dtfil%unddb,vrsddb,&
&         dtset%acell_orig(1:3,1),dtset%amu_orig(:,1),dtset%dilatmx,dtset%ecut,dtset%ecutsm,&
&         dtset%intxc,dtset%iscf,dtset%ixc,dtset%kpt,dtset%kptnrm,&
&         dtset%natom,dtset%nband,dtset%ngfft,dtset%nkpt,dtset%nspden,dtset%nspinor,&
&         dtset%nsppol,dtset%nsym,dtset%ntypat,occ_pert,dtset%occopt,dtset%pawecutdg,&
&         dtset%rprim_orig(1:3,1:3,1),dtset%sciss,dtset%spinat,dtset%symafm,dtset%symrel,&
&         dtset%tnons,tolwfr,dtset%tphysel,dtset%tsmear,&
&         dtset%typat,dtset%usepaw,dtset%wtk,xred,psps%ziontypat,dtset%znucl)
         nblok=1 ; fullinit=1 ; choice=2
         call psddb8 (choice,psps%dimekb,psps%ekb,fullinit,psps%indlmn,&
&         psps%lmnmax,nblok,ntypat,dtfil%unddb,pawtab,&
&         psps%pspso,psps%usepaw,psps%useylm,vrsddb)
         call outbsd(bdeigrf,dtset,eig2nkq,dtset%natom,nkpt_rbz,unitout)
!        print _EIGI2D file for this perturbation
         if(smdelta>0) then
           unitout = dtfil%unddb
           vrsddb=100401
           dscrpt=' Note : temporary (transfer) database '
!          tolwfr must be initialized here, but it is a dummy value
           tolwfr=1.0_dp
           call ioddb8_out (dscrpt,dtfil%fnameabo_eigi2d,dtset%natom,dtset%mband,&
&           dtset%nkpt,dtset%nsym,dtset%ntypat,dtfil%unddb,vrsddb,&
&           dtset%acell_orig(1:3,1),dtset%amu_orig(:,1),dtset%dilatmx,dtset%ecut,dtset%ecutsm,&
&           dtset%intxc,dtset%iscf,dtset%ixc,dtset%kpt,dtset%kptnrm,&
&           dtset%natom,dtset%nband,dtset%ngfft,dtset%nkpt,dtset%nspden,dtset%nspinor,&
&           dtset%nsppol,dtset%nsym,dtset%ntypat,occ_pert,dtset%occopt,dtset%pawecutdg,&
&           dtset%rprim_orig(1:3,1:3,1),dtset%sciss,dtset%spinat,dtset%symafm,dtset%symrel,&
&           dtset%tnons,tolwfr,dtset%tphysel,dtset%tsmear,&
&           dtset%typat,dtset%usepaw,dtset%wtk,xred,psps%ziontypat,dtset%znucl)
           nblok=1 ; fullinit=1 ; choice=2
           call psddb8 (choice,psps%dimekb,psps%ekb,fullinit,psps%indlmn,&
&           psps%lmnmax,nblok,ntypat,dtfil%unddb,pawtab,&
&           psps%pspso,psps%usepaw,psps%useylm,vrsddb)
           call outbsd(bdeigrf,dtset,eigbrd,dtset%natom,nkpt_rbz,unitout)
         end if !smdelta

!        Output of the EIGR2D.nc file.
         fname = strcat(dtfil%filnam_ds(4),"_EIGR2D.nc")
!        Crystalline structure.
         remove_inv=.false.
         if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true.

         call crystal_init(Crystal,dtset%spgroup,dtset%natom,dtset%npsp,psps%ntypat, &
&         dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
&         dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr0%title,&
&         dtset%symrel,dtset%tnons,dtset%symafm)

         if (isalchemical(Crystal)) then
           MSG_WARNING("Alchemical pseudos are not supported by ETSF-IO, EIGR2D file won't be produced")
         else
!          Electronic band energies.
           bantot= dtset%mband*dtset%nkpt*dtset%nsppol
           call ebands_init(bantot,Bands,dtset%nelect,doccde,eigen0_copy,hdr0%istwfk,hdr0%kptns,&
&           hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,&
&           hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk)

!Second order derivative EIG2D (real and Im)
           call eigr2d_init(eig2nkq,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
#ifdef HAVE_TRIO_NETCDF
           NCF_CHECK(ncfile_create(ncf,fname,NF90_CLOBBER),"Creating EIGR2D file")
           call crystal_ncwrite(Crystal,ncf%ncid)
           call ebands_ncwrite(Bands,dtset%nshiftk_orig,dtset%shiftk_orig,dtset%nshiftk,dtset%shiftk,&
&               dtset%ngkpt,dtset%kptrlatt,ncf%ncid)
           call eigr2d_ncwrite(eigr2d,dtset%qptn(:),ncf%ncid)
           NCF_CHECK(ncfile_close(ncf),"Closing EIGR2D file")
#else
           ABI_UNUSED(ncf%ncid)
#endif
         end if
       end if
     end if !ieig2rf==1.or.ieig2rf==2
   else
     write(message,'(3a)')&
&     'K point grids must be the same for every perturbation: eig2stern not called',ch10,&
&     'Action: Put kptopt=3 '
     MSG_WARNING(message)
   end if !kptopt
   ABI_DEALLOCATE(gh1c_pert)
   ABI_DEALLOCATE(gh0c1_pert)
   ABI_DEALLOCATE(cg1_pert)
   ABI_DEALLOCATE(istwfk_pert)
   ABI_DEALLOCATE(npwarr_pert)
   ABI_DEALLOCATE(npwar1_pert)
   ABI_DEALLOCATE(npwtot_pert)
   ABI_DEALLOCATE(occ_pert)
 end if  !if dtset%ieig2rf
!Free memory.
 ABI_DEALLOCATE(eigen0_copy)
 call hdr_free(hdr0)
 call crystal_free(Crystal)
 call ebands_free(Bands)
 call eigr2d_free(eigr2d)

 call timab(147,2,tsec)
!######################################################################################

!Get ddk file information, for later use in dyout3
 ddkfil(:)=0
 do idir=1,3
   ddkcase=idir+dtset%natom*3
   call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk)
!  Check that ddk file exists
   inquire(file=fiwfddk,iostat=t_iostat,exist=t_exist)
!  If the file exists set ddkfil to a non-zero value
   if (t_exist) then
     ddkfil(idir)=20+idir
   end if
 end do

 ABI_DEALLOCATE(ph1d)
 ABI_DEALLOCATE(ph1df)
 ABI_DEALLOCATE(pert_calc)

!destroy dtset_tmp
 if (dtset%prepgkk /= 0) then ! .and. dtset%use_nonscf_gkk == 1) then !Later uncomment this - in scf case rhor1_save is used below only for testing
   ABI_DEALLOCATE(rhor1_save)
   ABI_DEALLOCATE(blkflg_save)
   call dtset_free (dtset_tmp)
   ABI_DATATYPE_DEALLOCATE(dtset_tmp)
 end if

!In paral_pert-case some array's have to be reconstructed
 if(mpi_enreg%paral_pert==1) then
   ABI_ALLOCATE(buffer1,(2,3,mpert,3,mpert*(2+psps%usepaw)))
   buffer1(:,:,:,:,1:mpert)=d2lo(:,:,:,:,:)
   buffer1(:,:,:,:,1+mpert:2*mpert)=d2nl(:,:,:,:,:)
   if (psps%usepaw==1) then
     buffer1(:,:,:,:,1+2*mpert:3*mpert)=d2ovl(:,:,:,:,:)
   end if
   call xmpi_sum(buffer1,mpi_enreg%comm_pert,ierr)
   call xmpi_sum(blkflg,mpi_enreg%comm_pert,ierr)
   if(dtset%prtbbb==1) then
     call xmpi_sum(d2bbb,mpi_enreg%comm_pert,ierr)
   end if
   d2lo(:,:,:,:,:)=buffer1(:,:,:,:,1:mpert)
   d2nl(:,:,:,:,:)=buffer1(:,:,:,:,1+mpert:2*mpert)
   if (psps%usepaw==1) then
     d2ovl(:,:,:,:,:)=buffer1(:,:,:,:,1+2*mpert:3*mpert)
   end if
   ABI_DEALLOCATE(buffer1)
 end if

 call timab(141,2,tsec)

 DBG_EXIT("COLL")

end subroutine loper3
!!***

