TABLE OF CONTENTS


ABINIT/dfpt_init_mag1 [ Functions ]

[ Top ] [ Functions ]

NAME

  dfpt_init_mag1

FUNCTION

  Initial guess of the first order magnetization/density for magnetic field perturbation.
  The first order magnetization is set so as to zero out the first order XC magnetic field, which
  should minimize the second order XC energy (without taking self-consistency into account).

INPUTS

  ipert = perturbtation type (works only for ipert==natom+5)
  idir  = direction of the applied magnetic field
  cplex = complex or real first order density and magnetization
  nfft  = dimension of the fft grid
  nspden= number of density matrix components
  nkxc  = number of kxc components
  vxc0(nfft,nspden)  = GS XC potential
  kxc0(nfft,nspden)  = GS XC derivatives
  rhor0(nfft,nspden) = GS density matrix

OUTPUT

  rhor1(cplex*nfft) = first order density magnetization guess

PARENTS

      dfpt_looppert

CHILDREN

SOURCE

3197 subroutine dfpt_init_mag1(ipert,idir,rhor1,rhor0,cplex,nfft,nspden,vxc0,kxc0,nkxc)
3198 
3199 
3200 !This section has been created automatically by the script Abilint (TD).
3201 !Do not modify the following lines by hand.
3202 #undef ABI_FUNC
3203 #define ABI_FUNC 'dfpt_init_mag1'
3204 !End of the abilint section
3205 
3206  implicit none
3207 
3208 !Arguments ------------------------------------
3209  integer , intent(in)    :: ipert,idir,cplex,nfft,nspden,nkxc
3210  real(dp), intent(in)    :: vxc0(nfft,nspden),rhor0(nfft,nspden)
3211  real(dp), intent(in)    :: kxc0(nfft,nkxc)
3212  real(dp), intent(out)   :: rhor1(cplex*nfft,nspden)
3213 
3214 !Local variables-------------------------------
3215  integer  :: ipt
3216  real(dp) :: bxc0,bxc1
3217  real(dp) :: m1_norm,m0_norm
3218  real(dp) :: f_dot_m
3219  real(dp) :: mdir(3),fdir(3)
3220 
3221 ! *************************************************************************
3222 
3223  if (nspden==2) then
3224 
3225    if(cplex==1) then
3226      do ipt=1,nfft
3227        bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2)) ! d/dm Bxc
3228        !this overestimates the first order magnetization because of n1 not taken into account
3229        m1_norm=-half*(1/bxc1)
3230        rhor1(ipt,1)=zero             ! rho_up+rho_dwn    => charge density
3231        rhor1(ipt,2)=half*m1_norm     ! rho_up=1/2(rho+m) => half*m
3232      end do
3233    else
3234      do ipt=1,cplex*nfft
3235        rhor1(ipt,:)=zero
3236      end do
3237    end if
3238 
3239  else if(nspden==4) then
3240 
3241    fdir=zero
3242    fdir(idir)= 1.0d0
3243    do ipt=1,nfft
3244      m0_norm=sqrt(rhor0(ipt,2)**2+rhor0(ipt,3)**2+rhor0(ipt,4)**2)
3245      mdir(1)=rhor0(ipt,2)/m0_norm
3246      mdir(2)=rhor0(ipt,3)/m0_norm
3247      mdir(3)=rhor0(ipt,4)/m0_norm
3248      f_dot_m=fdir(1)*mdir(1)+fdir(2)*mdir(2)+fdir(3)*mdir(3) ! projection of the field direction on m0
3249 
3250      bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2))  ! d/dm Bxc
3251      m1_norm=(-half/bxc1)*f_dot_m                            ! get an estimate of the norm of m1
3252 
3253      bxc0=-sqrt((half*(vxc0(ipt,1)-vxc0(ipt,2)))**2+vxc0(ipt,3)**2+vxc0(ipt,4)**2)
3254      if(cplex==1) then
3255        rhor1(ipt,1)=zero       ! rho_up+rho_dwn    => charge density
3256        rhor1(ipt,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1))   ! m1x
3257        rhor1(ipt,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2))   ! m1y
3258        rhor1(ipt,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3))   ! m1z
3259        rhor1(ipt,:)=zero
3260      else
3261        rhor1(2*ipt-1,1)=zero       ! Re rho_up+rho_dwn
3262        rhor1(2*ipt-1,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1))   ! m1x
3263        rhor1(2*ipt-1,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2))   ! m1x
3264        rhor1(2*ipt-1,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3))   ! m1x
3265        rhor1(2*ipt  ,1)=zero       ! Im rho_up+rho_dwn
3266        rhor1(2*ipt  ,2)=zero
3267        rhor1(2*ipt  ,3)=zero
3268        rhor1(2*ipt  ,4)=zero
3269 
3270        rhor1(2*ipt-1,1)=zero; rhor1(2*ipt,1)=zero
3271        rhor1(2*ipt-1,2)=zero; rhor1(2*ipt,2)=zero
3272        rhor1(2*ipt-1,3)=zero; rhor1(2*ipt,3)=zero
3273        rhor1(2*ipt-1,4)=zero; rhor1(2*ipt,4)=zero
3274 
3275      end if
3276    end do
3277  end if
3278 
3279 end subroutine dfpt_init_mag1

ABINIT/dfpt_looppert [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_looppert

FUNCTION

 Loop over perturbations

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  codvsn=code version
  cpus=cpu time limit in seconds
  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
  dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D 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 tensor
  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
  eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor
  fermie=fermi energy (Hartree)
  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 rhotoxc.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 littlegroup_q
  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
  usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
  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,atom_gauss,crystal_free,crystal_init,ctocprj,ddb_free
      ddb_from_file,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write,dfpt_atm2fft
      dfpt_init_mag1,dfpt_mkcore,dfpt_mkrho,dfpt_prtene,dfpt_scfcv
      dfpt_vlocal,disable_timelimit,distrb2,dtset_copy,dtset_free,ebands_free
      ebands_init,efmas_main,efmas_analysis,eig2stern,eigen_meandege,eigr2d_free,eigr2d_init
      eigr2d_ncwrite,exit_check,fourdp,getcgqphase,getcut,getmpw,getnel,getph
      gkk_free,gkk_init,gkk_ncwrite,hdr_copy,hdr_free,hdr_init,hdr_update
      initmpi_band,initylmg,inwffil,kpgio,littlegroup_pert,localfilnam
      localrdfile,localredirect,localwrfile,metric,mkrdim,outbsd,outgkk,outwf
      pawang_free,pawang_init,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawcprj_getdim,pawrhoij_alloc,pawrhoij_copy,pawrhoij_free
      pawrhoij_nullify,prteigrs,put_eneocc_vect,read_rhor,rf2_getidirs
      rotate_rho,set_pert_comm,set_pert_paw,setsym,setsym_ylm,status,symkpt
      timab,transgrid,unset_pert_comm,unset_pert_paw,vlocalstr,wffclose
      wfk_open_read,wfk_read_eigenvalues,wrtout,xmpi_sum

SOURCE

 217 subroutine dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,&
 218 &  ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,&
 219 &  dyfr_cplex,dyfr_nondiag,d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,&
 220 &  eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
 221 &  etotal,fermie,iexit,indsym,kxc,&
 222 &  mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,&
 223 &  nfftf,nhat,nkpt,nkxc,nspden,nsym,occ,&
 224 &  paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
 225 &  pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,&
 226 &  usecprj,usevdw,vtrial,vxc,vxcavg,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,&
 227 &  eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0)
 228 
 229 
 230 !This section has been created automatically by the script Abilint (TD).
 231 !Do not modify the following lines by hand.
 232 #undef ABI_FUNC
 233 #define ABI_FUNC 'dfpt_looppert'
 234 !End of the abilint section
 235 
 236  implicit none
 237 
 238 !Arguments ------------------------------------
 239  integer, intent(in) :: dim_eigbrd,dim_eig2nkq,dyfr_cplex,dyfr_nondiag,mk1mem,mkmem,mkqmem,mpert
 240  integer, intent(in) :: nfftf,nkpt,nkxc,nspden,nsym,prtbbb,timrev,usecprj,usevdw
 241  integer, intent(out) :: iexit
 242  integer, intent(inout) :: my_natom
 243  real(dp), intent(in) :: cpus,vxcavg
 244  real(dp), intent(inout) :: fermie
 245  real(dp), intent(inout) :: etotal
 246  character(len=6), intent(in) :: codvsn
 247  type(MPI_type), intent(inout) :: mpi_enreg
 248  type(datafiles_type), intent(in) :: dtfil
 249  type(dataset_type), intent(in), target :: dtset
 250  type(pawang_type),intent(in) :: pawang
 251  type(pawfgr_type),intent(in) :: pawfgr
 252  type(pseudopotential_type), intent(inout) :: psps
 253  integer, intent(in) :: atindx(dtset%natom),indsym(4,nsym,dtset%natom)
 254  integer, intent(in) :: nattyp(dtset%ntypat),pertsy(3,dtset%natom+6)
 255  integer, intent(in) :: rfpert(mpert),rf2_dirs_from_rfpert_nl(3,3),symq(4,2,nsym),symrec(3,3,nsym)
 256  integer, intent(out) :: ddkfil(3)
 257  integer, intent(inout) :: blkflg(3,mpert,3,mpert)
 258  integer, intent(out) :: clflg(3,mpert)
 259  real(dp), intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol)
 260  real(dp), intent(in) :: dyew(2,3,dtset%natom,3,dtset%natom)
 261  real(dp), intent(in) :: dyfrlo(3,3,dtset%natom)
 262  real(dp), intent(in) :: dyfrnl(dyfr_cplex,3,3,dtset%natom,1+(dtset%natom-1)*dyfr_nondiag)
 263  real(dp), intent(in) :: dyfrx1(2,3,dtset%natom,3,dtset%natom)
 264  real(dp), intent(in) :: dyfrx2(3,3,dtset%natom),dyvdw(2,3,dtset%natom,3,dtset%natom*usevdw)
 265  real(dp), intent(in) :: eltcore(6,6),elteew(6+3*dtset%natom,6),eltfrhar(6,6)
 266  real(dp), intent(in) :: eltfrkin(6,6),eltfrloc(6+3*dtset%natom,6)
 267  real(dp), intent(in) :: eltfrnl(6+3*dtset%natom,6)
 268  real(dp), intent(in) :: eltfrxc(6+3*dtset%natom,6),eltvdw(6+3*dtset%natom,6*usevdw)
 269  real(dp), intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden)
 270  real(dp), intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol)
 271  real(dp), intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),vxc(nfftf,nspden)
 272  real(dp), intent(in) :: vtrial(nfftf,nspden)
 273  real(dp), intent(inout) :: xred(3,dtset%natom)
 274  real(dp), intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)!vz_i
 275  real(dp), intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i
 276  real(dp), intent(inout) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) !vz_i
 277  real(dp), intent(out) :: eigbrd(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eigbrd)
 278  real(dp), intent(out) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eig2nkq)
 279  type(efmasdeg_type),allocatable,intent(in) :: efmasdeg(:)
 280  type(efmasval_type),allocatable,intent(inout) :: efmasval(:,:)
 281  type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:)
 282  type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:)
 283  type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:)
 284  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 285  type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:)
 286  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 287  real(dp),pointer :: eigen1_pert(:,:,:)
 288  real(dp),intent(out) :: occ_rbz_pert(:),eigen0_pert(:),eigenq_pert(:)
 289  real(dp),pointer :: eigenq_fine(:,:,:)
 290  integer, intent(out) :: nkpt_rbz
 291  type(hdr_type),intent(out) :: hdr0,hdr_fine
 292 
 293 !Local variables-------------------------------
 294 !scalars
 295  integer,parameter :: level=11,response=1,formeig1=1,master=0,fake_unit=-666
 296  integer :: ask_accurate,band_index,bantot,bantot_rbz,bdeigrf,bdtot1_index,nsppol,nspinor,band2tot_index
 297  integer :: bdtot_index,choice,cplex,cplex_rhoij,dim_eig2rf,formeig
 298  integer :: gscase,iband,iblok,icase,icase_eq,idir,idir0,idir1,idir2,idir_eq,idir_dkdk,ierr
 299  integer :: ii,ikpt,ikpt1,jband,initialized,iorder_cprj,ipert,ipert_cnt,ipert_eq,ipert_me,ireadwf0
 300  integer :: iscf_mod,iscf_mod_save,isppol,istr,isym,mcg,mcgq,mcg1,mcprj,mcprjq,mband
 301  integer :: mcgmq,mcg1mq,mpw1_mq !+/-q duplicates
 302  integer :: maxidir,me,mgfftf,mkmem_rbz,mk1mem_rbz,mkqmem_rbz,mpw,mpw1,my_nkpt_rbz
 303  integer :: n3xccc,nband_k,ncpgr,ndir,nkpt_eff,nkpt_max,nline_save,nmatel,npert_io,npert_me,nspden_rhoij
 304  integer :: nstep_save,nsym1,ntypat,nwffile,nylmgr,nylmgr1,old_comm_atom,openexit,option,optorth,optthm,pertcase
 305  integer :: rdwr,rdwrpaw,spaceComm,smdelta,timrev_pert,timrev_kpt,to_compute_this_pert
 306  integer :: unitout,useylmgr,useylmgr1,use_rhoijim,vrsddb,dfpt_scfcv_retcode,optn2
 307 #ifdef HAVE_NETCDF
 308  integer :: ncerr,ncid
 309 #endif
 310  real(dp) :: boxcut,dosdeltae,eberry,ecore,ecut_eff,ecutf,edocc,eei,eeig0,eew,efrhar,efrkin,efrloc
 311  real(dp) :: efrnl,efrx1,efrx2,ehart,ehart01,ehart1,eii,ek,ek0,ek1,ek2,eloc0
 312  real(dp) :: elpsp1,enl,enl0,enl1,entropy,enxc,eovl1,epaw1,evdw,exc1,fsum,gsqcut,maxocc,nelectkq
 313  real(dp) :: residm,tolwfr,tolwfr_save,toldfe_save,toldff_save,tolrff_save,tolvrs_save
 314  real(dp) :: ucvol, eig1_r, eig1_i
 315  real(dp) :: residm_mq !+/-q duplicates
 316  logical,parameter :: paral_pert_inplace=.true.,remove_inv=.false.
 317  logical :: first_entry,found_eq_gkk,t_exist,paral_atom,write_1wfk,init_rhor1
 318  logical :: kramers_deg
 319  character(len=fnlen) :: dscrpt,fiden1i,fiwf1i,fiwf1o,fiwfddk,fnamewff(4),gkkfilnam,fname,filnam
 320  character(len=500) :: message
 321  type(crystal_t) :: crystal,ddb_crystal
 322  type(dataset_type), pointer :: dtset_tmp
 323  type(ebands_t) :: ebands_k,ebands_kq,gkk_ebands
 324  type(ebands_t) :: ebands_kmq !+/-q duplicates
 325  type(eigr2d_t)  :: eigr2d,eigi2d
 326  type(gkk_t)     :: gkk2d
 327  type(hdr_type) :: hdr,hdr_den,hdr_tmp
 328  type(ddb_hdr_type) :: ddb_hdr
 329  type(pawang_type) :: pawang1
 330  type(wffile_type) :: wff1,wffgs,wffkq,wffnow,wfftgs,wfftkq
 331  type(wfk_t) :: ddk_f(4)
 332  type(wvl_data) :: wvl
 333 !arrays
 334  integer :: eq_symop(3,3),ngfftf(18),file_index(4),rfdir(9),rf2dir(9),rf2_dir1(3),rf2_dir2(3)
 335  integer,allocatable :: blkflg_save(:,:,:,:),dimcprj_srt(:),dummy(:),dyn(:),indkpt1(:),indkpt1_tmp(:)
 336  integer,allocatable :: indsy1(:,:,:),irrzon1(:,:,:),istwfk_rbz(:),istwfk_pert(:,:,:)
 337  integer,allocatable :: kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:),npwarr(:),npwtot(:)
 338  integer,allocatable :: kg1_mq(:,:),npwar1_mq(:),npwtot1_mq(:) !+q/-q duplicates
 339  integer,allocatable :: npwtot1(:),npwar1_pert(:,:),npwarr_pert(:,:),npwtot_pert(:,:)
 340  integer,allocatable :: pert_calc(:,:),pert_tmp(:,:)
 341  integer,allocatable :: symaf1(:),symaf1_tmp(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:)
 342  integer, pointer :: old_atmtab(:)
 343  real(dp) :: dielt(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2)
 344  real(dp),allocatable :: buffer1(:,:,:,:,:),cg(:,:),cg1(:,:),cg1_active(:,:),cg0_pert(:,:)
 345  real(dp),allocatable :: cg1_pert(:,:,:,:),cgq(:,:),gh0c1_pert(:,:,:,:)
 346  real(dp),allocatable :: doccde_rbz(:),docckqde(:)
 347  real(dp),allocatable :: gh1c_pert(:,:,:,:),eigen0(:),eigen0_copy(:),eigen1(:),eigen1_mean(:)
 348  real(dp),allocatable :: eigenq(:),gh1c_set(:,:),gh0c1_set(:,:),kpq(:,:)
 349  real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_pert(:),occ_rbz(:),occkq(:),kpt_rbz_pert(:,:)
 350  real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons1(:,:,:),resid(:),rhog1(:,:)
 351  real(dp),allocatable :: rhor1_save(:,:,:)
 352  real(dp),allocatable :: rhor1(:,:),rho1wfg(:,:),rho1wfr(:,:),tnons1(:,:),tnons1_tmp(:,:)
 353  real(dp),allocatable :: rhor1_pq(:,:),rhor1_mq(:,:),rhog1_pq(:,:),rhog1_mq(:,:)          !+q/-q duplicates
 354  real(dp),allocatable :: cg_mq(:,:),cg1_mq(:,:),resid_mq(:)                   !
 355  real(dp),allocatable :: cg1_active_mq(:,:),occk_mq(:)                 !
 356  real(dp),allocatable :: kmq(:,:),kmq_rbz(:,:),gh0c1_set_mq(:,:)        !
 357  real(dp),allocatable :: eigen_mq(:),gh1c_set_mq(:,:),docckde_mq(:),eigen1_mq(:)          !
 358  real(dp),allocatable :: vpsp1(:),work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:)
 359  real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:),zeff(:,:,:)
 360  real(dp),allocatable :: phasecg(:,:),gauss(:,:)
 361  real(dp),allocatable :: gkk(:,:,:,:,:)
 362  type(pawcprj_type),allocatable :: cprj(:,:),cprjq(:,:)
 363  type(paw_ij_type),pointer :: paw_ij_pert(:)
 364  type(paw_an_type),pointer :: paw_an_pert(:)
 365  type(pawfgrtab_type),pointer :: pawfgrtab_pert(:)
 366  type(pawrhoij_type),allocatable :: pawrhoij1(:)
 367  type(pawrhoij_type),pointer :: pawrhoij_pert(:)
 368  type(ddb_type) :: ddb
 369 
 370 ! ***********************************************************************
 371 
 372  DBG_ENTER("COLL")
 373 
 374  _IBM6("In dfpt_looppert")
 375 
 376  call timab(141,1,tsec)
 377 
 378  call status(0,dtfil%filstat,iexit,level,'init          ')
 379 
 380 !Structured debugging if prtvol==-level
 381  if(dtset%prtvol==-level)then
 382    write(message,'(80a,a,a)')  ('=',ii=1,80),ch10,' dfpt_looppert : enter , debug mode '
 383    call wrtout(std_out,message,'COLL')
 384  end if
 385 
 386  dfpt_scfcv_retcode = -1
 387  nsppol = dtset%nsppol; nspinor = dtset%nspinor
 388 
 389  kramers_deg=.true.
 390  if (dtset%tim1rev==0) then
 391    kramers_deg=.false.
 392  end if
 393 
 394 !Obtain dimensional translations in reciprocal space gprimd,
 395 !metrics and unit cell volume, from rprimd. Also output rprimd, gprimd and ucvol
 396  call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd)
 397  call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol)
 398 
 399  call crystal_init(dtset%amu_orig(:,1),crystal,dtset%spgroup,dtset%natom,dtset%npsp,&
 400 & psps%ntypat,dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
 401 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,psps%title,&
 402 & symrel=dtset%symrel,tnons=dtset%tnons,symafm=dtset%symafm)
 403 
 404 !Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90
 405  if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
 406    mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:)
 407    ecutf=dtset%pawecutdg
 408  else
 409    mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:)
 410    ecutf=dtset%ecut
 411  end if
 412  ecut_eff=dtset%ecut*(dtset%dilatmx)**2
 413 
 414 !Compute large sphere cut-off gsqcut
 415  if (psps%usepaw==1) then
 416    call wrtout(std_out,ch10//' FFT (fine) grid used for densities/potentials:','COLL')
 417  end if
 418  call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfftf)
 419 
 420 !Various initializations/allocations
 421  iscf_mod=dtset%iscf
 422  ntypat=psps%ntypat
 423  nkpt_max=50;if (xmpi_paral==1) nkpt_max=-1
 424  paral_atom=(dtset%natom/=my_natom)
 425  cplex=2-timrev !cplex=2 ! DEBUG: impose cplex=2
 426  first_entry=.true.
 427  initialized=0
 428  ecore=zero ; ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero ; eii=zero
 429  clflg(:,:)=0 ! Array on calculated perturbations for eig2rf
 430  if (psps%usepaw==1) then
 431    ABI_ALLOCATE(dimcprj_srt,(dtset%natom))
 432    call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
 433  end if
 434 
 435 !Save values of SCF cycle parameters
 436  iscf_mod_save = iscf_mod
 437  nstep_save = dtset%nstep
 438  nline_save = dtset%nline
 439  tolwfr_save = dtset%tolwfr
 440  toldfe_save = dtset%toldfe
 441  toldff_save = dtset%toldff
 442  tolrff_save = dtset%tolrff
 443  tolvrs_save = dtset%tolvrs
 444 
 445 !This dtset will be used in dfpt_scfcv to force non scf calculations for equivalent perturbations
 446  nullify(dtset_tmp)
 447  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
 448    ABI_DATATYPE_ALLOCATE(dtset_tmp,)
 449    call dtset_copy(dtset_tmp, dtset)
 450  else
 451    dtset_tmp => dtset
 452  end if
 453 
 454 !Generate the 1-dimensional phases
 455  ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
 456  ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
 457  call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred)
 458  if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
 459    call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
 460  else
 461    ph1df(:,:)=ph1d(:,:)
 462  end if
 463 
 464 !!Determine existence of pertubations and of pertubation symmetries
 465 !!Create array with pertubations which have to be calculated
 466 ! ABI_ALLOCATE(pert_tmp,(3*mpert))
 467 ! ipert_cnt=0
 468 ! do ipert=1,mpert
 469 !   do idir=1,3
 470 !     if( rfpert(ipert)==1 .and. dtset%rfdir(idir) == 1 )then
 471 !       if ((pertsy(idir,ipert)==1).or.&
 472 !&       ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.&
 473 !&       ((dtset%prepgkk == 1).and.(ipert <= dtset%natom))  ) then
 474 !         ipert_cnt = ipert_cnt+1;
 475 !         pert_tmp(ipert_cnt) = idir+(ipert-1)*3
 476 !       else
 477 !         write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
 478 !&         ' The perturbation idir=',idir,'  ipert=',ipert,' is',ch10,&
 479 !&         ' symmetric of a previously calculated perturbation.',ch10,&
 480 !&         ' So, its SCF calculation is not needed.',ch10
 481 !         call wrtout(std_out,message,'COLL')
 482 !         call wrtout(ab_out,message,'COLL')
 483 !       end if ! Test of existence of symmetry of perturbation
 484 !     end if ! Test of existence of perturbation
 485 !   end do
 486 ! end do
 487 ! ABI_ALLOCATE(pert_calc,(ipert_cnt))
 488 ! do icase=1,ipert_cnt
 489 !   pert_calc(icase)=pert_tmp(icase)
 490 ! end do
 491 ! ABI_DEALLOCATE(pert_tmp)
 492 
 493 !Initialize rf2dir :
 494  rf2_dir1(1:3)=dtset%rf2_pert1_dir(1:3)
 495  rf2_dir2(1:3)=dtset%rf2_pert2_dir(1:3)
 496  if (sum(rf2_dir1)==3.and.sum(rf2_dir2)==3.and.dtset%prepanl==1) then
 497 !  Diagonal terms :
 498    rf2dir(1) = rf2_dirs_from_rfpert_nl(1,1)
 499    rf2dir(2) = rf2_dirs_from_rfpert_nl(2,2)
 500    rf2dir(3) = rf2_dirs_from_rfpert_nl(3,3)
 501 !  Upper triangular terms :
 502    rf2dir(4) = rf2_dirs_from_rfpert_nl(2,3)
 503    rf2dir(5) = rf2_dirs_from_rfpert_nl(1,3)
 504    rf2dir(6) = rf2_dirs_from_rfpert_nl(1,2)
 505 !  Lower triangular terms :
 506    rf2dir(7) = rf2_dirs_from_rfpert_nl(3,2)
 507    rf2dir(8) = rf2_dirs_from_rfpert_nl(3,1)
 508    rf2dir(9) = rf2_dirs_from_rfpert_nl(2,1)
 509  else
 510 !  Diagonal terms :
 511    rf2dir(1) = rf2_dir1(1)*rf2_dir2(1)
 512    rf2dir(2) = rf2_dir1(2)*rf2_dir2(2)
 513    rf2dir(3) = rf2_dir1(3)*rf2_dir2(3)
 514 !  Upper triangular terms :
 515    rf2dir(4) = rf2_dir1(2)*rf2_dir2(3)
 516    rf2dir(5) = rf2_dir1(1)*rf2_dir2(3)
 517    rf2dir(6) = rf2_dir1(1)*rf2_dir2(2)
 518 !  Lower triangular terms :
 519    rf2dir(7) = rf2_dir1(3)*rf2_dir2(2)
 520    rf2dir(8) = rf2_dir1(3)*rf2_dir2(1)
 521    rf2dir(9) = rf2_dir1(2)*rf2_dir2(1)
 522  end if
 523 
 524 !Determine existence of pertubations and of pertubation symmetries
 525 !Create array with pertubations which have to be calculated
 526  ABI_ALLOCATE(pert_tmp,(3,3*(dtset%natom+6)+18))
 527  ipert_cnt=0
 528  do ipert=1,mpert
 529    if (ipert<dtset%natom+10) then
 530      maxidir = 3
 531      rfdir(1:3) = dtset%rfdir(:)
 532      rfdir(4:9) = 0
 533    else
 534      maxidir = 9
 535      rfdir(1:9) = rf2dir(:)
 536    end if
 537    do idir=1,maxidir
 538      to_compute_this_pert = 0
 539      if(ipert<dtset%natom+10 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then
 540        if ((pertsy(idir,ipert)==1).or.&
 541 &       ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.&
 542 &       ((dtset%prepgkk == 1).and.(ipert <= dtset%natom))  ) then
 543          to_compute_this_pert = 1
 544        else
 545          write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
 546 &         ' The perturbation idir=',idir,'  ipert=',ipert,' is',ch10,&
 547 &         ' symmetric of a previously calculated perturbation.',ch10,&
 548 &         ' So, its SCF calculation is not needed.',ch10
 549          call wrtout(std_out,message,'COLL')
 550          call wrtout(ab_out,message,'COLL')
 551        end if ! Test of existence of symmetry of perturbation
 552      else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then
 553        to_compute_this_pert = 1
 554      else if (ipert==dtset%natom+10 .and. rfpert(ipert)==1 .and. idir <= 6) then
 555        if (idir<=3) then
 556          if (rfdir(idir) == 1) to_compute_this_pert = 1
 557        else
 558          if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1
 559        end if
 560      else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1) then
 561        if (idir <= 3 .and. rfdir(idir) == 1) then
 562          to_compute_this_pert = 1
 563        else if (idir>=4.and.idir<=6) then
 564          if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1
 565        else if (idir>=7.and.idir<=9) then
 566          if (rfdir(idir) == 1 .or. rfdir(idir-3) == 1) to_compute_this_pert = 1
 567        end if
 568      end if
 569      if (to_compute_this_pert /= 0) then
 570        ipert_cnt = ipert_cnt+1;
 571        pert_tmp(1,ipert_cnt) = ipert
 572        pert_tmp(2,ipert_cnt) = idir
 573 !      Store "pertcase" in pert_tmp(3,ipert_cnt)
 574        if (ipert<dtset%natom+10) then
 575          pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3
 576        else
 577          pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3
 578        end if
 579      end if
 580    end do ! idir
 581  end do !ipert
 582 ! do ipert=1,mpert
 583 !   if (ipert<dtset%natom+10) then
 584 !     maxidir = 3
 585 !     rfdir(1:3) = dtset%rfdir(:)
 586 !     rfdir(4:9) = 0
 587 !   else
 588 !     maxidir = 9
 589 !     rfdir(1:9) = rf2dir(:)
 590 !   end if
 591 !   do idir=1,maxidir
 592 !     if( rfpert(ipert)==1 .and. rfdir(idir) == 1 )then
 593 !       to_compute_this_pert = 0
 594 !       if (ipert>=dtset%natom+10) then
 595 !         to_compute_this_pert = 1
 596 !       else if ((pertsy(idir,ipert)==1).or.&
 597 !&         ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.&
 598 !&         ((dtset%prepgkk == 1).and.(ipert <= dtset%natom))  ) then
 599 !         to_compute_this_pert = 1
 600 !       end if
 601 !       if (to_compute_this_pert /= 0) then
 602 !         ipert_cnt = ipert_cnt+1;
 603 !         pert_tmp(1,ipert_cnt) = ipert
 604 !         pert_tmp(2,ipert_cnt) = idir
 605 !!        Store "pertcase" in pert_tmp(3,ipert_cnt)
 606 !         if (ipert<dtset%natom+10) then
 607 !           pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3
 608 !         else
 609 !           pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3
 610 !         end if
 611 !       else
 612 !         write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
 613 !&         ' The perturbation idir=',idir,'  ipert=',ipert,' is',ch10,&
 614 !&         ' symmetric of a previously calculated perturbation.',ch10,&
 615 !&         ' So, its SCF calculation is not needed.',ch10
 616 !         call wrtout(std_out,message,'COLL')
 617 !         call wrtout(ab_out,message,'COLL')
 618 !       end if ! Test of existence of symmetry of perturbation
 619 !     end if ! Test of existence of perturbation
 620 !   end do
 621 ! end do
 622  ABI_ALLOCATE(pert_calc,(3,ipert_cnt))
 623  do icase=1,ipert_cnt
 624    pert_calc(:,icase)=pert_tmp(:,icase)
 625  end do
 626  ABI_DEALLOCATE(pert_tmp)
 627 
 628  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
 629    ABI_ALLOCATE(rhor1_save,(cplex*nfftf,nspden,ipert_cnt))
 630    rhor1_save=zero
 631    ABI_ALLOCATE(blkflg_save,(3,mpert,3,mpert))
 632  end if
 633 
 634 ! Initialize quantities for netcdf print
 635  ABI_ALLOCATE(eigen0_copy,(dtset%mband*nkpt*dtset%nsppol))
 636  eigen0_copy(:)=zero
 637 
 638 ! SP : Retreval of the DDB information and computing of effective charge and
 639 ! dielectric tensor
 640  ABI_ALLOCATE(zeff,(3,3,dtset%natom))
 641  if (dtset%getddb .ne. 0 .or. dtset%irdddb .ne. 0 ) then
 642    filnam = dtfil%filddbsin  !'test_DDB'
 643    ABI_ALLOCATE(dummy,(dtset%natom))
 644    call ddb_from_file(ddb,filnam,1,dtset%natom,0,dummy,ddb_crystal,mpi_enreg%comm_world)
 645 !  Get Dielectric Tensor and Effective Charges
 646 !  (initialized to one_3D and zero if the derivatives are not available in
 647 !  the DDB file)
 648    iblok = ddb_get_dielt_zeff(ddb,ddb_crystal,1,0,0,dielt,zeff)
 649    call crystal_free(ddb_crystal)
 650    call ddb_free(ddb)
 651    ABI_DEALLOCATE(dummy)
 652  end if
 653 
 654 !%%%% Parallelization over perturbations %%%%%
 655 !*Define file output/log file names
 656  npert_io=ipert_cnt;if (dtset%nppert<=1) npert_io=0
 657  call localfilnam(mpi_enreg%comm_pert,mpi_enreg%comm_cell_pert,mpi_enreg%comm_world,dtfil%filnam_ds,'_PRT',npert_io)
 658 !Compute the number of perturbation done by the current cpu
 659  if(mpi_enreg%paral_pert==1) then
 660    npert_me = 0 ; ipert_me = 0
 661    do icase=1,ipert_cnt
 662      if (mpi_enreg%distrb_pert(icase)==mpi_enreg%me_pert) npert_me=npert_me +1
 663    end do
 664  end if
 665 
 666 !*Redefine communicators
 667  call set_pert_comm(mpi_enreg,dtset%nppert)
 668 
 669 !*Redistribute PAW on-site data
 670  nullify(old_atmtab,pawfgrtab_pert,pawrhoij_pert,paw_an_pert,paw_ij_pert)
 671  if (paral_pert_inplace) then
 672    call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
 673 &   paw_an,paw_ij,pawfgrtab,pawrhoij)
 674    pawfgrtab_pert=>pawfgrtab ; pawrhoij_pert=>pawrhoij
 675    paw_an_pert   =>paw_an    ; paw_ij_pert  =>paw_ij
 676 
 677  else
 678    call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
 679 &   paw_an,paw_ij,pawfgrtab,pawrhoij,&
 680 &   paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,&
 681 &   pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert)
 682 
 683  end if
 684 
 685  ! We can handle the time limit in dfpt_scfcv in a robust manner only if we have one perturbation.
 686  if (ipert_cnt > 1 .or. mpi_enreg%paral_pert == 1) call disable_timelimit()
 687 
 688 !Loop on perturbations
 689 !==========================================================================
 690  do icase=1,ipert_cnt
 691    _IBM6("In loop on perts")
 692 
 693 !  %%%% Parallelization over perturbations %%%%%
 694 !  Select the perturbations treated by curent processor
 695    if(mpi_enreg%paral_pert==1) then
 696      if (mpi_enreg%distrb_pert(icase)/=mpi_enreg%me_pert) cycle
 697    end if
 698 
 699 !  *Redefine output/log files
 700    call localwrfile(mpi_enreg%comm_cell,icase,npert_io,mpi_enreg%paral_pert,0)
 701 
 702 !!  Retrieve type and direction of the perturbation
 703 !   if (pert_calc(icase) <= dtset%natom*3) then
 704 !     idir = mod(pert_calc(icase),3)
 705 !     if (idir==0) idir=3
 706 !     ipert=( (pert_calc(icase)-idir) / 3 + 1)
 707 !   else if (pert_calc(icase) <= dtset%natom*3+4) then
 708 !     ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1
 709 !     idir = mod(pert_calc(icase),3)
 710 !     if (idir==0) idir=3
 711 !   else
 712 !     ipert = dtset%natom + ((pert_calc(icase) - 3*(dtset%natom+4) - 1) / 9) + 1
 713 !   end if
 714 !   pertcase=idir+(ipert-1)*3
 715 
 716 !  Retrieve type and direction of the perturbation
 717    ipert = pert_calc(1,icase)
 718    idir = pert_calc(2,icase)
 719    istr=idir
 720    pertcase = pert_calc(3,icase)
 721 
 722 !  Init MPI communicator
 723    spaceComm=mpi_enreg%comm_cell
 724    me=mpi_enreg%me_cell
 725 
 726 !  ===== Describe the perturbation in output/log file
 727    _IBM6("IBM6 before print perturbation")
 728 
 729    write(message, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,&
 730 &   ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:)
 731    call wrtout(std_out,message,'COLL')
 732    call wrtout(ab_out,message,'COLL')
 733    if(ipert>=1 .and. ipert<=dtset%natom)then
 734      write(message, '(a,i4,a,i4)' )' Perturbation : displacement of atom',ipert,'   along direction',idir
 735      call wrtout(std_out,message,'COLL')
 736      call wrtout(ab_out,message,'COLL')
 737      if(iscf_mod == -3)then
 738        write(message, '(a,a,a,a,a,a,a,a)' )ch10,&
 739 &       ' dfpt_looppert : COMMENT -',ch10,&
 740 &       '  The first-order density is imposed to be zero (iscf=-3).',ch10,&
 741 &       '  Although this is strange in the case of phonons,',ch10,&
 742 &       '  you are allowed to do so.'
 743        call wrtout(std_out,message,'COLL')
 744        call wrtout(ab_out,message,'COLL')
 745      end if
 746    else if(ipert==dtset%natom+1)then
 747      write(message,'(a,i4)')' Perturbation : derivative vs k along direction',idir
 748      call wrtout(std_out,message,'COLL')
 749      call wrtout(ab_out,message,'COLL')
 750      if( iscf_mod /= -3 )then
 751        write(message, '(4a)' )ch10,&
 752 &       ' dfpt_looppert : COMMENT -',ch10,&
 753 &       '  In a d/dk calculation, iscf is set to -3 automatically.'
 754        call wrtout(std_out,message,'COLL')
 755        call wrtout(ab_out,message,'COLL')
 756        iscf_mod=-3
 757      end if
 758      if( abs(dtset%dfpt_sciss) > 1.0d-8 )then
 759        write(message, '(a,a,a,a,f14.8,a,a)' )ch10,&
 760 &       ' dfpt_looppert : WARNING -',ch10,&
 761 &       '  Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,&
 762 &       '  Scissor with d/dk calculation : you are using a "naive" approach !'
 763        call wrtout(std_out,message,'COLL')
 764        call wrtout(ab_out,message,'COLL')
 765      end if
 766    else if(ipert==dtset%natom+2)then
 767      write(message, '(a,i4)' )' Perturbation : homogeneous electric field along direction',idir
 768      call wrtout(std_out,message,'COLL')
 769      call wrtout(ab_out,message,'COLL')
 770      if( iscf_mod == -3 )then
 771        write(message, '(a,a,a,a,a,a)' )ch10,&
 772 &       ' dfpt_looppert : COMMENT -',ch10,&
 773 &       '  The first-order density is imposed to be zero (iscf=-3).',ch10,&
 774 &       '  This corresponds to a calculation without local fields.'
 775        call wrtout(std_out,message,'COLL')
 776        call wrtout(ab_out,message,'COLL')
 777      end if
 778    else if(ipert==dtset%natom+10.or.ipert==dtset%natom+11)then
 779      call rf2_getidirs(idir,idir1,idir2)
 780      if(ipert==dtset%natom+10)then
 781        write(message,'(2(a,i1))') ' Perturbation : 2nd derivative wrt k, idir1 = ',idir1,&
 782 &       ' idir2 = ',idir2
 783      else
 784        write(message,'(2(a,i1),a)') ' Perturbation : 2nd derivative wrt k (idir1 =',idir1,&
 785 &       ') and Efield (idir2 =',idir2,')'
 786      end if
 787      call wrtout(std_out,message,'COLL')
 788      call wrtout(ab_out,message,'COLL')
 789      if( iscf_mod /= -3 )then
 790        write(message, '(4a)' )ch10,&
 791 &       ' dfpt_looppert : COMMENT -',ch10,&
 792 &       '  In this case, iscf is set to -3 automatically.'
 793        call wrtout(std_out,message,'COLL')
 794        call wrtout(ab_out,message,'COLL')
 795        iscf_mod=-3
 796      end if
 797      if( abs(dtset%dfpt_sciss) > 1.0d-8 )then
 798        write(message, '(a,a,a,a,f14.8,a,a)' )ch10,&
 799 &       ' dfpt_looppert : WARNING -',ch10,&
 800 &       '  Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,&
 801 &       '  Scissor with d/dk calculation : you are using a "naive" approach !'
 802        call wrtout(std_out,message,'COLL')
 803        call wrtout(ab_out,message,'COLL')
 804      end if
 805      ABI_ALLOCATE(occ_pert,(dtset%mband*nkpt*dtset%nsppol))
 806      occ_pert(:) = occ(:) - occ(1)
 807      maxocc = maxval(abs(occ_pert))
 808      if (maxocc>1.0d-6.and.abs(maxocc-occ(1))>1.0d-6) then ! True if non-zero occupation numbers are not equal
 809        write(message, '(3a)' ) ' ipert=natom+10 or 11 does not work for a metallic system.',ch10,&
 810        ' This perturbation will not be computed.'
 811        MSG_WARNING(message)
 812        ABI_DEALLOCATE(occ_pert)
 813        cycle
 814      end if
 815      ABI_DEALLOCATE(occ_pert)
 816    else if(ipert>dtset%natom+11 .or. ipert<=0 )then
 817      write(message, '(a,i4,a,a,a)' ) &
 818 &     '  ipert=',ipert,' is outside the [1,dtset%natom+11] interval.',ch10,&
 819 &     '  This perturbation is not (yet) allowed.'
 820      MSG_BUG(message)
 821    end if
 822 !  Initialize the diverse parts of energy :
 823    eew=zero ; evdw=zero ; efrloc=zero ; efrnl=zero ; efrx1=zero ; efrx2=zero
 824    efrhar=zero ; efrkin=zero
 825    if(ipert<=dtset%natom)then
 826      eew=dyew(1,idir,ipert,idir,ipert)
 827      if (usevdw==1) evdw=dyvdw(1,idir,ipert,idir,ipert)
 828      efrloc=dyfrlo(idir,idir,ipert)
 829      if (dyfr_nondiag==0) efrnl=dyfrnl(1,idir,idir,ipert,1)
 830      if (dyfr_nondiag/=0) efrnl=dyfrnl(1,idir,idir,ipert,ipert)
 831      efrx1=dyfrx1(1,idir,ipert,idir,ipert)
 832      efrx2=dyfrx2(idir,idir,ipert)
 833    else if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
 834 !    istr = 1,2,...,6 and indicates the cartesian strain component
 835      if(ipert==dtset%natom+4) istr=idir+3
 836      eii=eltcore(istr,istr)
 837      eew=elteew(istr,istr)
 838      if (usevdw==1) evdw=eltvdw(istr,istr)
 839      efrhar=eltfrhar(istr,istr)
 840      efrkin=eltfrkin(istr,istr)
 841      efrloc=eltfrloc(istr,istr)
 842      efrnl=eltfrnl(istr,istr)
 843      efrx1=eltfrxc(istr,istr)
 844    end if
 845 
 846 !  Determine the subset of symmetry operations (nsym1 operations)
 847 !  that leaves the perturbation invariant, and initialize corresponding arrays
 848 !  symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW)..
 849    ABI_ALLOCATE(symaf1_tmp,(nsym))
 850    ABI_ALLOCATE(symrl1_tmp,(3,3,nsym))
 851    ABI_ALLOCATE(tnons1_tmp,(3,nsym))
 852    if (dtset%prepanl/=1.and.&
 853 &   dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
 854 &   dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17) then
 855      call littlegroup_pert(gprimd,idir,indsym,ab_out,ipert,dtset%natom,nsym,nsym1,2,&
 856 &     dtset%symafm,symaf1_tmp,symq,symrec,dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp)
 857    else
 858      nsym1 = 1
 859      symaf1_tmp(1) = 1
 860      symrl1_tmp(:,:,1) = dtset%symrel(:,:,1)
 861      tnons1_tmp(:,1) = 0_dp
 862    end if
 863    ABI_ALLOCATE(indsy1,(4,nsym1,dtset%natom))
 864    ABI_ALLOCATE(symrc1,(3,3,nsym1))
 865    ABI_ALLOCATE(symaf1,(nsym1))
 866    ABI_ALLOCATE(symrl1,(3,3,nsym1))
 867    ABI_ALLOCATE(tnons1,(3,nsym1))
 868    symaf1(1:nsym1)=symaf1_tmp(1:nsym1)
 869    symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1)
 870    tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1)
 871    ABI_DEALLOCATE(symaf1_tmp)
 872    ABI_DEALLOCATE(symrl1_tmp)
 873    ABI_DEALLOCATE(tnons1_tmp)
 874 
 875 !  Set up corresponding symmetry data
 876    ABI_ALLOCATE(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4)))
 877    ABI_ALLOCATE(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4)))
 878    call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,&
 879 &   nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)
 880    if (psps%usepaw==1) then
 881 !    Allocate/initialize only zarot in pawang1 datastructure
 882      call pawang_init(pawang1,0,pawang%l_max-1,0,nsym1,0,1,0,0,0)
 883      call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot)
 884    end if
 885 
 886 !  Initialize k+q array
 887    ABI_ALLOCATE(kpq,(3,nkpt))
 888    if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
 889      kpq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran
 890    else
 891      do ikpt=1,nkpt
 892        kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt)
 893      end do
 894    end if
 895 !  In case wf1 at +q and -q are not related by time inversion symmetry, compute k-q as well
 896 !  for initializations
 897    if (.not.kramers_deg) then
 898      ABI_ALLOCATE(kmq,(3,nkpt))
 899      if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
 900        kmq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran
 901      else
 902        do ikpt=1,nkpt
 903          kmq(:,ikpt)=-dtset%qptn(:)+dtset%kptns(:,ikpt) ! kmq <= k-q
 904        end do
 905      end if
 906    end if
 907 
 908 !  Determine the subset of k-points needed in the "reduced Brillouin zone",
 909 !  and initialize other quantities
 910    ABI_ALLOCATE(indkpt1_tmp,(nkpt))
 911    ABI_ALLOCATE(wtk_folded,(nkpt))
 912    indkpt1_tmp(:)=0 ; optthm=0
 913    timrev_pert=timrev
 914    if(dtset%ieig2rf>0) then
 915      timrev_pert=0
 916      call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,&
 917 &     1,symrc1,timrev_pert,dtset%wtk,wtk_folded)
 918    else
 919 !    For the time being, the time reversal symmetry is not used
 920 !    for ddk, elfd, mgfd perturbations.
 921      timrev_pert=timrev
 922      if(ipert==dtset%natom+1.or.ipert==dtset%natom+2.or.&
 923 &     ipert==dtset%natom+10.or.ipert==dtset%natom+11.or. &
 924 &     dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.  &
 925 &     dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17.or.  &
 926 &     ipert==dtset%natom+5.or.dtset%prtfull1wf==1) timrev_pert=0
 927      timrev_kpt = timrev_pert
 928 !    The time reversal symmetry is not used for the BZ sampling when kptopt=3 or 4
 929      if (dtset%kptopt==3.or.dtset%kptopt==4) timrev_kpt = 0
 930      call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,&
 931      nsym1,symrc1,timrev_kpt,dtset%wtk,wtk_folded)
 932    end if
 933 
 934    ABI_ALLOCATE(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
 935    ABI_ALLOCATE(indkpt1,(nkpt_rbz))
 936    ABI_ALLOCATE(istwfk_rbz,(nkpt_rbz))
 937    ABI_ALLOCATE(kpq_rbz,(3,nkpt_rbz))
 938    if (.not.kramers_deg) then
 939      ABI_ALLOCATE(kmq_rbz,(3,nkpt_rbz))
 940    end if
 941    ABI_ALLOCATE(kpt_rbz,(3,nkpt_rbz))
 942    ABI_ALLOCATE(nband_rbz,(nkpt_rbz*dtset%nsppol))
 943    ABI_ALLOCATE(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
 944    ABI_ALLOCATE(wtk_rbz,(nkpt_rbz))
 945    indkpt1(:)=indkpt1_tmp(1:nkpt_rbz)
 946    do ikpt=1,nkpt_rbz
 947      istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt))
 948      kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt))
 949      kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt))
 950      wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt))
 951    end do
 952    if (.not.kramers_deg) then
 953      do ikpt=1,nkpt_rbz
 954        kmq_rbz(:,ikpt)=kmq(:,indkpt1(ikpt))
 955      end do
 956    end if
 957    ABI_DEALLOCATE(indkpt1_tmp)
 958    ABI_DEALLOCATE(wtk_folded)
 959 
 960 !  Transfer occ to occ_rbz and doccde to doccde_rbz :
 961 !  this is a more delicate issue
 962 !  NOTE : this takes into account that indkpt1 is ordered
 963 !  MG: What about using occ(band,kpt,spin) ???
 964    bdtot_index=0;bdtot1_index=0
 965    do isppol=1,dtset%nsppol
 966      ikpt1=1
 967      do ikpt=1,nkpt
 968        nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
 969 !      Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1)
 970        if(ikpt1/=nkpt_rbz+1)then
 971          if(ikpt==indkpt1(ikpt1))then
 972            nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k
 973            occ_rbz(1+bdtot1_index:nband_k+bdtot1_index)    = occ(1+bdtot_index:nband_k+bdtot_index)
 974            doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index) = doccde(1+bdtot_index:nband_k+bdtot_index)
 975            ikpt1=ikpt1+1
 976            bdtot1_index=bdtot1_index+nband_k
 977          end if
 978        end if
 979        bdtot_index=bdtot_index+nband_k
 980      end do
 981    end do
 982 
 983    _IBM6("IBM6 in dfpt_looppert before getmpw")
 984 
 985 !  Compute maximum number of planewaves at k
 986    call timab(142,1,tsec)
 987    call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz)
 988    call timab(142,2,tsec)
 989 
 990 !  Allocate some k-dependent arrays at k
 991    ABI_ALLOCATE(kg,(3,mpw*nkpt_rbz))
 992    ABI_ALLOCATE(npwarr,(nkpt_rbz))
 993    ABI_ALLOCATE(npwtot,(nkpt_rbz))
 994 
 995 !  Determine distribution of k-points/bands over MPI processes
 996    if (allocated(mpi_enreg%my_kpttab)) then
 997      ABI_DEALLOCATE(mpi_enreg%my_kpttab)
 998    end if
 999    ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz))
1000    if(xmpi_paral==1) then
1001      ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol))
1002      call distrb2(dtset%mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg)
1003    else
1004      mpi_enreg%my_kpttab(:)=(/(ii,ii=1,nkpt_rbz)/)
1005    end if
1006    my_nkpt_rbz=maxval(mpi_enreg%my_kpttab)
1007    call initmpi_band(mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol)
1008    mkmem_rbz =my_nkpt_rbz ; mkqmem_rbz=my_nkpt_rbz ; mk1mem_rbz=my_nkpt_rbz
1009    ABI_UNUSED((/mkmem,mk1mem,mkqmem/))
1010 
1011    _IBM6("IBM6 before kpgio")
1012 
1013 !  Set up the basis sphere of planewaves at k
1014    call timab(143,1,tsec)
1015    call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,&
1016 &   kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw,npwarr,npwtot,dtset%nsppol)
1017    call timab(143,2,tsec)
1018 
1019 !  Set up the spherical harmonics (Ylm) at k
1020    useylmgr=0; option=0 ; nylmgr=0
1021    if (psps%useylm==1.and. &
1022 &   (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. &
1023 &   (psps%usepaw==1.and.ipert==dtset%natom+2))) then
1024      useylmgr=1; option=1 ; nylmgr=3
1025    else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then
1026      useylmgr=1; option=2 ; nylmgr=9
1027    end if
1028    ABI_ALLOCATE(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
1029    ABI_ALLOCATE(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
1030    if (psps%useylm==1) then
1031      call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,&
1032 &     npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
1033    end if
1034 
1035    _IBM6("Before ieig2rf > 0")
1036 
1037 !  Set up occupations for this perturbation
1038    if (dtset%ieig2rf>0) then
1039      if (.not.allocated(istwfk_pert)) then
1040        ABI_ALLOCATE(istwfk_pert,(nkpt,3,mpert))
1041        ABI_ALLOCATE(occ_pert,(dtset%mband*nkpt*dtset%nsppol))
1042        istwfk_pert(:,:,:)=0 ; occ_pert(:)=zero
1043      end if
1044      istwfk_pert(:,idir,ipert)=istwfk_rbz(:)
1045      occ_pert(:)=occ_rbz(:)
1046    end if
1047    if (dtset%efmas>0) then
1048      if (.not.allocated(istwfk_pert)) then
1049        ABI_ALLOCATE(istwfk_pert,(nkpt,3,mpert))
1050        istwfk_pert(:,:,:)=0
1051      end if
1052      istwfk_pert(:,idir,ipert)=istwfk_rbz(:)
1053    end if
1054 
1055 !  Print a separator in output file
1056    write(message,'(3a)')ch10,'--------------------------------------------------------------------------------',ch10
1057    call wrtout(ab_out,message,'COLL')
1058 
1059 !  Initialize band structure datatype at k
1060    bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol))
1061    ABI_ALLOCATE(eigen0,(bantot_rbz))
1062    eigen0(:)=zero
1063    call ebands_init(bantot_rbz,ebands_k,dtset%nelect,doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,&
1064 &   nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,&
1065 &   dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1066 &   dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1067    ABI_DEALLOCATE(eigen0)
1068 
1069 !  Initialize header, update it with evolving variables
1070    gscase=0 ! A GS WF file is read
1071    call hdr_init(ebands_k,codvsn,dtset,hdr0,pawtab,gscase,psps,wvl%descr,&
1072 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1073 
1074    call hdr_update(hdr0,bantot_rbz,etotal,fermie,&
1075 &   residm,rprimd,occ_rbz,pawrhoij_pert,xred,dtset%amu_orig(:,1),&
1076 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1077 
1078    _IBM6("before inwffil")
1079 
1080 !  Initialize GS wavefunctions at k
1081    ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0
1082    mcg=mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol
1083    if (one*mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol > huge(1)) then
1084      write (message,'(4a, 5(a,i0), 2a)')&
1085 &     "Default integer is not wide enough to store the size of the GS wavefunction array (WF0, mcg).",ch10,&
1086 &     "Action: increase the number of processors. Consider also OpenMP threads.",ch10,&
1087 &     "nspinor: ",dtset%nspinor, "mpw: ",mpw, "mband: ",dtset%mband, "mkmem_rbz: ",&
1088 &     mkmem_rbz, "nsppol: ",dtset%nsppol,ch10,&
1089 &     'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode'
1090      MSG_ERROR(message)
1091    end if
1092    ABI_STAT_ALLOCATE(cg,(2,mcg), ierr)
1093    ABI_CHECK(ierr==0, "out-of-memory in cg")
1094 
1095    ABI_ALLOCATE(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol))
1096    call timab(144,1,tsec)
1097    call status(0,dtfil%filstat,iexit,level,'call inwffil-k')
1098    call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,&
1099 &   formeig,hdr0,ireadwf0,istwfk_rbz,kg,&
1100 &   kpt_rbz,dtset%localrdwf,dtset%mband,mcg,&
1101 &   mkmem_rbz,mpi_enreg,mpw,nband_rbz,dtset%ngfft,nkpt_rbz,npwarr,&
1102 &   dtset%nsppol,nsym,occ_rbz,optorth,dtset%symafm,&
1103 &   dtset%symrel,dtset%tnons,dtfil%unkg,wffgs,wfftgs,&
1104 &   dtfil%unwffgs,dtfil%fnamewffk,wvl)
1105    call timab(144,2,tsec)
1106 !  Close wffgs%unwff, if it was ever opened (in inwffil)
1107    if (ireadwf0==1) then
1108      call WffClose(wffgs,ierr)
1109    end if
1110    ! Update energies GS energies at k
1111    call put_eneocc_vect(ebands_k, "eig", eigen0)
1112 
1113 !  PAW: compute on-site projections of GS wavefunctions (cprj) (and derivatives) at k
1114    ncpgr=0
1115    ABI_DATATYPE_ALLOCATE(cprj,(0,0))
1116    if (psps%usepaw==1) then
1117      ncpgr=3 ! Valid for ipert<=natom (phonons), ipert=natom+2 (elec. field)
1118              ! or for ipert==natom+10,11
1119      if (ipert==dtset%natom+1) ncpgr=1
1120      if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) ncpgr=1
1121      if (usecprj==1) then
1122        mcprj=dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol
1123        ABI_DATATYPE_DEALLOCATE(cprj)
1124        ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj))
1125        call pawcprj_alloc(cprj,ncpgr,dimcprj_srt)
1126        if (ipert<=dtset%natom) then
1127          choice=2; iorder_cprj=0; idir0=0
1128        else if (ipert==dtset%natom+1) then
1129          choice=5; iorder_cprj=0; idir0=idir
1130        else if (ipert==dtset%natom+2) then
1131          choice=5; iorder_cprj=0; idir0=0
1132        else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
1133          choice=3; iorder_cprj=0; idir0=istr
1134        else if (ipert==dtset%natom+10.or.ipert==dtset%natom+11) then
1135          choice=5; iorder_cprj=0; idir0=0 ! Compute all first derivatives
1136        else
1137          choice=1; iorder_cprj=0; idir0=0
1138        end if
1139        call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,-1,idir0,iorder_cprj,istwfk_rbz,&
1140 &       kg,kpt_rbz,mcg,mcprj,dtset%mgfft,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,&
1141 &       dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,&
1142 &       npwarr,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,&
1143 &       rmet,dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr)
1144      end if
1145    end if
1146 
1147 !  Compute maximum number of planewaves at k+q
1148 !  Will be useful for both GS wfs at k+q and RF wavefunctions
1149    call timab(143,1,tsec)
1150    call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpq_rbz,mpi_enreg,mpw1,nkpt_rbz)
1151    if (.not.kramers_deg) then
1152      call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kmq_rbz,mpi_enreg,mpw1_mq,nkpt_rbz)
1153      !number of plane waves at k+q and k-q should be in principle the same to reconstruct rhor1_pq (?)
1154      !mpw1=max(mpw1,mpw1_tmp)
1155    else
1156      mpw1_mq=0
1157    end if
1158    call timab(143,2,tsec)
1159 
1160 !  Allocate some arrays at k+q
1161    ABI_ALLOCATE(kg1,(3,mpw1*mk1mem_rbz))
1162    ABI_ALLOCATE(npwar1,(nkpt_rbz))
1163    ABI_ALLOCATE(npwtot1,(nkpt_rbz))
1164 !  In case Kramers degeneracy is broken, do the same for k-q
1165    if (.not.kramers_deg) then
1166      ABI_ALLOCATE(kg1_mq,(3,mpw1_mq*mk1mem_rbz))
1167      ABI_ALLOCATE(npwar1_mq,(nkpt_rbz))
1168      ABI_ALLOCATE(npwtot1_mq,(nkpt_rbz))
1169    end if
1170 
1171 !  Set up the basis sphere of planewaves at k+q
1172 !  Will be useful for both GS wfs at k+q and RF wavefunctions
1173    call timab(142,1,tsec)
1174    call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1,&
1175 &   kpq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1,&
1176 &   npwar1,npwtot1,dtset%nsppol)
1177    if (.not.kramers_deg) then
1178      call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1_mq,&
1179 &     kmq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1_mq,&
1180 &     npwar1_mq,npwtot1_mq,dtset%nsppol)
1181    end if
1182    call timab(142,2,tsec)
1183 
1184 !  Set up the spherical harmonics (Ylm) at k+q
1185    useylmgr1=0; option=0 ; ; nylmgr1=0
1186    if (psps%useylm==1.and. &
1187 &   (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. &
1188 &   (psps%usepaw==1.and.ipert==dtset%natom+2))) then
1189      useylmgr1=1; option=1 ;  ; nylmgr1=3
1190    else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then
1191      useylmgr1=1; option=2 ;  ; nylmgr1=9
1192    end if
1193    ABI_ALLOCATE(ylm1,(mpw1*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
1194    ABI_ALLOCATE(ylmgr1,(mpw1*mk1mem_rbz,nylmgr1,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
1195    if (psps%useylm==1) then
1196      call initylmg(gprimd,kg1,kpq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1,nband_rbz,nkpt_rbz,&
1197 &     npwar1,dtset%nsppol,option,rprimd,ylm1,ylmgr1)
1198    end if
1199 !  SPr: do the same for k-q if kramers_deg=.false.
1200 
1201 !  Print a separator in output file
1202    write(message, '(a,a)' )'--------------------------------------------------------------------------------',ch10
1203    call wrtout(ab_out,message,'COLL')
1204 
1205 !  Initialize band structure datatype at k+q
1206    ABI_ALLOCATE(eigenq,(bantot_rbz))
1207    eigenq(:)=zero
1208    call ebands_init(bantot_rbz,ebands_kq,dtset%nelect,doccde_rbz,eigenq,istwfk_rbz,kpq_rbz,&
1209 &   nband_rbz,nkpt_rbz,npwar1,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,&
1210 &   dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1211 &   dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1212    if (.not.kramers_deg) then
1213      eigenq(:)=zero
1214      call ebands_init(bantot_rbz,ebands_kmq,dtset%nelect,doccde_rbz,eigenq,istwfk_rbz,kmq_rbz,&
1215 &     nband_rbz,nkpt_rbz,npwar1_mq,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,&
1216 &     dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1217 &     dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1218    end if
1219    ABI_DEALLOCATE(eigenq)
1220 
1221 !  Initialize header
1222    call hdr_init(ebands_kq,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, &
1223 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
1224    if (.not.kramers_deg) then
1225      call hdr_init(ebands_kmq,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, &
1226 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
1227    end if
1228 
1229 !  Initialize wavefunctions at k+q
1230 !  MG: Here it is possible to avoid the extra reading if the same k mesh can be used.
1231    ireadwf0=1 ; formeig=0 ; ask_accurate=1 ; optorth=0
1232    mcgq=mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol
1233    !SPr: verified until here, add mcgq for -q
1234    if (one*mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol > huge(1)) then
1235      write (message,'(4a, 5(a,i0), 2a)')&
1236 &     "Default integer is not wide enough to store the size of the GS wavefunction array (WFKQ, mcgq).",ch10,&
1237 &     "Action: increase the number of processors. Consider also OpenMP threads.",ch10,&
1238 &     "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mkqmem_rbz: ",&
1239 &     mkqmem_rbz, "nsppol: ",dtset%nsppol,ch10,&
1240 &     'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode'
1241      MSG_ERROR(message)
1242    end if
1243 
1244    ABI_STAT_ALLOCATE(cgq,(2,mcgq), ierr)
1245    ABI_CHECK(ierr==0, "out-of-memory in cgq")
1246    ABI_ALLOCATE(eigenq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1247    if (.not.kramers_deg) then
1248      !ABI_STAT_ALLOCATE(cg_pq,(2,mcgq), ierr)
1249      !ABI_CHECK(ierr==0, "out-of-memory in cgmq")
1250      !ABI_ALLOCATE(eigen_pq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1251      mcgmq=mpw1_mq*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol
1252      ABI_STAT_ALLOCATE(cg_mq,(2,mcgmq), ierr)
1253      ABI_CHECK(ierr==0, "out-of-memory in cgmq")
1254      ABI_ALLOCATE(eigen_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1255    end if
1256 
1257    !if (sum(dtset%qptn(1:3)**2)>=1.d-14) then ! non-zero q
1258    if (dtfil%fnamewffq == dtfil%fnamewffk .and. sum(dtset%qptn(1:3)**2) < 1.d-14) then
1259      call wrtout(std_out, "qpt is Gamma, psi_k+q initialized from psi_k in memory")
1260      cgq = cg
1261      eigenq = eigen0
1262    else
1263      call timab(144,1,tsec)
1264      call status(0,dtfil%filstat,iexit,level,'call inwffilkq')
1265      call inwffil(ask_accurate,cgq,dtset,dtset%ecut,ecut_eff,eigenq,dtset%exchn2n3d,&
1266 &     formeig,hdr,&
1267 &     ireadwf0,istwfk_rbz,kg1,kpq_rbz,dtset%localrdwf,dtset%mband,mcgq,&
1268 &     mkqmem_rbz,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,&
1269 &     dtset%nsppol,nsym,occ_rbz,optorth,&
1270 &     dtset%symafm,dtset%symrel,dtset%tnons,&
1271 &     dtfil%unkg1,wffkq,wfftkq,dtfil%unwffkq,dtfil%fnamewffq,wvl)
1272      call timab(144,2,tsec)
1273 !    Close dtfil%unwffkq, if it was ever opened (in inwffil)
1274      if (ireadwf0==1) then
1275        call WffClose(wffkq,ierr)
1276      end if
1277      if (.not.kramers_deg) then
1278        !SPr: later "make" a separate WFQ file for "-q"
1279        call timab(144,1,tsec)
1280        call status(0,dtfil%filstat,iexit,level,'call inwffilkq')
1281        call inwffil(ask_accurate,cg_mq,dtset,dtset%ecut,ecut_eff,eigen_mq,dtset%exchn2n3d,&
1282 &       formeig,hdr,&
1283 &       ireadwf0,istwfk_rbz,kg1_mq,kmq_rbz,dtset%localrdwf,dtset%mband,mcgmq,&
1284 &       mkqmem_rbz,mpi_enreg,mpw1_mq,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1_mq,&
1285 &       dtset%nsppol,nsym,occ_rbz,optorth,&
1286 &       dtset%symafm,dtset%symrel,dtset%tnons,&
1287 &       dtfil%unkg1,wffkq,wfftkq,dtfil%unwffkq,dtfil%fnamewffq,wvl)
1288        call timab(144,2,tsec)
1289 !      Close dtfil%unwffkq, if it was ever opened (in inwffil)
1290        if (ireadwf0==1) then
1291          call WffClose(wffkq,ierr)
1292        end if
1293      end if
1294    end if
1295    ! Update energies GS energies at k + q
1296    call put_eneocc_vect(ebands_kq, "eig", eigenq)
1297    if (.not.kramers_deg) then
1298      call put_eneocc_vect(ebands_kmq, "eig", eigen_mq)
1299    end if
1300 
1301 !  PAW: compute on-site projections of GS wavefunctions (cprjq) (and derivatives) at k+q
1302    ABI_DATATYPE_ALLOCATE(cprjq,(0,0))
1303    if (psps%usepaw==1) then
1304      if (usecprj==1) then
1305        call status(0,dtfil%filstat,iexit,level,'call ctocprjq')
1306        mcprjq=dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol
1307        ABI_DATATYPE_DEALLOCATE(cprjq)
1308        ABI_DATATYPE_ALLOCATE(cprjq,(dtset%natom,mcprjq))
1309        call pawcprj_alloc(cprjq,0,dimcprj_srt)
1310        if (ipert<=dtset%natom.and.(sum(dtset%qptn(1:3)**2)>=1.d-14)) then ! phonons at non-zero q
1311          choice=1 ; iorder_cprj=0 ; idir0=0
1312          call ctocprj(atindx,cgq,choice,cprjq,gmet,gprimd,-1,idir0,0,istwfk_rbz,&
1313 &         kg1,kpq_rbz,mcgq,mcprjq,dtset%mgfft,mkqmem_rbz,mpi_enreg,psps%mpsang,mpw1,&
1314 &         dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,&
1315 &         npwar1,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,&
1316 &         psps,rmet,dtset%typat,ucvol,dtfil%unpawq,xred,ylm1,ylmgr1)
1317        else if (mcprjq>0) then
1318          call pawcprj_copy(cprj,cprjq)
1319        end if
1320      end if
1321    end if
1322 
1323 !  ===== Report on eigenq values
1324    if (dtset%ieig2rf>0.and.icase==ipert_cnt) then
1325      eigen0_pert(:) = eigen0(:)
1326      eigenq_pert(:) = eigenq(:)
1327      occ_rbz_pert(:) = occ_rbz(:)
1328    end if
1329    if (dtset%efmas>0.and.icase==ipert_cnt) then
1330      eigen0_pert(:) = eigen0(:)
1331    end if
1332    call wrtout(std_out,ch10//' dfpt_looppert : eigenq array',"COLL")
1333    nkpt_eff=nkpt
1334    if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. nkpt>nkpt_max ) nkpt_eff=nkpt_max
1335    band_index=0
1336    do isppol=1,dtset%nsppol
1337      do ikpt=1,nkpt_rbz
1338        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
1339        if(ikpt<=nkpt_eff)then
1340          write(message, '(a,i2,a,i5)' )'  isppol=',isppol,', k point number',ikpt
1341          call wrtout(std_out,message,'COLL')
1342          do iband=1,nband_k,4
1343            write(message, '(a,4es16.6)')'  ',eigenq(iband+band_index:min(iband+3,nband_k)+band_index)
1344            call wrtout(std_out,message,'COLL')
1345          end do
1346        else if(ikpt==nkpt_eff+1)then
1347          write(message,'(a,a)' )'  respfn : prtvol=0, 1 or 2, stop printing eigenq.',ch10
1348          call wrtout(std_out,message,'COLL')
1349        end if
1350        band_index=band_index+nband_k
1351      end do
1352    end do
1353 
1354 !  Generate occupation numbers for the reduced BZ at k+q
1355    ABI_ALLOCATE(docckqde,(dtset%mband*nkpt_rbz*dtset%nsppol))
1356    ABI_ALLOCATE(occkq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1357    if (.not.kramers_deg) then
1358      ABI_ALLOCATE(docckde_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1359      ABI_ALLOCATE(occk_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1360    end if
1361 
1362    if(0<=dtset%occopt .and. dtset%occopt<=2)then
1363 !    Same occupation numbers at k and k+q (usually, insulating)
1364      occkq(:)=occ_rbz(:)
1365      docckqde(:)=zero  ! docckqde is irrelevant in this case
1366      if(.not.kramers_deg) then
1367        occk_mq(:)=occ_rbz(:)
1368        docckde_mq(:)=zero
1369      end if
1370    else
1371 !    Metallic occupation numbers
1372      option=1
1373      dosdeltae=zero ! the DOS is not computed with option=1
1374      maxocc=two/(dtset%nspinor*dtset%nsppol)
1375      call getnel(docckqde,dosdeltae,eigenq,entropy,fermie,maxocc,dtset%mband,&
1376 &     nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occkq,dtset%occopt,option,&
1377 &     dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz)
1378 !    Compare nelect at k and nelelect at k+q
1379      write(message, '(a,a,a,es16.6,a,es16.6,a)')&
1380 &     ' dfpt_looppert : total number of electrons, from k and k+q',ch10,&
1381 &     '  fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.'
1382      call wrtout(ab_out,message,'COLL')
1383      call wrtout(std_out,message,'COLL')
1384      if (.not.kramers_deg) then
1385        call getnel(docckde_mq,dosdeltae,eigen_mq,entropy,fermie,maxocc,dtset%mband,&
1386 &       nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occk_mq,dtset%occopt,option,&
1387 &       dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz)
1388 !      Compare nelect at k and nelelect at k-q
1389        write(message, '(a,a,a,es16.6,a,es16.6,a)')&
1390 &       ' dfpt_looppert : total number of electrons, from k and k-q',ch10,&
1391 &       '  fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.'
1392        call wrtout(ab_out,message,'COLL')
1393        call wrtout(std_out,message,'COLL')
1394      end if
1395    end if
1396 
1397 !  Debug message
1398    if(dtset%prtvol==-level)then
1399      call wrtout(std_out,'dfpt_looppert : initialisation of q part done. ','COLL')
1400    end if
1401 
1402 !  Initialisation of first-order wavefunctions
1403    write(message,'(3a,i4)')&
1404 &   ' Initialisation of the first-order wave-functions :',ch10,&
1405 &   '  ireadwf=',dtfil%ireadwf
1406    call wrtout(std_out,message,'COLL')
1407    call wrtout(ab_out,message,'COLL')
1408    call appdig(pertcase,dtfil%fnamewff1,fiwf1i)
1409    call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o)
1410 
1411 !  Allocate 1st-order PAW occupancies (rhoij1)
1412    if (psps%usepaw==1) then
1413      ABI_DATATYPE_ALLOCATE(pawrhoij1,(my_natom))
1414      call pawrhoij_nullify(pawrhoij1)
1415      cplex_rhoij=max(cplex,dtset%pawcpxocc)
1416      nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
1417      use_rhoijim = 1
1418 !    use_rhoijim = 0 ; if (sum(dtset%qptn(1:3)**2)>1.d-14) use_rhoijim = 1
1419      call pawrhoij_alloc(pawrhoij1,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
1420 &     dtset%typat,use_rhoijim=use_rhoijim,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,&
1421 &     mpi_atmtab=mpi_enreg%my_atmtab)
1422      if (cplex_rhoij/=hdr%pawrhoij(1)%cplex) then
1423 !      Eventually reallocate hdr%pawrhoij
1424        call pawrhoij_free(hdr%pawrhoij)
1425        call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
1426 &       dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1427      end if
1428    else
1429      ABI_DATATYPE_ALLOCATE(pawrhoij1,(0))
1430    end if
1431 
1432 !  Initialize 1st-order wavefunctions
1433    formeig=1; ask_accurate=0; optorth=0
1434 ! NB: 4 Sept 2013: this was introducing a bug - for ieig2rf ==0 the dim was being set to 1 and passed to dfpt_vtowfk
1435    dim_eig2rf=0
1436    if ((dtset%ieig2rf > 0 .and. dtset%ieig2rf/=2) .or. dtset%efmas > 0) then
1437      dim_eig2rf=1
1438    end if
1439    mcg1=mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol
1440    if (one*mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol > huge(1)) then
1441      write (message,'(4a, 5(a,i0), 2a)')&
1442 &     "Default integer is not wide enough to store the size of the GS wavefunction array (WFK1, mcg1).",ch10,&
1443 &     "Action: increase the number of processors. Consider also OpenMP threads.",ch10,&
1444 &     "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mk1mem_rbz: ",&
1445 &     mk1mem_rbz, "nsppol: ",dtset%nsppol,ch10,&
1446 &     'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode'
1447      MSG_ERROR(message)
1448    end if
1449    ABI_STAT_ALLOCATE(cg1,(2,mcg1), ierr)
1450    ABI_CHECK(ierr==0, "out of memory in cg1")
1451    if (.not.kramers_deg) then
1452      mcg1mq=mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol
1453      ABI_STAT_ALLOCATE(cg1_mq,(2,mcg1mq), ierr)
1454      ABI_CHECK(ierr==0, "out of memory in cg1_mq")
1455    end if
1456 
1457    ABI_ALLOCATE(cg1_active,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1458    ABI_ALLOCATE(gh1c_set,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1459    ABI_ALLOCATE(gh0c1_set,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1460    if (.not.kramers_deg) then
1461      ABI_ALLOCATE(cg1_active_mq,(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1462      ABI_ALLOCATE(gh1c_set_mq,(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1463      ABI_ALLOCATE(gh0c1_set_mq,(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1464    end if
1465 !  XG090606 This is needed in the present 5.8.2 , for portability for the pathscale machine.
1466 !  However, it is due to a bug to be corrected by Paul Boulanger. When the bug will be corrected,
1467 !  this line should be removed.
1468    if(mk1mem_rbz/=0 .and. dtset%ieig2rf/=0)then
1469      cg1_active=zero
1470      gh1c_set=zero
1471      gh0c1_set=zero
1472      if (.not.kramers_deg) then
1473        cg1_active_mq=zero
1474        gh1c_set_mq=zero
1475        gh0c1_set_mq=zero
1476      end if
1477    end if
1478    ABI_ALLOCATE(eigen1,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol))
1479    ABI_ALLOCATE(resid,(dtset%mband*nkpt_rbz*dtset%nsppol))
1480    call timab(144,1,tsec)
1481    call status(pertcase,dtfil%filstat,iexit,level,'call inwffil  ')
1482    call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
1483 &   formeig,hdr,&
1484 &   dtfil%ireadwf,istwfk_rbz,kg1,kpq_rbz,dtset%localrdwf,&
1485 &   dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,&
1486 &   dtset%nsppol,nsym1,occ_rbz,optorth,&
1487 &   symaf1,symrl1,tnons1,dtfil%unkg1,wff1,wffnow,dtfil%unwff1,&
1488 &   fiwf1i,wvl)
1489    call timab(144,2,tsec)
1490 !  Close wff1 (filename fiwf1i), if it was ever opened (in inwffil)
1491    if (dtfil%ireadwf==1) then
1492      call WffClose(wff1,ierr)
1493    end if
1494    if(.not.kramers_deg) then
1495      ABI_ALLOCATE(eigen1_mq,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol))
1496      ABI_ALLOCATE(resid_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1497      !initialize cg1_mq:
1498      call timab(144,1,tsec)
1499      call status(pertcase,dtfil%filstat,iexit,level,'call inwffil  ')
1500      call inwffil(ask_accurate,cg1_mq,dtset,dtset%ecut,ecut_eff,eigen1_mq,dtset%exchn2n3d,&
1501 &     formeig,hdr,&
1502 &     dtfil%ireadwf,istwfk_rbz,kg1_mq,kmq_rbz,dtset%localrdwf,&
1503 &     dtset%mband,mcg1mq,mk1mem_rbz,mpi_enreg,mpw1_mq,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1_mq,&
1504 &     dtset%nsppol,nsym1,occ_rbz,optorth,&
1505 &     symaf1,symrl1,tnons1,dtfil%unkg1,wff1,wffnow,dtfil%unwff1,&
1506 &     fiwf1i,wvl)
1507      call timab(144,2,tsec)
1508 !    Close wff1 (filename fiwf1i), if it was ever opened (in inwffil)
1509      if (dtfil%ireadwf==1) then
1510        call WffClose(wff1,ierr)
1511      end if
1512    end if
1513 
1514 !  Eventually reytrieve 1st-order PAW occupancies from file header
1515    if (psps%usepaw==1.and.dtfil%ireadwf/=0) then
1516      call pawrhoij_copy(hdr%pawrhoij,pawrhoij1,comm_atom=mpi_enreg%comm_atom , &
1517 &     mpi_atmtab=mpi_enreg%my_atmtab)
1518    end if
1519 
1520 !  In case of electric field, or 2nd order perturbation : open the ddk (or dE) wf file(s)
1521    if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and. &
1522 &   (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and. &
1523 &   dtset%berryopt/= 7.and.dtset%berryopt/=14.and. &
1524 &   dtset%berryopt/=16.and.dtset%berryopt/=17)) .or. &
1525 &   ipert==dtset%natom+10.or.ipert==dtset%natom+11) then
1526      if (ipert<dtset%natom+10) then
1527        ! 1st order or berry phase, one direction
1528        nwffile=1 ! one direction needed
1529        file_index(1)=idir+dtset%natom*3
1530        fnamewff(1)=dtfil%fnamewffddk
1531      else if (ipert==dtset%natom+10) then
1532        ! 2nd order (k,k)
1533        if (idir<=3) then ! one direction needed
1534          nwffile=1
1535          file_index(1)=idir+dtset%natom*3
1536          fnamewff(1)=dtfil%fnamewffddk
1537        else ! two directions needed
1538          nwffile=2
1539          file_index(1)=idir1+dtset%natom*3
1540          file_index(2)=idir2+dtset%natom*3
1541          fnamewff(1)=dtfil%fnamewffddk
1542          fnamewff(2)=dtfil%fnamewffddk
1543        end if
1544      else if (ipert==dtset%natom+11) then
1545        ! 2nd order (k,E)
1546        nwffile=3 ! dk, dE and dkdk
1547        idir_dkdk = idir
1548        if(idir_dkdk>6) idir_dkdk = idir_dkdk - 3
1549        file_index(1)=idir_dkdk +(dtset%natom+6)*3 ! dkdk
1550        file_index(2)=idir2+(dtset%natom+1)*3 ! defld file (dir2)
1551        file_index(3)=idir1+dtset%natom*3     ! ddk file (dir1)
1552        fnamewff(1)=dtfil%fnamewffdkdk
1553        fnamewff(2)=dtfil%fnamewffdelfd
1554        fnamewff(3)=dtfil%fnamewffddk
1555        if (idir>3) then
1556          nwffile=4
1557          file_index(4)=idir2+dtset%natom*3   ! ddk file (dir2)
1558          fnamewff(4)=dtfil%fnamewffddk
1559        end if
1560      end if
1561      do ii=1,nwffile
1562        call appdig(file_index(ii),fnamewff(ii),fiwfddk)
1563        ! Checking the existence of data file
1564        if (.not. file_exists(fiwfddk)) then
1565          ! Trick needed to run Abinit test suite in netcdf mode.
1566          if (file_exists(nctk_ncify(fiwfddk))) then
1567            write(message,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name."
1568            call wrtout(std_out,message,'COLL')
1569            fiwfddk = nctk_ncify(fiwfddk)
1570          end if
1571          if (.not. file_exists(fiwfddk)) then
1572            MSG_ERROR('Missing file: '//TRIM(fiwfddk))
1573          end if
1574        end if
1575        write(message,'(2a)')'-dfpt_looppert : read the wavefunctions from file: ',trim(fiwfddk)
1576        call wrtout(std_out,message,'COLL')
1577        call wrtout(ab_out,message,'COLL')
1578 !      Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50)
1579        call wfk_open_read(ddk_f(ii),fiwfddk,formeig1,dtset%iomode,dtfil%unddk+(ii-1),spaceComm)
1580      end do
1581    end if
1582 
1583 !  Get first-order local potentials and 1st-order core correction density change
1584 !  (do NOT include xccc3d1 in vpsp1 : this will be done in dfpt_scfcv because vpsp1
1585 !  might become spin-polarized)
1586 
1587    n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
1588    ABI_ALLOCATE(xccc3d1,(cplex*n3xccc))
1589    ABI_ALLOCATE(vpsp1,(cplex*nfftf))
1590 
1591 !  PAW: compute Vloc(1) and core(1) together in reciprocal space
1592 !  --------------------------------------------------------------
1593    if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
1594      ndir=1
1595      call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,istr,ipert,&
1596 &     mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
1597 &     ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
1598 &     atmrhor1=xccc3d1,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl)
1599    else
1600 
1601 !    Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
1602 !    ------------------------------------------------------------------------------
1603 
1604      if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
1605 !      Section for strain perturbation
1606        call vlocalstr(gmet,gprimd,gsqcut,istr,mgfftf,mpi_enreg,&
1607 &       psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,dtset%paral_kgb,ph1df,psps%qgrid_vl,&
1608 &       ucvol,psps%vlspl,vpsp1)
1609      else
1610        call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir,ipert,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
1611 &       nattyp,nfftf,ngfftf,ntypat,ngfftf(1),ngfftf(2),ngfftf(3),dtset%paral_kgb,ph1df,psps%qgrid_vl,&
1612 &       dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
1613        !SPr: need vpsp1 for -q as well, but for magnetic field it's zero, to be done later
1614      end if
1615 
1616      if(psps%n1xccc/=0)then
1617        call dfpt_mkcore(cplex,idir,ipert,dtset%natom,ntypat,ngfftf(1),psps%n1xccc,&
1618 &       ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
1619        !SPr: same here, need xccc3d1 for -q as well for phonon pert.. to be done later
1620      end if ! psps%n1xccc/=0
1621    end if ! usepaw
1622 
1623    eigen1(:)=zero; resid(:)=zero
1624    if(.not.kramers_deg) then
1625      eigen1_mq(:)=zero
1626      resid_mq(:)=zero
1627    end if
1628 !  Get starting charge density and Hartree + xc potential
1629    ABI_ALLOCATE(rhor1,(cplex*nfftf,nspden))
1630    ABI_ALLOCATE(rhog1,(2,nfftf))
1631 
1632    if(.not.kramers_deg) then
1633    !Case when first order spinors at both +q and -q are not related by symmetry (time and/or space inversion)
1634      ABI_ALLOCATE(rhor1_pq,(cplex*nfftf,nspden))
1635      ABI_ALLOCATE(rhog1_pq,(2,nfftf))
1636 
1637      ABI_ALLOCATE(rhor1_mq,(cplex*nfftf,nspden))
1638      ABI_ALLOCATE(rhog1_mq,(2,nfftf))
1639    end if
1640 
1641 !  can we get this set of gkk matrices from previously calculated rhog1 through a non-scf calculation?
1642    found_eq_gkk=.false.
1643    if (dtset%prepgkk == 1 .and. ipert <= dtset%natom) then
1644 !    NOTE: this does not take into account combinations e.g. x+y -> z
1645 !    if rhor1 add linearly this could be done...
1646      do icase_eq = 1, icase-1
1647        idir_eq = mod(icase_eq,3)
1648        if (idir_eq==0) idir_eq=3
1649        ipert_eq = ( (icase_eq-idir_eq) / 3 + 1)
1650 
1651 ! find sym which links old perturbation to present one
1652        do isym=1, nsym
1653          ! check that isym preserves qpt to begin with
1654          if (symq(4,1,isym) /= 1 .or. &
1655 &         symq(1,1,isym) /= 0 .or. &
1656 &         symq(2,1,isym) /= 0 .or. &
1657 &         symq(3,1,isym) /= 0      ) cycle
1658 
1659          eq_symop = dtset%symrel(:,:,isym)
1660          if (indsym(4,isym,ipert) == ipert_eq .and. &
1661 &         abs(eq_symop(idir,idir_eq)) == 1 .and. &
1662 &         sum(abs(eq_symop(:,idir_eq))) == 1) then
1663            found_eq_gkk = .true.
1664            exit
1665          end if
1666        end do ! isym
1667        if (found_eq_gkk) exit
1668      end do ! icase_eq
1669    end if ! check for prepgkk with symmetric pert
1670 
1671    if (found_eq_gkk) then
1672      write (message, '(a,l6,i6,2a,3i6,2a,3i6,a)')  &
1673 &     ' found_eq_gkk,isym = ', found_eq_gkk, isym, ch10, &
1674 &     ' idir,  ipert,  icase   =  ', idir, ipert, icase, ch10, &
1675 &     ' idireq iperteq icaseeq =  ', idir_eq, ipert_eq, icase_eq, ch10
1676      call wrtout(std_out,message,'COLL')
1677 !
1678 !    Make density for present perturbation, which is symmetric of 1 or more previous perturbations:
1679 !    rotate 1DEN arrays with symrel(isym) to produce rhog1_eq rhor1_eq
1680 !
1681      if (dtset%use_nonscf_gkk == 1) then
1682        call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, &
1683 &       rhor1_save(:,:,icase_eq), rhog1, rhor1, eq_symop, dtset%tnons(:,isym))
1684 
1685        rhor1 = rhor1 * eq_symop(idir,idir_eq)
1686 
1687 ! TODO: rotate rhoij in PAW case
1688 
1689        blkflg_save = blkflg
1690        dtset_tmp%iscf = -2
1691        iscf_mod = -2
1692        dtset_tmp%nstep = 1
1693        dtset_tmp%nline = 1
1694        if (abs(dtset_tmp%tolwfr) < 1.e-24) dtset_tmp%tolwfr = 1.e-24
1695        dtset_tmp%toldfe = zero
1696        dtset_tmp%toldff = zero
1697        dtset_tmp%tolrff = zero
1698        dtset_tmp%tolvrs = zero
1699        write (message, '(a,i6,a)') ' NOTE: doing GKK calculation for icase ', icase, ' with non-SCF calculation'
1700        call wrtout(std_out,message,'COLL')
1701        !call wrtout(ab_out,message,'COLL') ! decomment and update output files
1702 
1703      else ! do not use non-scf shortcut, but save rotated 1DEN for comparison
1704 ! saves the rotated rho, for later comparison with the full SCF rhor1: comment lines below for iscf = -2
1705        call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, &
1706 &       rhor1_save(:,:,icase_eq), rhog1, rhor1_save(:,:,icase), eq_symop, &
1707 &       dtset%tnons(:,isym))
1708        rhor1_save(:,:,icase) = rhor1_save(:,:,icase) * eq_symop(idir,idir_eq)
1709 
1710      end if ! force non scf calculation of other gkk, or not
1711    end if ! found an equiv perturbation for the gkk
1712 
1713    if ( (dtfil%ireadwf==0 .and. iscf_mod/=-4 .and. dtset%get1den==0 .and. dtset%ird1den==0) &
1714    .or. (iscf_mod== -3 .and. ipert/=dtset%natom+11) ) then
1715 !    NOTE : For ipert==natom+11, we want to read the 1st order density from a previous calculation
1716      rhor1(:,:)=zero ; rhog1(:,:)=zero
1717 !    PAW: rhoij have been set to zero in call to pawrhoij_alloc above
1718 
1719      init_rhor1 = ((ipert>=1 .and. ipert<=dtset%natom).or.ipert==dtset%natom+5)
1720      ! This section is needed in order to maintain the old behavior and pass the automatic tests
1721      if (psps%usepaw == 0) then
1722        init_rhor1 = init_rhor1 .and. all(psps%nctab(:ntypat)%has_tvale)
1723      else
1724        init_rhor1 = .False.
1725      end if
1726 
1727      if (init_rhor1) then
1728 
1729        if(ipert/=dtset%natom+5) then
1730 
1731          ! Initialize rhor1 and rhog1 from the derivative of atomic densities/gaussians.
1732          ndir = 1; optn2 = 3
1733          if (psps%usepaw==1) then
1734            ! FIXME: Here there's a bug because has_tvale == 0 or 1 instead of 2
1735            if (all(pawtab(:ntypat)%has_tvale/=0)) optn2=2
1736          else if (psps%usepaw==0) then
1737            if (all(psps%nctab(:ntypat)%has_tvale)) optn2=2
1738          end if
1739 
1740          if (optn2 == 3) then
1741            call wrtout(std_out," Initializing rhor1 from atom-centered gaussians", "COLL")
1742            ABI_ALLOCATE(gauss,(2,ntypat))
1743            call atom_gauss(ntypat, dtset%densty, psps%ziontypat, psps%znucltypat, gauss)
1744 
1745            call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
1746            mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
1747            ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
1748            atmrhor1=rhor1,optn_in=1,optn2_in=3,gauss=gauss)
1749 
1750            ABI_FREE(gauss)
1751          else
1752            call wrtout(std_out," Initializing rhor1 from valence densities taken from pseudopotential files", "COLL")
1753            call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
1754            mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
1755            ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
1756            atmrhor1=rhor1,optn_in=1,optn2_in=2)
1757          end if
1758 
1759        else
1760          !Magnetic field perturbation
1761          call wrtout(std_out," Initializing rhor1 guess based on the ground state XC magnetic field", "COLL")
1762 
1763          call dfpt_init_mag1(ipert,idir,rhor1,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc)
1764 
1765          if(.not.kramers_deg) then
1766            rhor1_pq=rhor1
1767            call dfpt_init_mag1(ipert,idir,rhor1_mq,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc)
1768          end if
1769        end if
1770 
1771        call fourdp(cplex,rhog1,rhor1,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1772        if (.not.kramers_deg) then
1773          !call fourdp(cplex,rhog1_pq,rhor1_pq,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1774          rhog1_pq=rhog1
1775          call fourdp(cplex,rhog1_mq,rhor1_mq,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1776        end if
1777      end if
1778 
1779    else  ! rhor1 not being forced to 0.0
1780      if(iscf_mod>0) then
1781 !      cplex=2 gets the complex density, =1 only real part
1782        if (psps%usepaw==1) then
1783 !        Be careful: in PAW, rho does not include the 1st-order compensation
1784 !        density (to be added in dfpt_scfcv.F90) !
1785          ABI_ALLOCATE(rho1wfg,(2,dtset%nfft))
1786          ABI_ALLOCATE(rho1wfr,(dtset%nfft,nspden))
1787          call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,&
1788 &         kg,kg1,dtset%mband,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,&
1789 &         dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
1790 &         occ_rbz,dtset%paral_kgb,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,ucvol,&
1791 &         wtk_rbz)
1792          call transgrid(cplex,mpi_enreg,nspden,+1,1,1,dtset%paral_kgb,pawfgr,rho1wfg,rhog1,rho1wfr,rhor1)
1793          ABI_DEALLOCATE(rho1wfg)
1794          ABI_DEALLOCATE(rho1wfr)
1795        else
1796          !SPr: need to modify dfpt_mkrho to taken into account q,-q and set proper formulas when +q and -q spinors are
1797          !     related
1798          call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,&
1799 &         kg,kg1,dtset%mband,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,&
1800 &         dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
1801 &         occ_rbz,dtset%paral_kgb,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,ucvol,&
1802 &         wtk_rbz)
1803        end if
1804 
1805      else if (.not. found_eq_gkk) then
1806        ! negative iscf_mod and no symmetric rotation of rhor1
1807        ! Read rho1(r) from a disk file and broadcast data.
1808        rdwr=1;rdwrpaw=psps%usepaw;if(dtfil%ireadwf/=0) rdwrpaw=0
1809        if (ipert/=dtset%natom+11) then
1810          call appdig(pertcase,dtfil%fildens1in,fiden1i)
1811        else ! For ipert==natom+11, we want to read the 1st order density from a previous calculation
1812          call appdig(idir2+(dtset%natom+1)*3,dtfil%fildens1in,fiden1i)
1813        end if
1814 !       call appdig(pertcase,dtfil%fildens1in,fiden1i)
1815        call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor1, &
1816        hdr_den, pawrhoij1, spaceComm, check_hdr=hdr)
1817        etotal = hdr_den%etot; call hdr_free(hdr_den)
1818 
1819 !      Compute up+down rho1(G) by fft
1820        ABI_ALLOCATE(work,(cplex*nfftf))
1821        work(:)=rhor1(:,1)
1822        call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1823        ABI_DEALLOCATE(work)
1824      end if ! rhor1 generated or read in from file
1825 
1826    end if ! rhor1 set to 0 or read in from file
1827 
1828 !  Check whether exiting was required by the user.
1829 !  If found then do not start minimization steps
1830    openexit=1 ; if(dtset%chkexit==0) openexit=0
1831    call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
1832 !  If immediate exit, and wavefunctions were not read, must zero eigenvalues
1833    if (iexit/=0) eigen1(:)=zero
1834    if (iexit/=0.and.(.not.kramers_deg)) eigen1_mq(:)=zero
1835 
1836    if (iexit==0) then
1837      _IBM6("before dfpt_scfcv")
1838 
1839 !    Main calculation: get 1st-order wavefunctions from Sternheimer equation (SCF cycle)
1840 !    if ipert==natom+10 or natom+11 : get 2nd-order wavefunctions
1841      if (kramers_deg) then
1842        call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,&
1843 &       dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,&
1844 &       d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
1845 &       ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,&
1846 &       enl0,enl1,eovl1,epaw1,etotal,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,&
1847 &       indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
1848 &       kg,kg1,kpt_rbz,kxc,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,&
1849 &       mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,&
1850 &       nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,&
1851 &       npwarr,npwar1,nspden,&
1852 &       nsym1,n3xccc,occkq,occ_rbz,&
1853 &       paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,&
1854 &       pertcase,phnons1,ph1d,ph1df,prtbbb,psps,&
1855 &       dtset%qptn,resid,residm,rhog,rhog1,&
1856 &       rhor,rhor1,rprimd,symaf1,symrc1,symrl1,&
1857 &       usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,&
1858 &       wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,&
1859 &       kramers_deg)
1860      else
1861        call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,&
1862 &       dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,&
1863 &       d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
1864 &       ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,&
1865 &       enl0,enl1,eovl1,epaw1,etotal,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,&
1866 &       indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
1867 &       kg,kg1,kpt_rbz,kxc,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,&
1868 &       mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,&
1869 &       nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,&
1870 &       npwarr,npwar1,nspden,&
1871 &       nsym1,n3xccc,occkq,occ_rbz,&
1872 &       paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,&
1873 &       pertcase,phnons1,ph1d,ph1df,prtbbb,psps,&
1874 &       dtset%qptn,resid,residm,rhog,rhog1,&
1875 &       rhor,rhor1,rprimd,symaf1,symrc1,symrl1,&
1876 &       usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,&
1877 &       wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,&
1878 &       kramers_deg,&
1879 &       cg_mq=cg_mq,cg1_mq=cg1_mq,cg1_active_mq=cg1_active_mq,docckde_mq=docckde_mq,eigen_mq=eigen_mq,&
1880 &       eigen1_mq=eigen1_mq,gh0c1_set_mq=gh0c1_set_mq,gh1c_set_mq=gh1c_set_mq,&
1881 &       kg1_mq=kg1_mq,npwar1_mq=npwar1_mq,occk_mq=occk_mq,resid_mq=resid_mq,residm_mq=residm_mq,&
1882 &       rhog1_pq=rhog1_pq,rhog1_mq=rhog1_mq,rhor1_pq=rhor1_pq,rhor1_mq=rhor1_mq)
1883      end if
1884 
1885      _IBM6("after dfpt_scfcv")
1886 
1887 !    2nd-order eigenvalues stuff
1888      if (dtset%ieig2rf>0) then
1889        if (first_entry) then
1890          nullify(eigen1_pert)
1891          first_entry = .false.
1892        end if
1893        if (.not.associated(eigen1_pert)) then
1894          ABI_ALLOCATE(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert))
1895          ABI_STAT_ALLOCATE(cg1_pert,(2,mpw1*nspinor*dtset%mband*mk1mem_rbz*nsppol*dim_eig2rf,3,mpert),ierr)
1896          ABI_CHECK(ierr==0, "out-of-memory in cg1_pert")
1897          ABI_ALLOCATE(gh0c1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
1898          ABI_ALLOCATE(gh1c_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
1899          ABI_ALLOCATE(kpt_rbz_pert,(3,nkpt_rbz))
1900          ABI_ALLOCATE(npwarr_pert,(nkpt_rbz,mpert))
1901          ABI_ALLOCATE(npwar1_pert,(nkpt_rbz,mpert))
1902          ABI_ALLOCATE(npwtot_pert,(nkpt_rbz,mpert))
1903          eigen1_pert(:,:,:) = zero
1904          cg1_pert(:,:,:,:) = zero
1905          gh0c1_pert(:,:,:,:) = zero
1906          gh1c_pert(:,:,:,:) = zero
1907          npwar1_pert (:,:) = 0
1908          npwarr_pert (:,:) = 0
1909          kpt_rbz_pert = kpt_rbz
1910        end if
1911        clflg(idir,ipert)=1
1912        eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:)
1913        if(dtset%ieig2rf==1.or.dtset%ieig2rf==3.or.dtset%ieig2rf==4.or.dtset%ieig2rf==5) then
1914          cg1_pert(:,:,idir,ipert)=cg1_active(:,:)
1915          gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:)
1916          gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:)
1917        end if
1918        npwarr_pert(:,ipert)=npwarr(:)
1919        npwar1_pert(:,ipert)=npwar1(:)
1920        npwtot_pert(:,ipert)=npwtot(:)
1921      end if
1922 !    2nd-order eigenvalues stuff for EFMAS
1923      if (dtset%efmas>0) then
1924        if (first_entry) then
1925          nullify(eigen1_pert)
1926          first_entry = .false.
1927        end if
1928        if (.not.associated(eigen1_pert)) then
1929          ABI_ALLOCATE(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert))
1930          ABI_ALLOCATE(cg1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
1931          ABI_ALLOCATE(gh0c1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
1932          ABI_ALLOCATE(gh1c_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
1933          ABI_ALLOCATE(kpt_rbz_pert,(3,nkpt_rbz))
1934          ABI_ALLOCATE(npwarr_pert,(nkpt_rbz,mpert))
1935          eigen1_pert(:,:,:) = zero
1936          cg1_pert(:,:,:,:) = zero
1937          gh0c1_pert(:,:,:,:) = zero
1938          gh1c_pert(:,:,:,:) = zero
1939          npwarr_pert (:,:) = 0
1940          kpt_rbz_pert = kpt_rbz
1941          ABI_ALLOCATE(cg0_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1942          cg0_pert = cg
1943        end if
1944        eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:)
1945        cg1_pert(:,:,idir,ipert)=cg1_active(:,:)
1946        gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:)
1947        gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:)
1948        npwarr_pert(:,ipert)=npwarr(:)
1949      end if
1950      ABI_DEALLOCATE(gh1c_set)
1951      ABI_DEALLOCATE(gh0c1_set)
1952      ABI_DEALLOCATE(cg1_active)
1953      if(.not.kramers_deg) then
1954        ABI_DEALLOCATE(gh1c_set_mq)
1955        ABI_DEALLOCATE(gh0c1_set_mq)
1956        ABI_DEALLOCATE(cg1_active_mq)
1957      end if
1958 
1959    end if ! End of the check of hasty exit
1960 
1961    call timab(146,1,tsec)
1962 
1963 !  Print out message at the end of the iterations
1964    write(message, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,&
1965 &   ' ----iterations are completed or convergence reached----',ch10
1966    call wrtout(ab_out,message,'COLL')
1967    call wrtout(std_out,  message,'COLL')
1968 
1969 !  Print _gkk file for this perturbation
1970    if (dtset%prtgkk == 1) then
1971      call appdig(3*(ipert-1)+idir,dtfil%fnameabo_gkk,gkkfilnam)
1972      nmatel = dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol
1973      ABI_ALLOCATE(phasecg, (2, nmatel))
1974      call getcgqphase(dtset, timrev, cg,  mcg,  cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg)
1975      phasecg(1,:) = one
1976      phasecg(2,:) = zero
1977 ! NB: phasecg not actually used in outgkk for the moment (2013/08/15)
1978      call outgkk(bantot_rbz, nmatel,gkkfilnam,eigen0,eigen1,hdr0,hdr,mpi_enreg,phasecg)
1979      ABI_DEALLOCATE(phasecg)
1980 
1981 #ifdef HAVE_NETCDF
1982      ! Reshape eigen1 into gkk for netCDF output
1983      ABI_STAT_ALLOCATE(gkk,(2*dtset%mband*dtset%nsppol,dtset%nkpt,1,1,dtset%mband), ierr)
1984      ABI_CHECK(ierr==0, "out-of-memory in gkk")
1985      gkk(:,:,:,:,:) = zero
1986      mband = dtset%mband
1987      band_index = 0
1988      band2tot_index = 0
1989      do isppol=1,dtset%nsppol
1990        do ikpt =1,nkpt_rbz
1991          do iband=1,dtset%mband
1992            do jband=1,dtset%mband
1993              eig1_r = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index)
1994              eig1_i = eigen1(2*jband+(iband-1)*2*mband+band2tot_index)
1995              gkk(2*iband-1+2*band_index,ikpt,1,1,jband) = &
1996 &             gkk(2*iband-1+2*band_index,ikpt,1,1,jband) + eig1_r
1997              gkk(2*iband+2*band_index,ikpt,1,1,jband) = &
1998 &             gkk(2*iband+2*band_index,ikpt,1,1,jband) + eig1_i
1999            end do !jband
2000          end do !iband
2001          band2tot_index = band2tot_index + 2*mband**2
2002        end do !ikpt
2003        band_index = band_index + mband
2004      end do !isppol
2005 
2006      ! Initialize ggk_ebands to write in the GKK.nc file
2007      ! MG FIXME: Here there's a bug because eigen0 is dimensioned with nkpt_rbz i.e. IBZ(q)
2008      ! but the ebands_t object is constructed with dimensions taken from hdr0 i.e. the IBZ(q=0).
2009      bantot= dtset%mband*dtset%nkpt*dtset%nsppol
2010      call ebands_init(bantot,gkk_ebands,dtset%nelect,doccde,eigen0,hdr0%istwfk,hdr0%kptns,&
2011 &     hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,&
2012 &     hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,&
2013 &     hdr0%charge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, &
2014 &     hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk)
2015 
2016      ! Init a gkk_t object
2017      call gkk_init(gkk,gkk2d,dtset%mband,dtset%nsppol,nkpt_rbz,1,1)
2018 
2019      ! Write the netCDF file.
2020      fname = strcat(gkkfilnam,".nc")
2021      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file")
2022      NCF_CHECK(crystal_ncwrite(crystal, ncid))
2023      NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid))
2024      call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid)
2025      NCF_CHECK(nf90_close(ncid))
2026 
2027      ! Free memory
2028      ABI_DEALLOCATE(gkk)
2029      call gkk_free(gkk2d)
2030      call ebands_free(gkk_ebands)
2031 #endif
2032    end if
2033 
2034    if (dtset%prepgkk == 1 .and. found_eq_gkk) then
2035      if (dtset%use_nonscf_gkk == 1) then
2036 !      Restore old values of SCF cycle parameters
2037        iscf_mod = iscf_mod_save
2038        dtset_tmp%iscf = iscf_mod_save
2039        dtset_tmp%nstep = nstep_save
2040        dtset_tmp%nline = nline_save
2041        dtset_tmp%tolwfr = tolwfr_save
2042        dtset_tmp%toldfe = toldfe_save
2043        dtset_tmp%toldff = toldff_save
2044        dtset_tmp%tolrff = tolrff_save
2045        dtset_tmp%tolvrs = tolvrs_save
2046        blkflg = blkflg_save ! this ensures we do not use the (unconverged) 2DTE from this non scf run
2047 !      Save density for present perturbation, for future use in symmetric perturbations
2048        rhor1_save(:,:,icase) = rhor1
2049      else
2050        write (message, '(a,3E20.10)') 'norm diff = ', sum(abs(rhor1_save(:,:,icase) - rhor1)), &
2051 &       sum(abs(rhor1)), sum(abs(rhor1_save(:,:,icase) - rhor1))/sum(abs(rhor1))
2052        call wrtout(std_out,  message,'COLL')
2053      end if
2054    end if
2055 
2056    ! Write wavefunctions file only if convergence was not achieved.
2057    write_1wfk = .True.
2058    if (dtset%prtwf==-1 .and. dfpt_scfcv_retcode == 0) then
2059      write_1wfk = .False.
2060      call wrtout(ab_out," dfpt_looppert: DFPT cycle converged with prtwf=-1. Will skip output of the 1st-order WFK file.","COLL")
2061    end if
2062 
2063    if (write_1wfk) then
2064      ! Output 1st-order wavefunctions in file
2065      call status(pertcase,dtfil%filstat,iexit,level,'call outwf    ')
2066      call outwf(cg1,dtset,psps,eigen1,fiwf1o,hdr,kg1,kpt_rbz,&
2067 &     dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,dtset%natom,nband_rbz,&
2068 &     nkpt_rbz,npwar1,dtset%nsppol,&
2069 &     occ_rbz,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr)
2070    end if
2071 
2072 #ifdef HAVE_NETCDF
2073   ! Output DDK file in netcdf format.
2074    if (me == master .and. ipert == dtset%natom + 1) then
2075      fname = strcat(dtfil%filnam_ds(4), "_EVK.nc")
2076      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
2077     ! Have to build hdr on k-grid with info about perturbation.
2078      call hdr_copy(hdr0, hdr_tmp)
2079      hdr_tmp%kptopt = dtset%kptopt
2080      hdr_tmp%pertcase = pertcase
2081      hdr_tmp%qptn = dtset%qptn(1:3)
2082      NCF_CHECK(hdr_ncwrite(hdr_tmp, ncid, 43, nc_define=.True.))
2083      call hdr_free(hdr_tmp)
2084      NCF_CHECK(crystal_ncwrite(crystal, ncid))
2085      NCF_CHECK(ebands_ncwrite(ebands_k, ncid))
2086      ncerr = nctk_def_arrays(ncid, [ &
2087      nctkarr_t('h1_matrix_elements', "dp", &
2088      "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
2089      NCF_CHECK(ncerr)
2090      NCF_CHECK(nctk_set_datamode(ncid))
2091      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), eigen1, &
2092      count=[2, dtset%mband, dtset%mband, nkpt_rbz, dtset%nsppol])
2093      NCF_CHECK(ncerr)
2094      NCF_CHECK(nf90_close(ncid))
2095    end if
2096 #endif
2097 
2098 
2099 !  If the perturbation is d/dk, evaluate the f-sum rule.
2100    if (ipert==dtset%natom+1 )then
2101 !    Note : the factor of two is related to the difference
2102 !    between Taylor expansion and perturbation expansion
2103 !    Note : this expression should be modified for ecutsm.
2104 !    Indeed, the present one will NOT tend to 1.0_dp.
2105      ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect
2106      fsum=-ek1/ek2
2107      if(dtset%ecutsm<tol6)then
2108        write(message, '(a,es20.10,a,a,es20.10)' ) &
2109 &       ' dfpt_looppert : ek2=',ek2,ch10,&
2110 &       '          f-sum rule ratio=',fsum
2111      else
2112        write(message, '(a,es20.10,a,a,es20.10,a)' ) &
2113 &       ' dfpt_looppert : ek2=',ek2,ch10,&
2114 &       '          f-sum rule ratio=',fsum,' (note : ecutsm/=0)'
2115      end if
2116      call wrtout(std_out,message,'COLL')
2117      call wrtout(ab_out,message,'COLL')
2118 !    Write the diagonal elements of the dH/dk operator, after averaging over degenerate states
2119      ABI_ALLOCATE(eigen1_mean,(dtset%mband*nkpt_rbz*dtset%nsppol))
2120      call eigen_meandege(eigen0,eigen1,eigen1_mean,dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,1)
2121      option=4
2122      if (me == master) then
2123        call prteigrs(eigen1_mean,dtset%enunit,fermie,dtfil%fnametmp_1wf1_eig,ab_out,iscf_mod,kpt_rbz,dtset%kptopt,&
2124 &       dtset%mband,nband_rbz,nkpt_rbz,dtset%nnsclo,dtset%nsppol,occ_rbz,dtset%occopt,&
2125 &       option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk_rbz)
2126      end if
2127      ABI_DEALLOCATE(eigen1_mean)
2128    end if
2129 
2130 !  Print the energies
2131    if (dtset%nline/=0 .or. dtset%nstep/=0)then
2132      call dfpt_prtene(dtset%berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
2133 &     ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,eovl1,epaw1,evdw,exc1,ab_out,&
2134 &     ipert,dtset%natom,psps%usepaw,usevdw)
2135    end if
2136 
2137    if(mpi_enreg%paral_pert==1) then
2138      if (ipert_me < npert_me -1) then
2139        call hdr_free(hdr0)
2140      else
2141        eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
2142      end if
2143      ipert_me = ipert_me +1
2144    else
2145      if (icase == ipert_cnt) then
2146        eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
2147      else
2148        call hdr_free(hdr0)
2149      end if
2150    end if
2151 
2152 !Release the temporary arrays (for k, k+q and 1st-order)
2153    ABI_DEALLOCATE(cg)
2154    ABI_DEALLOCATE(cgq)
2155    ABI_DEALLOCATE(cg1)
2156    ABI_DEALLOCATE(docckqde)
2157    if(.not.kramers_deg) then
2158      ABI_DEALLOCATE(cg_mq)
2159      ABI_DEALLOCATE(cg1_mq)
2160      ABI_DEALLOCATE(docckde_mq)
2161    end if
2162    ABI_DEALLOCATE(doccde_rbz)
2163    ABI_DEALLOCATE(eigen0)
2164    ABI_DEALLOCATE(eigenq)
2165    ABI_DEALLOCATE(eigen1)
2166    ABI_DEALLOCATE(kpq)
2167    if(.not.kramers_deg) then
2168      ABI_DEALLOCATE(eigen_mq)
2169      ABI_DEALLOCATE(eigen1_mq)
2170      ABI_DEALLOCATE(kmq)
2171    end if
2172    ABI_DEALLOCATE(indkpt1)
2173    ABI_DEALLOCATE(indsy1)
2174    ABI_DEALLOCATE(istwfk_rbz)
2175    ABI_DEALLOCATE(irrzon1)
2176    ABI_DEALLOCATE(kg)
2177    ABI_DEALLOCATE(kg1)
2178    ABI_DEALLOCATE(kpq_rbz)
2179    if(.not.kramers_deg) then
2180      ABI_DEALLOCATE(kg1_mq)
2181      ABI_DEALLOCATE(kmq_rbz)
2182    end if
2183    ABI_DEALLOCATE(kpt_rbz)
2184    ABI_DEALLOCATE(nband_rbz)
2185    ABI_DEALLOCATE(npwarr)
2186    ABI_DEALLOCATE(npwar1)
2187    ABI_DEALLOCATE(npwtot)
2188    ABI_DEALLOCATE(npwtot1)
2189    ABI_DEALLOCATE(occkq)
2190    ABI_DEALLOCATE(occ_rbz)
2191    ABI_DEALLOCATE(phnons1)
2192    ABI_DEALLOCATE(resid)
2193    ABI_DEALLOCATE(rhog1)
2194    ABI_DEALLOCATE(rhor1)
2195    ABI_DEALLOCATE(symaf1)
2196    ABI_DEALLOCATE(symrc1)
2197    ABI_DEALLOCATE(symrl1)
2198    ABI_DEALLOCATE(tnons1)
2199    ABI_DEALLOCATE(wtk_rbz)
2200    ABI_DEALLOCATE(xccc3d1)
2201    ABI_DEALLOCATE(vpsp1)
2202    ABI_DEALLOCATE(ylm)
2203    ABI_DEALLOCATE(ylm1)
2204    ABI_DEALLOCATE(ylmgr)
2205    ABI_DEALLOCATE(ylmgr1)
2206    if(.not.kramers_deg) then
2207      ABI_DEALLOCATE(npwar1_mq)
2208      ABI_DEALLOCATE(npwtot1_mq)
2209      ABI_DEALLOCATE(occk_mq)
2210      ABI_DEALLOCATE(resid_mq)
2211      ABI_DEALLOCATE(rhor1_pq)
2212      ABI_DEALLOCATE(rhor1_mq)
2213      ABI_DEALLOCATE(rhog1_pq)
2214      ABI_DEALLOCATE(rhog1_mq)
2215    end if
2216    if (psps%usepaw==1) then
2217      call pawang_free(pawang1)
2218      call pawrhoij_free(pawrhoij1)
2219      if (usecprj==1) then
2220        call pawcprj_free(cprj)
2221        call pawcprj_free(cprjq)
2222      end if
2223    end if
2224    ABI_DATATYPE_DEALLOCATE(pawrhoij1)
2225    ABI_DATATYPE_DEALLOCATE(cprjq)
2226    ABI_DATATYPE_DEALLOCATE(cprj)
2227    if(xmpi_paral==1)  then
2228      ABI_DEALLOCATE(mpi_enreg%proc_distrb)
2229      ABI_DEALLOCATE(mpi_enreg%my_kpttab)
2230    end if
2231    call hdr_free(hdr)
2232 
2233    ! Clean band structure datatypes (should use it more in the future !)
2234    call ebands_free(ebands_k)
2235    call ebands_free(ebands_kq)
2236    if(.not.kramers_deg) then
2237      call ebands_free(ebands_kmq)
2238    end if
2239 
2240 !  %%%% Parallelization over perturbations %%%%%
2241 !  *Redefine output/log files
2242    call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,npert_io,&
2243 &   mpi_enreg%paral_pert,0)
2244 
2245    call timab(146,2,tsec)
2246    if(iexit/=0) exit
2247  end do ! End loop on perturbations
2248  ABI_DEALLOCATE(zeff)
2249 
2250 !%%%% Parallelization over perturbations %%%%%
2251 !*Restore default communicators
2252  call unset_pert_comm(mpi_enreg)
2253 !*Gather output/log files
2254  ABI_ALLOCATE(dyn,(npert_io))
2255  if (npert_io>0) dyn=1
2256  call localrdfile(mpi_enreg%comm_pert,mpi_enreg%comm_world,.true.,npert_io,&
2257 & mpi_enreg%paral_pert,0,dyn)
2258  ABI_DEALLOCATE(dyn)
2259 !*Restore PAW on-site data
2260  if (paral_pert_inplace) then
2261    call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
2262 &   paw_an,paw_ij,pawfgrtab,pawrhoij)
2263  else
2264    call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
2265 &   paw_an,paw_ij,pawfgrtab,pawrhoij,&
2266 &   paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,&
2267 &   pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert)
2268  end if
2269 
2270 !#################################################################################
2271 !Calculate the second-order eigenvalues for a wavevector Q
2272 
2273  call timab(147,1,tsec)
2274  smdelta = dtset%smdelta
2275  bdeigrf = dtset%bdeigrf
2276  if(dtset%bdeigrf == -1) bdeigrf = dtset%mband
2277 
2278  if(dtset%ieig2rf > 0) then
2279 
2280    if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
2281 &   (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
2282      call wrtout(std_out,'Reading the dense grid WF file',"COLL")
2283 !    We get the Abinit header of the file hdr_fine as ouput
2284 !    We get eigenq_fine(mband,hdr_fine%nkpt,hdr_fine%nsppol) as ouput
2285      fname = dtfil%fnameabi_wfkfine
2286      if (dtset%iomode == IO_MODE_ETSF) fname = nctk_ncify(dtfil%fnameabi_wfkfine)
2287 
2288      call wfk_read_eigenvalues(fname,eigenq_fine,hdr_fine,mpi_enreg%comm_world)
2289      ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband")
2290    end if
2291 ! DBSP ==> Has been changed to be able to make Bandstructure calculation
2292 !   if(dtset%kptopt==3 .or. dtset%kptopt==0)then
2293    if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%kptopt < -4 .or. dtset%nsym==1) then
2294 !END
2295      if (dtset%nsym > 1) then ! .and. dtset%efmas==0) then
2296        MSG_ERROR("Symmetries are not implemented for temperature dependence calculations")
2297      end if
2298      write(std_out,*) 'Entering: eig2stern'
2299      if(smdelta>0)then
2300        if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
2301 &       (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
2302          call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2303 &         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2304 &         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2305 &         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd,eigenq_fine,hdr_fine,hdr0)
2306        else
2307          call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2308 &         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2309 &         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2310 &         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd)
2311        end if
2312      else
2313        if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
2314 &       (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
2315          call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2316 &         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2317 &         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2318 &         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd,eigenq_fine,hdr_fine,hdr0)
2319        else
2320          call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2321 &         eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2322 &         dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2323 &         npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset)
2324        end if
2325      end if
2326      call wrtout(std_out, 'Leaving: eig2stern', "COLL")
2327 
2328      if (dtset%ieig2rf==1.or.dtset%ieig2rf==2) then
2329        if (me==master) then
2330 !        print _EIGR2D file for this perturbation
2331          dscrpt=' Note : temporary (transfer) database '
2332          unitout = dtfil%unddb
2333          vrsddb=100401
2334 
2335          call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,&
2336 &         1,xred=xred,occ=occ_pert)
2337 
2338          call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigr2d, dtfil%unddb)
2339 
2340          call outbsd(bdeigrf,dtset,eig2nkq,dtset%natom,nkpt_rbz,unitout)
2341 !        print _EIGI2D file for this perturbation
2342          if(smdelta>0) then
2343 
2344            unitout = dtfil%unddb
2345            call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigi2d, unitout)
2346 
2347            call outbsd(bdeigrf,dtset,eigbrd,dtset%natom,nkpt_rbz,unitout)
2348          end if !smdelta
2349 
2350          call ddb_hdr_free(ddb_hdr)
2351 
2352 !        Output of the EIGR2D.nc file.
2353          fname = strcat(dtfil%filnam_ds(4),"_EIGR2D.nc")
2354 
2355 !        Electronic band energies.
2356          bantot= dtset%mband*dtset%nkpt*dtset%nsppol
2357          call ebands_init(bantot,gkk_ebands,dtset%nelect,doccde,eigen0_pert,hdr0%istwfk,hdr0%kptns,&
2358 &         hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,&
2359 &         hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,&
2360 &         hdr0%charge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, &
2361 &         hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk)
2362 
2363 !        Second order derivative EIGR2D (real and Im)
2364          call eigr2d_init(eig2nkq,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
2365 #ifdef HAVE_NETCDF
2366          NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGR2D file")
2367          NCF_CHECK(crystal_ncwrite(crystal, ncid))
2368          NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid))
2369          call eigr2d_ncwrite(eigr2d,dtset%qptn(:),dtset%wtq,ncid)
2370          NCF_CHECK(nf90_close(ncid))
2371 #endif
2372          if(smdelta>0) then
2373 !          Output of the EIGI2D.nc file.
2374            fname = strcat(dtfil%filnam_ds(4),"_EIGI2D.nc")
2375 !          Broadening EIGI2D (real and Im)
2376            call eigr2d_init(eigbrd,eigi2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom)
2377 #ifdef HAVE_NETCDF
2378            NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGI2D file")
2379            NCF_CHECK(crystal_ncwrite(crystal, ncid))
2380            NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid))
2381            call eigr2d_ncwrite(eigi2d,dtset%qptn(:),dtset%wtq,ncid)
2382            NCF_CHECK(nf90_close(ncid))
2383 #endif
2384          end if
2385        end if
2386      end if !ieig2rf==1.or.ieig2rf==2
2387    else
2388      write(message,'(3a)')&
2389 &     'K point grids must be the same for every perturbation: eig2stern not called',ch10,&
2390 &     'Action: Put kptopt=3 '
2391      MSG_WARNING(message)
2392    end if !kptopt
2393    ABI_DEALLOCATE(gh1c_pert)
2394    ABI_DEALLOCATE(gh0c1_pert)
2395    ABI_DEALLOCATE(cg1_pert)
2396    ABI_DEALLOCATE(kpt_rbz_pert)
2397    ABI_DEALLOCATE(istwfk_pert)
2398    ABI_DEALLOCATE(npwarr_pert)
2399    ABI_DEALLOCATE(npwar1_pert)
2400    ABI_DEALLOCATE(npwtot_pert)
2401    ABI_DEALLOCATE(occ_pert)
2402  end if  !if dtset%ieig2rf
2403 
2404  !Calculation of effective masses.
2405  if(dtset%efmas == 1) then
2406    call efmas_main(cg0_pert,cg1_pert,dim_eig2rf,dtset,efmasdeg,efmasval,eigen0_pert,&
2407 &   eigen1_pert,gh0c1_pert,gh1c_pert,istwfk_pert,mpert,mpi_enreg,nkpt_rbz,npwarr_pert,rprimd)
2408    ABI_DEALLOCATE(gh1c_pert)
2409    ABI_DEALLOCATE(gh0c1_pert)
2410    ABI_DEALLOCATE(cg1_pert)
2411    ABI_DEALLOCATE(istwfk_pert)
2412    ABI_DEALLOCATE(npwarr_pert)
2413    ABI_DEALLOCATE(cg0_pert)
2414 
2415    if(dtset%prtefmas==1)then
2416      fname = strcat(dtfil%filnam_ds(4),"_EFMAS.nc")
2417      write(std_out,*)' dfpt_looppert : will write ',fname
2418 #ifdef HAVE_NETCDF
2419      NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EFMAS file")
2420      NCF_CHECK(crystal_ncwrite(crystal, ncid))
2421 !    NCF_CHECK(ebands_ncwrite(ebands_k, ncid)) ! At this stage, ebands_k is not available
2422      call print_efmas(efmasdeg,efmasval,kpt_rbz_pert,ncid)
2423      NCF_CHECK(nf90_close(ncid))
2424 #endif
2425    endif
2426 
2427    call efmas_analysis(dtset,efmasdeg,efmasval,kpt_rbz_pert,mpi_enreg,nkpt_rbz,rprimd)
2428 
2429    ABI_DEALLOCATE(kpt_rbz_pert)
2430 
2431  end if
2432 
2433 
2434 !Free memory.
2435  if(dtset%ieig2rf /= 3 .and. dtset%ieig2rf /= 4 .and. dtset%ieig2rf /= 5) then
2436    call hdr_free(hdr0)
2437  end if
2438  ABI_DEALLOCATE(eigen0_copy)
2439  call crystal_free(crystal)
2440 
2441  ! GKK stuff (deprecated)
2442  call ebands_free(gkk_ebands)
2443  call eigr2d_free(eigr2d)
2444  call eigr2d_free(eigi2d)
2445 
2446  call timab(147,2,tsec)
2447 !######################################################################################
2448 
2449 !Get ddk file information, for later use in dfpt_dyout
2450  ddkfil(:)=0
2451  do idir=1,3
2452    file_index(1)=idir+dtset%natom*3
2453    call appdig(file_index(1),dtfil%fnamewffddk,fiwfddk)
2454 !  Check that ddk file exists
2455    t_exist = file_exists(fiwfddk)
2456    if (.not. t_exist) then
2457      ! Trick needed to run Abinit test suite in netcdf mode.
2458      t_exist = file_exists(nctk_ncify(fiwfddk))
2459      if (t_exist) then
2460        write(message,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name."
2461        call wrtout(std_out,message,'COLL')
2462        fiwfddk = nctk_ncify(fiwfddk)
2463      end if
2464    end if
2465 
2466 !  If the file exists set ddkfil to a non-zero value
2467    if (t_exist) ddkfil(idir)=20+idir
2468  end do
2469 
2470  ABI_DEALLOCATE(ph1d)
2471  ABI_DEALLOCATE(ph1df)
2472  ABI_DEALLOCATE(pert_calc)
2473  if (psps%usepaw==1) then
2474    ABI_DEALLOCATE(dimcprj_srt)
2475  end if
2476 
2477 !destroy dtset_tmp
2478  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
2479    ABI_DEALLOCATE(rhor1_save)
2480    ABI_DEALLOCATE(blkflg_save)
2481    call dtset_free (dtset_tmp)
2482    ABI_DATATYPE_DEALLOCATE(dtset_tmp)
2483  end if
2484 
2485 !In paral_pert-case some array's have to be reconstructed
2486  if(mpi_enreg%paral_pert==1) then
2487    ABI_ALLOCATE(buffer1,(2,3,mpert,3,mpert*(2+psps%usepaw)))
2488    buffer1(:,:,:,:,1:mpert)=d2lo(:,:,:,:,:)
2489    buffer1(:,:,:,:,1+mpert:2*mpert)=d2nl(:,:,:,:,:)
2490    if (psps%usepaw==1) then
2491      buffer1(:,:,:,:,1+2*mpert:3*mpert)=d2ovl(:,:,:,:,:)
2492    end if
2493    call xmpi_sum(buffer1,mpi_enreg%comm_pert,ierr)
2494    call xmpi_sum(blkflg,mpi_enreg%comm_pert,ierr)
2495    if(dtset%prtbbb==1) then
2496      call xmpi_sum(d2bbb,mpi_enreg%comm_pert,ierr)
2497    end if
2498    d2lo(:,:,:,:,:)=buffer1(:,:,:,:,1:mpert)
2499    d2nl(:,:,:,:,:)=buffer1(:,:,:,:,1+mpert:2*mpert)
2500    if (psps%usepaw==1) then
2501      d2ovl(:,:,:,:,:)=buffer1(:,:,:,:,1+2*mpert:3*mpert)
2502    end if
2503    ABI_DEALLOCATE(buffer1)
2504  end if
2505 
2506  if ( associated(old_atmtab)) then
2507    ABI_DEALLOCATE(old_atmtab)
2508    nullify(old_atmtab)
2509  end if
2510 
2511  call status(0,dtfil%filstat,iexit,level,'exit')
2512 
2513  call timab(141,2,tsec)
2514 
2515  DBG_EXIT("COLL")
2516 
2517 end subroutine dfpt_looppert

ABINIT/dfpt_prtene [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_prtene

FUNCTION

 Print components of second derivative of total energy in nice format

INPUTS

 eberry=energy associated with Berry phase
 edocc=correction to 2nd-order total energy coming from changes of occupation
 eeig0=0th-order eigenenergies part of 2nd-order total energy
 eew=Ewald part of 2nd-order total energy
 efrhar=hartree frozen-wavefunction part of 2nd-order tot. en.
 efrkin=kinetic frozen-wavefunction part of 2nd-order tot. en.
 efrloc=local psp. frozen-wavefunction part of 2nd-order tot. en.
 efrnl=nonlocal psp. frozen-wavefunction part of 2nd-order tot. en
 efrx1=xc core corr.(1) frozen-wavefunction part of 2nd-order tot. en
 efrx2=xc core corr.(2) frozen-wavefunction part of 2nd-order tot. en
 ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
   for strain perturbation only (zero otherwise, and not used)
 ehart1=1st-order Hartree part of 2nd-order total energy
 eii=pseudopotential core part of 2nd-order total energy
 ek0=0th-order kinetic energy part of 2nd-order total energy.
 ek1=1st-order kinetic energy part of 2nd-order total energy.
 eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
 elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
 enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
 enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
 eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
       PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
 epaw1=1st-order PAW on-site part of 2nd-order total energy.
 evdw=DFT-D semi-empirical part of 2nd-order total energy
 exc1=1st-order exchange-correlation part of 2nd-order total energy
 iout=unit number to which output is written
 ipert=type of the perturbation
 natom=number of atoms in unit cell
 usepaw= 0 for non paw calculation; =1 for paw calculation
 usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use

OUTPUT

  (only writing)

NOTES

 all energies in Hartree

PARENTS

      dfpt_looppert

CHILDREN

      wrtout

SOURCE

2768 subroutine dfpt_prtene(berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
2769 &  ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,eovl1,epaw1,evdw,exc1,iout,ipert,natom,&
2770 &  usepaw,usevdw)
2771 
2772 
2773 !This section has been created automatically by the script Abilint (TD).
2774 !Do not modify the following lines by hand.
2775 #undef ABI_FUNC
2776 #define ABI_FUNC 'dfpt_prtene'
2777 !End of the abilint section
2778 
2779  implicit none
2780 
2781 !Arguments -------------------------------
2782 !scalars
2783  integer,intent(in) :: berryopt,iout,ipert,natom,usepaw,usevdw
2784  real(dp),intent(in) :: eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1
2785  real(dp),intent(in) :: efrx2,ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1
2786  real(dp),intent(in) :: eovl1,epaw1,evdw,exc1
2787 
2788 !Local variables -------------------------
2789 !scalars
2790  integer :: nn
2791  logical :: berry_activated
2792  real(dp) :: enl1_effective,erelax,etotal
2793  character(len=10) :: numb
2794  character(len=10),parameter :: numbstr(20) = &
2795 &  (/'One       ','Two       ','Three     ','Four      ','Five      ', &
2796 &    'Six       ','Seven     ','Eight     ','Nine      ','Ten       ', &
2797 &    'Eleven    ','Twelve    ','Thirteen  ','Fourteen  ','Fifteen   ', &
2798 &    'Sixteen   ','Seventeen ','Eighteen  ','Nineteen  ','Twenty    '/)
2799  character(len=500) :: message
2800 
2801 ! *********************************************************************
2802 
2803 !Count and print the number of components of 2nd-order energy
2804 !MT feb 2015: this number is wrong! Should change it but
2805 !             need to change a lot of ref. files
2806  berry_activated=(berryopt== 4.or.berryopt== 6.or.berryopt== 7.or. &
2807 & berryopt==14.or.berryopt==16.or.berryopt==17)
2808  if (ipert==natom+1) nn=8
2809  if (ipert==natom+5) nn=8
2810  if (ipert==natom+2) nn=7
2811  if (ipert>=1.and.ipert<=natom) nn=13
2812  if (ipert==natom+3.or.ipert==natom+4) nn=17
2813  if (ipert==natom+2.and.berry_activated) nn=nn+1
2814  if (ipert==natom+10.or.ipert==natom+11) nn=1 ! means nothing,
2815 ! because we do not compute derivatives of the energy in this case
2816  if (usepaw==1) nn=nn+1
2817  if (usevdw==1) nn=nn+1
2818  write(message, '(4a)' ) ch10,&
2819 & ' ',trim(numbstr(nn)),' components of 2nd-order total energy (hartree) are '
2820  call wrtout(iout,message,'COLL')
2821  call wrtout(std_out,message,'COLL')
2822 
2823  numb='1,2,3'
2824  write(message, '(3a)' )&
2825 & ' ',trim(numb),': 0th-order hamiltonian combined with 1st-order wavefunctions'
2826  call wrtout(iout,message,'COLL')
2827  call wrtout(std_out,message,'COLL')
2828  write(message, '(a,es17.8,a,es17.8,a,es17.8)' )&
2829 & '     kin0=',ek0,   ' eigvalue=',eeig0,'  local=',eloc0
2830  call wrtout(iout,message,'COLL')
2831  call wrtout(std_out,message,'COLL')
2832 
2833  numb='4,5,6';if( ipert==natom+3.or.ipert==natom+4) numb='4,5,6,7'
2834  write(message, '(3a)' )&
2835 & ' ',trim(numb),': 1st-order hamiltonian combined with 1st and 0th-order wfs'
2836  call wrtout(iout,message,'COLL')
2837  call wrtout(std_out,message,'COLL')
2838  if(ipert/=natom+1.and.ipert/=natom+2)then
2839    write(message, '(a,es17.8,a,es17.8,a,es17.8,a,a)' ) &
2840 &   ' loc psp =',elpsp1,'  Hartree=',ehart1,'     xc=',exc1,ch10,&
2841 &   ' note that "loc psp" includes a xc core correction that could be resolved'
2842  else if(ipert==natom+1) then
2843    write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2844 &   '     kin1=',ek1,   '  Hartree=',ehart1,'     xc=',exc1
2845  else if(ipert==natom+2) then
2846    write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2847 &   '    dotwf=',enl1,  '  Hartree=',ehart1,'     xc=',exc1
2848  end if
2849  if(ipert==natom+3 .or. ipert==natom+4) then
2850    write(message, '(a,es17.8,a,es17.8,a,es17.8,a,a,es17.8)' ) &
2851 &   ' loc psp =',elpsp1,'  Hartree=',ehart1,'     xc=',exc1,ch10,&
2852 &   '     kin1=',ek1
2853  end if
2854  call wrtout(iout,message,'COLL')
2855  call wrtout(std_out,message,'COLL')
2856 
2857  enl1_effective=enl1;if (ipert==natom+2) enl1_effective=zero
2858  numb='7,8,9';if( ipert==natom+3.or.ipert==natom+4) numb='8,9,10'
2859  write(message, '(5a,es17.8,a,es17.8,a,es17.8)' )&
2860 & ' ',trim(numb),': eventually, occupation + non-local contributions',ch10,&
2861 & '    edocc=',edocc,'     enl0=',enl0,'   enl1=',enl1_effective
2862  call wrtout(iout,message,'COLL')
2863  call wrtout(std_out,message,'COLL')
2864 
2865  if (usepaw==1) then
2866    numb='10';if( ipert==natom+3.or.ipert==natom+4) numb='11'
2867    write(message, '(3a,es17.8)' )&
2868 &   ' ',trim(numb),': eventually, PAW "on-site" Hxc contribution: epaw1=',epaw1
2869    call wrtout(iout,message,'COLL')
2870    call wrtout(std_out,message,'COLL')
2871  end if
2872 
2873  if(ipert/=natom+10 .and.ipert/=natom+11) then
2874    erelax=0.0_dp
2875    if(ipert>=1.and.ipert<=natom)then
2876      erelax=ek0+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1
2877    else if(ipert==natom+1.or.ipert==natom+2)then
2878      erelax=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+epaw1
2879    else if(ipert==natom+3.or.ipert==natom+4)then
2880      erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1
2881    else if(ipert==natom+5)then
2882      erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1
2883    end if
2884    enl1_effective=enl1
2885    if (ipert==natom+1.or.ipert==natom+2) then
2886      if (1.0_dp+enl1/10.0_dp==1.0_dp) enl1_effective=zero
2887    end if
2888 
2889    numb='1-9';if (usepaw==1) numb='1-10'
2890    if( ipert==natom+3.or.ipert==natom+4) then
2891      numb='1-10';if (usepaw==1) numb='1-11'
2892    end if
2893    write(message, '(5a,es17.8)' )&
2894 &   ' ',trim(numb),' gives the relaxation energy (to be shifted if some occ is /=2.0)',&
2895 &   ch10,'   erelax=',erelax
2896    call wrtout(iout,message,'COLL')
2897    call wrtout(std_out,message,'COLL')
2898  end if
2899 
2900  if(ipert>=1.and.ipert<=natom)then
2901 
2902    numb='10,11,12';if (usepaw==1) numb='11,12,13'
2903    write(message, '(4a)' )&
2904 &   ' ',trim(numb),' Non-relaxation  contributions : ',&
2905 &   'frozen-wavefunctions and Ewald'
2906    call wrtout(iout,message,'COLL')
2907    call wrtout(std_out,message,'COLL')
2908    write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2909 &   ' fr.local=',efrloc,' fr.nonlo=',efrnl,'  Ewald=',eew
2910    call wrtout(iout,message,'COLL')
2911    call wrtout(std_out,message,'COLL')
2912 
2913    write(message, '(a,es16.6)' )' dfpt_prtene : non-relax=',efrloc+efrnl+eew
2914    call wrtout(std_out,message,'COLL')
2915 
2916    numb='13,14';if (usepaw==1) numb='14,15'
2917    write(message, '(3a)' )&
2918 &   ' ',trim(numb),' Frozen wf xc core corrections (1) and (2)'
2919    call wrtout(iout,message,'COLL')
2920    call wrtout(std_out,message,'COLL')
2921    write(message, '(a,es17.8,a,es17.8)' ) &
2922 &   ' frxc 1  =',efrx1,'  frxc 2 =',efrx2
2923    call wrtout(iout,message,'COLL')
2924    call wrtout(std_out,message,'COLL')
2925    if (usepaw==1) then
2926      numb='16'
2927      write(message, '(5a,es17.8)' )&
2928 &     ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
2929 &     ch10,' eovl1 =',eovl1
2930      call wrtout(iout,message,'COLL')
2931      call wrtout(std_out,message,'COLL')
2932    end if
2933    if (usevdw==1) then
2934      numb='15';if (usepaw==1) numb='17'
2935      write(message, '(3a,es17.8)' )&
2936 &     ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw
2937      call wrtout(iout,message,'COLL')
2938      call wrtout(std_out,message,'COLL')
2939    end if
2940 
2941    write(message, '(a)' )' Resulting in : '
2942    call wrtout(iout,message,'COLL')
2943    call wrtout(std_out,message,'COLL')
2944    etotal=erelax+eew+efrloc+efrnl+efrx1+efrx2+evdw
2945    write(message, '(a,e20.10,a,e22.12,a)' ) &
2946 &   ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
2947    call wrtout(iout,message,'COLL')
2948    call wrtout(std_out,message,'COLL')
2949    write(message, '(a,es20.10,a,es20.10,a)' ) &
2950 &   '    (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)'
2951    call wrtout(iout,message,'COLL')
2952    call wrtout(std_out,message,'COLL')
2953    write(message, '(a,es20.10,a,a)' ) &
2954 &   '    (  non-var. 2DEtotal :',&
2955 &   0.5_dp*(elpsp1+enl1)+eovl1+eew+efrloc+efrnl+efrx1+efrx2+evdw,' Ha)',ch10
2956    call wrtout(iout,message,'COLL')
2957    call wrtout(std_out,message,'COLL')
2958 
2959  else if(ipert==natom+1.or.ipert==natom+2)then
2960    if (ipert==natom+1.and.usepaw==1) then
2961      numb='11'
2962      write(message, '(5a,es17.8)' )&
2963 &     ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
2964 &     ch10,' eovl1 =',eovl1
2965      call wrtout(iout,message,'COLL')
2966      call wrtout(std_out,message,'COLL')
2967    end if
2968    write(message,*)' No Ewald or frozen-wf contrib.:',&
2969 &   ' the relaxation energy is the total one'
2970    if(berry_activated) then
2971      numb='10';
2972      write(message,'(3a,es20.10)')' ',trim(numb),' Berry phase energy :',eberry
2973    end if
2974    call wrtout(iout,message,'COLL')
2975    call wrtout(std_out,message,'COLL')
2976    etotal=erelax
2977    write(message, '(a,e20.10,a,e22.12,a)' ) &
2978 &   ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
2979    call wrtout(iout,message,'COLL')
2980    call wrtout(std_out,message,'COLL')
2981    write(message, '(a,es20.10,a)' ) &
2982 &   '    (  non-var. 2DEtotal :',0.5_dp*(ek1+enl1_effective)+eovl1,' Ha)'
2983    call wrtout(iout,message,'COLL')
2984    call wrtout(std_out,message,'COLL')
2985 
2986  else if(ipert==natom+3 .or. ipert==natom+4) then
2987    numb='11,12,13';if (usepaw==1) numb='12,13,14'
2988    write(message, '(4a)' )&
2989 &   ' ',trim(numb),' Non-relaxation  contributions : ',&
2990 &   'frozen-wavefunctions and Ewald'
2991    call wrtout(iout,message,'COLL')
2992    call wrtout(std_out,message,'COLL')
2993    write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2994 &   '  fr.hart=',efrhar,'   fr.kin=',efrkin,' fr.loc=',efrloc
2995    call wrtout(iout,message,'COLL')
2996    call wrtout(std_out,message,'COLL')
2997 
2998    numb='14,15,16';if (usepaw==1) numb='15,16,17'
2999    write(message, '(4a)' )&
3000 &   ' ',trim(numb),' Non-relaxation  contributions : ',&
3001 &   'frozen-wavefunctions and Ewald'
3002    call wrtout(iout,message,'COLL')
3003    call wrtout(std_out,message,'COLL')
3004    write(message, '(a,es17.8,a,es17.8,a,es17.8)' ) &
3005 &   '  fr.nonl=',efrnl,'    fr.xc=',efrx1,'  Ewald=',eew
3006    call wrtout(iout,message,'COLL')
3007    call wrtout(std_out,message,'COLL')
3008 
3009    numb='17';if (usepaw==1) numb='18'
3010    write(message, '(4a)' )&
3011 &   ' ',trim(numb),' Non-relaxation  contributions : ',&
3012 &   'pseudopotential core energy'
3013    call wrtout(iout,message,'COLL')
3014    call wrtout(std_out,message,'COLL')
3015    write(message, '(a,es17.8)' ) &
3016 &   '  pspcore=',eii
3017    call wrtout(iout,message,'COLL')
3018    call wrtout(std_out,message,'COLL')
3019    if (usepaw==1) then
3020      numb='19'
3021      write(message, '(5a,es17.8)' )&
3022 &     ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
3023 &     ch10,' eovl1 =',eovl1
3024      call wrtout(iout,message,'COLL')
3025      call wrtout(std_out,message,'COLL')
3026    end if
3027    if (usevdw==1) then
3028      numb='18';if (usepaw==1) numb='20'
3029      write(message, '(3a,es17.8)' )&
3030 &     ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw
3031      call wrtout(iout,message,'COLL')
3032      call wrtout(std_out,message,'COLL')
3033    end if
3034 
3035    write(message, '(a,es16.6)' )' dfpt_prtene : non-relax=',&
3036 &   efrhar+efrkin+efrloc+efrnl+efrx1+eew+evdw
3037    call wrtout(std_out,message,'COLL')
3038    write(message, '(a)' )' Resulting in : '
3039    call wrtout(iout,message,'COLL')
3040    call wrtout(std_out,message,'COLL')
3041    etotal=erelax+efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw
3042    write(message, '(a,e20.10,a,e22.12,a)' ) &
3043 &   ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
3044    call wrtout(iout,message,'COLL')
3045    call wrtout(std_out,message,'COLL')
3046    write(message, '(a,es20.10,a,es20.10,a)' ) &
3047 &   '    (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)'
3048    call wrtout(iout,message,'COLL')
3049    call wrtout(std_out,message,'COLL')
3050    write(message, '(a,es20.10,a,a)' ) &
3051 &   '    (  non-var. 2DEtotal :',&
3052 &   0.5_dp*(elpsp1+enl1+ek1+ehart01)+eovl1+&
3053 &   efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw,' Ha)',ch10
3054    call wrtout(iout,message,'COLL')
3055    call wrtout(std_out,message,'COLL')
3056  end if
3057 
3058 end subroutine dfpt_prtene

ABINIT/eigen_meandege [ Functions ]

[ Top ] [ Functions ]

NAME

 eigen_meandege

FUNCTION

 This routine takes the mean values of the responses
 for the eigenstates that are degenerate in energy.

INPUTS

  eigenresp((3-option)*mband**(3-option)*nkpt*nsppol)= input eigenresp
       eigenrep(2*mband**2*nkpt*nsppol) for first-order derivatives of the eigenvalues
       eigenrep(mband*nkpt*nsppol) for Fan or Debye-Waller second-order derivatives of the eigenvalues
  mband= maximum number of bands
  natom= number of atoms in the unit cell
  nkpt= number of k-points
  nsppol= 1 for unpolarized, 2 for spin-polarized
  option= 1 for eigen(1), 2 for eigen(2) - Fan or Debye-Waller

OUTPUT

  eigenresp_mean(mband*nkpt*nsppol)= eigenresp, averaged over degenerate states

PARENTS

      dfpt_looppert,respfn

CHILDREN

SOURCE

3089 subroutine eigen_meandege(eigen0,eigenresp,eigenresp_mean,mband,nband,nkpt,nsppol,option)
3090 
3091 
3092 !This section has been created automatically by the script Abilint (TD).
3093 !Do not modify the following lines by hand.
3094 #undef ABI_FUNC
3095 #define ABI_FUNC 'eigen_meandege'
3096 !End of the abilint section
3097 
3098  implicit none
3099 
3100 !Arguments ------------------------------------
3101 !scalars
3102  integer,intent(in) :: mband,nkpt,nsppol,option
3103  integer,intent(in) :: nband(nkpt*nsppol)
3104 
3105 !arrays
3106  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
3107  real(dp),intent(in) :: eigenresp((3-option)*mband**(3-option)*nkpt*nsppol)
3108  real(dp),intent(out) :: eigenresp_mean(mband*nkpt*nsppol)
3109 
3110 !Local variables-------------------------------
3111 !scalars
3112  integer :: bdtot_index,bd2tot_index,iband,ii,ikpt,isppol,nband_k
3113  real(dp) :: eig0,mean
3114  character(len=500) :: message
3115 !arrays
3116 
3117 ! *********************************************************************
3118 
3119  if(option/=1 .and. option/=2)then
3120    write(message, '(a,i0)' )' The argument option should be 1 or 2, while it is found that option=',option
3121    MSG_BUG(message)
3122  end if
3123 
3124  bdtot_index=0 ; bd2tot_index=0
3125  do isppol=1,nsppol
3126    do ikpt=1,nkpt
3127      nband_k=nband(ikpt+(isppol-1)*nkpt)
3128      if(option==1)then
3129        do iband=1,nband_k
3130          eigenresp_mean(iband+bdtot_index)=&
3131 &         eigenresp(2*iband-1 + (iband-1)*2*nband_k + bd2tot_index)
3132        end do
3133      else if(option==2)then
3134        do iband=1,nband_k
3135          eigenresp_mean(iband+bdtot_index)=eigenresp(iband+bdtot_index)
3136        end do
3137      end if
3138 
3139 !    Treat the case of degeneracies : take the mean of degenerate states
3140      if(nband_k>1)then
3141        eig0=eigen0(1+bdtot_index)
3142        ii=1
3143        do iband=2,nband_k
3144          if(eigen0(iband+bdtot_index)-eig0<tol8)then
3145            ii=ii+1
3146          else
3147            mean=sum(eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index))/ii
3148            eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index)=mean
3149            ii=1
3150          end if
3151          eig0=eigen0(iband+bdtot_index)
3152          if(iband==nband_k)then
3153            mean=sum(eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index))/ii
3154            eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index)=mean
3155          end if
3156        end do
3157      end if
3158 
3159      bdtot_index=bdtot_index+nband_k
3160      bd2tot_index=bd2tot_index+2*nband_k**2
3161    end do
3162  end do
3163 
3164 end subroutine eigen_meandege

ABINIT/getcgqphase [ Functions ]

[ Top ] [ Functions ]

NAME

 getcgqphase

FUNCTION

 extract phases from wave functions, to cancel contributions to gkk matrix elements

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  timrev = flag for use of time reversal symmetry
  cg = input wavefunctions
  mcg = dimension of cg = nspinor*mband*mpw*mkmem
  cgq = input wavefunctions at k+q
  mcgq = dimension of cgq = nspinor*mband*mpw*mkmem
  mpi_enreg = datastructure for mpi communication
  nkpt_rbz = number of k-points in reduced zone for present q point
  npwarr = array of numbers of plane waves for each k-point
  npwar1 = array of numbers of plane waves for each k+q point

OUTPUT

  phasecg = phase of different wavefunction products <k,n | k+q,n'>

PARENTS

      dfpt_looppert

CHILDREN

      smatrix,wrtout,xmpi_barrier,xmpi_sum_master

SOURCE

2550 subroutine getcgqphase(dtset, timrev, cg,  mcg,  cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg)
2551 
2552 
2553 !This section has been created automatically by the script Abilint (TD).
2554 !Do not modify the following lines by hand.
2555 #undef ABI_FUNC
2556 #define ABI_FUNC 'getcgqphase'
2557 !End of the abilint section
2558 
2559  implicit none
2560 
2561 !Arguments -------------------------------
2562  ! scalars
2563  integer, intent(in) :: mcg, mcgq, timrev
2564  integer, intent(in) :: nkpt_rbz
2565  type(dataset_type), intent(in) :: dtset
2566 
2567  ! arrays
2568  integer, intent(in) :: npwarr(nkpt_rbz)
2569  integer, intent(in) :: npwar1(nkpt_rbz)
2570  real(dp), intent(in) :: cg(2,mcg)
2571  real(dp), intent(in) :: cgq(2,mcgq)
2572  type(MPI_type), intent(in) :: mpi_enreg
2573  real(dp),intent(out) :: phasecg(2, dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)
2574 
2575 !Local variables -------------------------
2576  ! local vars
2577  integer :: icg, icgq, isppol, ikpt, ipw
2578  integer :: istate, iband1, iband2
2579  integer :: npw_k, npw_q
2580  integer :: me, ierr, master, spaceComm, nprocs
2581  integer :: usepaw
2582  integer :: ddkflag, itrs, job, maxbd, mcg1_k, minbd, shiftbd
2583  real(dp) :: normsmat
2584 
2585  integer, allocatable :: sflag_k(:)
2586  integer, allocatable :: pwind_k(:)
2587 
2588  real(dp) :: cg1_dummy(1,1)
2589  real(dp) :: smat_inv_dummy(1,1,1)
2590  real(dp) :: smat_k_paw_dummy(1,1,1)
2591  real(dp) :: dtm_k_dummy(2)
2592  real(dp), allocatable :: smat_k(:,:,:)
2593  real(dp), allocatable :: pwnsfac_k(:,:)
2594  logical, allocatable :: my_kpt(:,:)
2595  character(len=500) :: message
2596 
2597 ! *********************************************************************
2598 
2599  ABI_ALLOCATE(smat_k,(2,dtset%mband,dtset%mband))
2600  ABI_ALLOCATE(sflag_k,(dtset%mband))
2601 
2602 !dummy use of timrev so abirules stops complaining.
2603  icg = timrev
2604 
2605 !!MPI data for future use
2606  spaceComm=mpi_enreg%comm_cell
2607  nprocs=xmpi_comm_size(spaceComm)
2608  master=0
2609  me=mpi_enreg%me_kpt
2610 
2611  ABI_ALLOCATE(my_kpt, (nkpt_rbz, dtset%nsppol))
2612  my_kpt = .true.
2613  if (mpi_enreg%nproc_kpt > 1) then
2614    do isppol = 1, dtset%nsppol
2615      do ikpt = 1, nkpt_rbz
2616        my_kpt(ikpt, isppol) = .not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,&
2617 &       dtset%nband(ikpt),isppol,me))
2618      end do
2619    end do
2620  end if
2621 
2622 
2623 !make trivial association of G vectors: we just want <psi_k| psi_k+q>
2624 !TODO: check this is correct wrt arrangement of kg vectors for k+q
2625 !looks ok : usually made in initberry, from the translations associated
2626 !to the symops, scalar product with the G vectors. The symop is the one
2627 !used to go from the irreducible k to the full zone k. In present context
2628 !we should be using only the reduced zone, and anyhow have the same k-grid
2629 !for the gkk matrix elements and for the cg here...
2630  ABI_ALLOCATE(pwind_k,(dtset%mpw))
2631  ABI_ALLOCATE(pwnsfac_k,(4,dtset%mpw))
2632  do ipw = 1, dtset%mpw
2633    pwind_k(ipw) = ipw
2634    pwnsfac_k(1,ipw) = one
2635    pwnsfac_k(2,ipw) = zero
2636    pwnsfac_k(3,ipw) = one
2637    pwnsfac_k(4,ipw) = zero
2638  end do
2639 
2640 !flags for call to smatrix
2641  usepaw = 0 ! for now
2642  ddkflag = 0
2643  itrs = 0
2644  job = 0
2645  maxbd = 1
2646  mcg1_k = 1
2647  minbd = 1
2648  shiftbd = 1
2649 
2650 !from overlap matrix for each wavefunction, extract phase
2651  icg = 0
2652  icgq = 0
2653  istate = 0
2654 
2655  phasecg = zero
2656  do isppol = 1, dtset%nsppol
2657    do ikpt = 1, nkpt_rbz
2658 !    each proc only has certain k
2659      if (.not. my_kpt(ikpt, isppol)) then
2660        istate = istate +  dtset%nband(ikpt)*dtset%nband(ikpt)
2661        cycle
2662      end if
2663 
2664      npw_k = npwarr(ikpt)
2665      npw_q= npwar1(ikpt)
2666 
2667 !    TODO: question: are the k-points in the ibz correctly ordered in cg and cgq? if not the icg below have to be adapted.
2668      sflag_k = 0 ! make sure all elements are calculated
2669      smat_k = zero
2670 
2671      call smatrix(cg, cgq, cg1_dummy, ddkflag, dtm_k_dummy, icg, icgq,&
2672 &     itrs, job, maxbd, mcg, mcgq, mcg1_k, minbd,dtset%mpw, dtset%mband, dtset%mband,&
2673 &     npw_k, npw_q, dtset%nspinor, pwind_k, pwnsfac_k, sflag_k, shiftbd,&
2674 &     smat_inv_dummy, smat_k, smat_k_paw_dummy, usepaw)
2675 
2676      icg  = icg  + npw_k*dtset%nspinor*dtset%nband(ikpt)
2677      icgq = icgq + npw_q*dtset%nspinor*dtset%nband(ikpt)
2678 
2679      do iband1 = 1, dtset%nband(ikpt)
2680        do iband2 = 1, dtset%nband(ikpt)
2681          istate = istate + 1
2682 !        normalise the overlap matrix element to get just the phase difference phi_k - phi_k+q
2683          normsmat = sqrt(smat_k(1,iband2, iband1)**2 &
2684 &         + smat_k(2,iband2, iband1)**2)
2685          if (normsmat > tol12) then
2686            phasecg(:,istate) = smat_k(:,iband2, iband1) / normsmat
2687 !          NOTE: 21/9/2011 these appear to be always 1, i, or -i, to within 1.e-5 at worst!
2688          end if
2689        end do
2690      end do
2691    end do
2692  end do
2693 
2694 !eventually do an mpi allreduce over the k-points for phasecg
2695  if (nprocs>1) then
2696    call xmpi_barrier(spaceComm)
2697    call xmpi_sum_master(phasecg,master,spaceComm,ierr)
2698    call xmpi_barrier(spaceComm)
2699    if (1==1) then
2700      write(message,'(a)') '  In getcgqphase - contributions to phasecg collected'
2701      call wrtout(std_out,message,'PERS')
2702    end if
2703  end if
2704 
2705  ABI_DEALLOCATE(sflag_k)
2706  ABI_DEALLOCATE(smat_k)
2707  ABI_DEALLOCATE(pwind_k)
2708  ABI_DEALLOCATE(pwnsfac_k)
2709  ABI_DEALLOCATE(my_kpt)
2710 
2711 end subroutine getcgqphase

ABINIT/m_dfpt_loopert [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_loopert

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XG, DRH, MB, XW, MT, SPr, MJV)
  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

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_dfpt_loopert
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_efmas_defs
34  use m_abicore
35  use m_xmpi
36  use m_errors
37  use m_wfk
38  use m_wffile
39  use m_io_redirect
40  use m_paral_pert
41  use m_abi_etsf
42  use m_nctk
43  use m_ddb
44 #ifdef HAVE_NETCDF
45  use netcdf
46 #endif
47  use m_hdr
48  use m_ebands
49 
50  use m_dtfil,      only : status
51  use m_occ,        only : getnel
52  use m_ddb_hdr,    only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write
53  use m_io_tools,   only : file_exists
54  use m_time,       only : timab
55  use m_fstrings,   only : strcat
56  use m_geometry,   only : mkrdim, metric, littlegroup_pert
57  use m_exit,       only : exit_check, disable_timelimit
58  use m_atomdata,   only : atom_gauss
59  use m_eig2d,      only : eigr2d_init,eigr2d_t, eigr2d_ncwrite,eigr2d_free, &
60                           gkk_t, gkk_init, gkk_ncwrite,gkk_free, outbsd, eig2stern
61  use m_crystal,    only : crystal_init, crystal_free, crystal_t
62  use m_crystal_io, only : crystal_ncwrite
63  use m_efmas,      only : efmas_main, efmas_analysis, print_efmas
64  use m_fft,        only : fourdp
65  use m_kg,         only : getcut, getmpw, kpgio, getph
66  use m_dtset,      only : dtset_copy, dtset_free
67  use m_iowf,       only : outwf
68  use m_ioarr,      only : read_rhor
69  use m_pawang,     only : pawang_type, pawang_init, pawang_free
70  use m_pawrad,     only : pawrad_type
71  use m_pawtab,     only : pawtab_type
72  use m_paw_an,     only : paw_an_type
73  use m_paw_ij,     only : paw_ij_type
74  use m_pawfgrtab,  only : pawfgrtab_type
75  use m_pawrhoij,   only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_bcast, pawrhoij_copy, &
76                           pawrhoij_nullify, pawrhoij_redistribute, pawrhoij_get_nspden
77  use m_pawcprj,    only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, pawcprj_getdim
78  use m_pawfgr,     only : pawfgr_type
79  use m_paw_sphharm,only : setsym_ylm
80  use m_rf2,        only : rf2_getidirs
81  use m_iogkk,      only : outgkk
82  use m_inwffil,    only : inwffil
83  use m_spacepar,   only : rotate_rho, setsym
84  use m_initylmg,   only : initylmg
85  use m_dfpt_scfcv, only : dfpt_scfcv
86  use m_dfpt_mkrho, only : dfpt_mkrho
87  use m_mpinfo,     only : initmpi_band, distrb2, proc_distrb_cycle
88  use m_atm2fft,    only : dfpt_atm2fft
89  use m_berrytk,    only : smatrix
90  use m_common,     only : prteigrs
91  use m_fourier_interpol, only : transgrid
92  use m_mkcore,     only : dfpt_mkcore
93  use m_mklocl,     only : dfpt_vlocal, vlocalstr
94  use m_cgprj,      only : ctocprj
95  use m_symkpt,     only : symkpt
96 
97  implicit none
98 
99  private