Logo Search packages:      
Sourcecode: abinit version File versions  Download package

loper3.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/loper3
!! NAME
!! loper3
!!
!! FUNCTION
!! Loop over perturbations
!!
!! COPYRIGHT
!! Copyright (C) 1999-2007 ABINIT group (XG, DRH, MB, XW)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
!!
!! INPUTS
!!  amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...)
!!  atindx(natom)=index table for atoms (see scfcv.f)
!!  atindx1(natom)=index table for atoms, inverse of atindx (see scfcv.f)
!!  codvsn=code version
!!  cpui=initial cpu time
!!  cpus=cpu time limit in seconds
!!  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(3,3,natom)=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
!!  eltcore(6,6)=core contribution to the elastic tensor
!!  elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor
!!  eltfrhar(6,6)=Hartree contribution to the elastic tensor
!!  eltfrkin(6,6)=kinetic contribution to the elastic tensor
!!  eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor
!!  eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor
!!  eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor
!!  fermie=fermi energy (Hartree)
!!  gsqcut_eff=Fourier cutoff on G^2 for "large sphere" of radius double
!!    that of the basis sphere--appropriate for charge density rho(G),
!!    Hartree potential, and pseudopotentials, corresponding to ecut_eff
!!  iexit=index of "exit" on first line of file (0 if not found)
!!  indsym(4,nsym,natom)=indirect indexing array for atom labels
!!  kxc(nfft,nkxc)=exchange and correlation kernel (see rhohxc.f)
!!  mgfft=maximum single fft dimension
!!  mkmem =max. number of k points which can fit in memory (GS data)  ; 0 if use disk
!!  mkqmem=max. number of k+q points which can fit in memory (GS data); 0 if use disk
!!  mk1mem=max. number of k points which can fit in memory (RF data)  ; 0 if use disk
!!  mpert=maximum number of ipert
!!  mpi_enreg=informations about MPI parallelization
!!  mpsang=1+maximum angular momentum for nonlocal pseudopotential
!!  nattyp(ntypat)= # atoms of each type.
!!  nfft=(effective) number of FFT grid points (for this processor)
!!  nkpt=number of k points
!!  nkxc=second dimension of the kxc array
!!  nspden=number of spin-density components
!!  nspinor=number of spinorial components of the wavefunctions
!!  nsym=number of symmetry elements in space group
!!  occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
!!  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,nfft)=array for Fourier transform of GS electron density
!!  rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
!!  symq(4,2,nsym)=1 if symmetry preserves present qpoint. From symq3
!!  symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
!!  timrev=1 if time-reversal preserves the q wavevector; 0 otherwise.
!!  tmpfil(11)=names for the temporary files based on dtfil%filnam_ds(5)
!!  vtrial(nfft,nspden)=GS potential (Hartree)
!!  vxcavg=average of vxc potential
!!  walli=initial wall clock time
!!  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
!!  etotal=total energy (sum of 7 contributions) (hartree)
!!
!! PARENTS
!!      respfn
!!
!! CHILDREN
!!      appdig,bstruct_clean,bstruct_init,chkexi,distrb2,fourdp,getmpw,getnel
!!      getph,handle_ncerr,hdr_clean,hdr_init,hdr_update,initmpi_respfn
!!      initylmg,insy3,inwffil,ioarr,kpgio,leave_new,leave_test,metric,mkcor3
!!      mkrdim,mkrho3,mpi_bcast,mpi_comm_free,outgkk,outwf,prteigrs,prtene3
!!      scfcv3,setsym,status,symkpt,timab,vloca3,vlocalstr,wffclose,wffdelete
!!      wffopen,wrtout,xcomm_init,xme_init
!!
!! SOURCE

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

subroutine loper3(amass,atindx,atindx1,blkflg,codvsn,cpui,cpus,doccde,&
&  ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,d2bbb,d2lo,d2nl,&
&  eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,&
&  etotal,fermie,gsqcut_eff,iexit,indsym,kxc,mgfft,&
&  mkmem,mkqmem,mk1mem,mpert,mpi_enreg,mpsang,nattyp,&
&  nfft,nkpt,nkxc,nspden,nspinor,nsym,occ,pawtab,pertsy,prtbbb,&
&  psps,rfpert,rhog,rhor,symq,symrec,timrev,tmpfil,&
&  vtrial,vxcavg,walli,xred)

 use defs_basis
 use defs_datatypes
#if defined HAVE_NETCDF
 use netcdf
#endif

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_01manage_mpi
 use interfaces_11util
 use interfaces_12ffts
 use interfaces_12geometry
 use interfaces_13io_mpi
 use interfaces_13recipspace
 use interfaces_14iowfdenpot
 use interfaces_14occeig
 use interfaces_15common
 use interfaces_16response
 use interfaces_18seqpar, except_this_one => loper3
 use interfaces_lib01hidempi
#endif
!End of the abilint section

 implicit none

#if defined MPI
           include 'mpif.h'
#endif
!Arguments ------------------------------------
 integer, intent(in) :: mgfft,mk1mem,mkmem,mkqmem,mpert,mpsang,nfft,nkpt
 integer, intent(in) :: nkxc,nspden,nsym,prtbbb,timrev
 integer, intent(inout) :: nspinor
 integer, intent(out) :: iexit
 real(dp), intent(in) :: cpui,cpus,gsqcut_eff,vxcavg,walli
 real(dp), intent(inout) :: fermie
 real(dp), intent(out) :: etotal
 character(len=6), intent(in) :: codvsn
 type(MPI_type), intent(inout) :: mpi_enreg
 type(datafiles_type), intent(in) :: dtfil
 type(dataset_type), intent(inout) :: dtset
 type(pseudopotential_type), intent(in) :: psps
 integer, intent(in) :: atindx(dtset%natom),atindx1(dtset%natom)
 integer, intent(in) :: indsym(4,nsym,dtset%natom)
 integer, intent(in) :: nattyp(dtset%ntypat),pertsy(3,mpert)
 integer, intent(in) :: rfpert(mpert),symq(4,2,nsym),symrec(3,3,nsym)
 integer, intent(out) :: blkflg(3,mpert,3,mpert),ddkfil(3)
 real(dp), intent(in) :: amass(dtset%natom),doccde(dtset%mband*nkpt*dtset%nsppol)
 real(dp), intent(in) :: dyew(2,3,dtset%natom,3,dtset%natom)
 real(dp), intent(in) :: dyfrlo(3,3,dtset%natom),dyfrnl(3,3,dtset%natom)
 real(dp), intent(in) :: dyfrx1(2,3,dtset%natom,3,dtset%natom)
 real(dp), intent(in) :: dyfrx2(3,3,dtset%natom),eltcore(6,6)
 real(dp), intent(in) :: elteew(6+3*dtset%natom,6),eltfrhar(6,6)
 real(dp), intent(in) :: eltfrkin(6,6),eltfrloc(6+3*dtset%natom,6)
 real(dp), intent(in) :: eltfrnl(6+3*dtset%natom,6)
 real(dp), intent(in) :: eltfrxc(6+3*dtset%natom,6),kxc(nfft,nkxc)
 real(dp), intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol)
 real(dp), intent(in) :: rhog(2,nfft),rhor(nfft,nspden)
 real(dp), intent(inout) :: vtrial(nfft,nspden),xred(3,dtset%natom)
 real(dp), intent(out) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)
 real(dp), intent(out) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert)
 character(len=fnlen), intent(in) :: tmpfil(12)
 type(pawtab_type), intent(in) :: pawtab(psps%ntypat*psps%usepaw)

!Local variables-------------------------------
!Define file format for different type of files. Presently,
!only one file format is supported for each type of files, but this might
!change soon ...
!1   for wavefunction file, old format (version prior to 2.0)
!2   for wavefunction file, new format (version 2.0 and after)    (fform)
!(51 or 52   for density rho(r)       (fformr)
!101 or 102 for potential V(r) file. (fformv)
 integer,parameter :: fform=2,fformv=102,level=14,response=1
 integer :: fformr=52
 integer :: accessfil,ask_accurate,band_index,bantot_rbz,bd2tot_index,bdtot1_index
 integer :: bdtot_index,choice,cplex,ddkcase,formeig,gscase,iatom,iband,idir
 integer :: ierr,ifft,ii,ikpt,ikpt1,index,initialized,ios,ipert,ipw,ir,ireadwf0
 integer :: iscf_mod,isppol,istr,itime,ixc,localrdwf,master,me,mpw,mpw1,mxfh,n1
 integer :: n2,n3,n3xccc,nband_k,nkpt_eff,nkpt_rbz,nqpt,nsym1,ntypat,nxfh
 integer :: openexit,option,optorth,optthm,pertcase,prtdos,prteig,prtvol,rdwr,rdwrpaw,spaceComm
 integer :: spaceworld,t_iostat,timrev_pert,tmkmem
 integer :: ncerr,ncid_hdr
 real(dp) :: dosdeltae,ecore,ecut_eff,edocc,eei,eeig0,eew,efrhar,efrkin,efrloc
 real(dp) :: efrnl,efrx1,efrx2,ehart,ehart01,ehart1,eig0,eii,ek,ek0,ek1,ek2,eloc0
 real(dp) :: elpsp1,enl,enl0,enl1,entropy,enxc,exc1,fsum,gsqcut,maxocc,mean,nelectkq
 real(dp) :: residm,rho1_dn,rho1im_dn,rho1re_dn,rhosum,tolmxf,tolwfr,ucvol
 logical :: t_exist
 character(len=fnlen) :: fiden1i,fiwf1i,fiwf1o,fiwfddk,wff2nm
 character(len=fnlen) :: gkkfilnam
 character(len=4) :: mode_par_psp
 character(len=6) :: tag
 character(len=500) :: message
 type(bandstructure_type) :: bs_rbz
 type(dens_sym_operator_type) :: densymop_rf
 type(hdr_type) :: hdr,hdr0
 type(wffile_type) :: wff1,wffddk,wffgs,wffkq,wffnew,wffnow,wfftgs,wfftkq
 type(wvl_wf_type) :: wfs
 integer,allocatable :: indkpt1(:),indkpt1_tmp(:),indsy1(:,:,:),irrzon1(:,:,:)
 integer,allocatable :: istwfk_rbz(:),kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:)
 integer,allocatable :: npwarr(:),npwtot(:),npwtot1(:),symaf1(:),symaf1_tmp(:)
 integer,allocatable :: symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:)
 real(dp) :: corstr(6),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2)
 real(dp),allocatable :: cg(:,:),cg1(:,:),cgq(:,:),doccde_rbz(:),docckqde(:)
 real(dp),allocatable :: eigen0(:),eigen1(:),eigen1_diag(:),eigenq(:),kpq(:,:)
 real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_rbz(:),occkq(:)
 real(dp),allocatable :: ph1d(:,:),phnons1(:,:,:),resid(:),rhog1(:,:)
 real(dp),allocatable :: rhor1(:,:),tnons1(:,:),tnons1_tmp(:,:),vpsp1(:)
 real(dp),allocatable :: work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:)
 real(dp),allocatable :: xfhist(:,:,:,:)
 real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:)
 type(pawrhoij_type),allocatable :: pawrhoij(:)

!BEGIN TF_CHANGES
 integer, allocatable :: pert_tmp(:), pert_calc(:)
 integer :: ipert_cnt, icase, igroup_cnt
 integer :: my_group,ngroup_respfn,old_spaceComm
!END TF_CHANGES

!no_abirules
#if !defined MPI
           integer,parameter :: nkpt_max=50
#endif
#if defined MPI
           integer,parameter :: nkpt_max=-1
#endif

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

!DEBUG
!write(6,*)' loper3: enter'
!write(6,*)' loper3 : xred=',xred
!stop
!ENDDEBUG

! write(6,*)' loper3 nkpt',nkpt
 call timab(141,1,tsec)

!BEGIN TF_CHANGES
!Init mpi_comm
 call xcomm_init(mpi_enreg,spaceComm)

#if defined MPI
 if(dtset%paral_rf == 1) then
  old_spaceComm=mpi_enreg%spaceComm
  call initmpi_respfn(mpi_enreg,spaceComm)
 end if
#endif


 master=0
!Init me
 call xme_init(mpi_enreg,me)
!Init mpi_comm
 call xcomm_init(mpi_enreg,spaceComm)
!END TF_CHANGES

 call status(0,dtfil%filstat,iexit,level,'enter         ')

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

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

 prteig=dtset%prteig

!Option input variables
 iscf_mod=dtset%iscf
 ixc=dtset%ixc
 localrdwf=dtset%localrdwf
 ntypat=psps%ntypat

 ecut_eff= (dtset%ecut) * (dtset%dilatmx)**2

 wff2nm=trim(dtfil%filnam_ds(4))//'_1WF'

 ecore=0.0_dp

 optorth=1;if (psps%usepaw==1) optorth=0

!Initialization
 initialized=0
!Are these needed ?
 ek=0.0_dp ; ehart=0.0_dp ; enxc=0.0_dp ; eei=0.0_dp ; enl=0.0_dp
 eii=0.0_dp

 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)

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

!Generate the 1-dimensional phases
 allocate(ph1d(2,3*(2*mgfft+1)*dtset%natom))
 call getph(atindx,dtset%natom,n1,n2,n3,ph1d,xred)
!TF_DBEUG
!write(6,*)"atindx:",atindx
!Some allocations
 cplex=2-timrev
!HERE : impose cplex=2
!cplex=2
!DEBUG
!write(6,*)' loper3 : for this perturbation, timrev,cplex=',timrev,cplex
!ENDDEBUG
 allocate(vpsp1(cplex*nfft))

 call status(0,dtfil%filstat,iexit,level,'before loop   ')

!BEGIN TF_CHANGES

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

!init groups for parallelization over perturbations if paral_rf=1
#if defined MPI
  if(mpi_enreg%paral_compil_respfn==1) then
   my_group=mpi_enreg%my_respfn_group
   ngroup_respfn=dtset%ngroup_rf
  end if
#endif
!END TF_CHANGES

!Variational part : loop on polarisations
!BEGIN TF_CHANGES
 do icase=1,ipert_cnt
     !calculate only private part of perturbations
     !when current perturbation is not part of my group then skip this loop
#if defined MPI
     if(dtset%paral_rf==1) then
       if(mpi_enreg%respfn_group(modulo(icase,ngroup_respfn)+1) /= my_group) then
         cycle
       end if
     end if
#endif

    if (pert_calc(icase) <= dtset%natom*3) then
     idir = mod(pert_calc(icase),3)
     if (idir==0) idir=3
     ipert=( (pert_calc(icase)-idir) / 3 + 1)
    else
     ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1
     idir = mod(pert_calc(icase),3)
     if (idir==0) idir=3
    end if
!END TF_CHANGES

     pertcase=idir+(ipert-1)*3
     call status(pertcase,dtfil%filstat,iexit,level,'enter loop    ')

     write(message, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,&
&     ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:)
     call wrtout(6,message,'COLL')
     call wrtout(ab_out,message,'COLL')

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

!    Determine the subset of symmetry operations (nsym1 operations)
!    that leaves the perturbation invariant,
!    and initialize corresponding arrays symaf1, symrl1 and tnons1.
     allocate(symaf1_tmp(nsym),symrl1_tmp(3,3,nsym),tnons1_tmp(3,nsym))
     if (dtset%prepanl /= 1 .and. dtset%berryopt /=4 ) then
      call status(pertcase,dtfil%filstat,iexit,level,'call insy3    ')
      call insy3(gprimd,idir,indsym,ipert,dtset%natom,nsym,nsym1,2,&
&      dtset%symafm,symaf1_tmp,symq,symrec,&
&      dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp)
     else
      nsym1 = 1
      symaf1_tmp(1) = 1
      symrl1_tmp(:,:,1) = dtset%symrel(:,:,1)
      tnons1_tmp(:,1) = 0_dp
     end if
     allocate(indsy1(4,nsym1,dtset%natom),symrc1(3,3,nsym1))
     allocate(symaf1(nsym1),symrl1(3,3,nsym1),tnons1(3,nsym1))
     symaf1(1:nsym1)=symaf1_tmp(1:nsym1)
     symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1)
     tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1)
     deallocate(symaf1_tmp,symrl1_tmp,tnons1_tmp)
!     write(6,*)' loper3 3 nkpt',nkpt

!    Set up corresponding symmetry data
     call status(pertcase,dtfil%filstat,iexit,level,'call setsym   ')
     allocate(irrzon1(nfft**(1-1/nsym1),2,nspden/dtset%nsppol))
     allocate(phnons1(2,nfft**(1-1/nsym1),nspden/dtset%nsppol))
     call setsym(densymop_rf,indsy1,irrzon1,1,dtset%natom,&
&     nfft,dtset%ngfft,nspden,dtset%nsppol,nsym1,&
&     phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)

!    Initialize k+q array
     allocate(kpq(3,nkpt))
!    write(6,*)' loper3 4 nkpt',nkpt
     if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
      kpq(:,:)=dtset%kptns(:,:)
     else
      do ikpt=1,nkpt
       kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt)
      end do
     end if

!    Determine the subset of k-points needed in the "reduced Brillouin zone",
!    and initialize other quantities
!    write(6,*)'loper3 1 nkpt',nkpt
     allocate(indkpt1_tmp(nkpt),wtk_folded(nkpt))
     indkpt1_tmp(:)=0
     optthm=0 ; option=1 ; timrev_pert=timrev
!    For the time being, the time reversal symmetry is not used
!    for ddk and elfd perturbations.
     if(ipert==dtset%natom+1 .or. ipert==dtset%natom+2)timrev_pert=0
     call status(pertcase,dtfil%filstat,iexit,level,'call symkpt   ')
     call symkpt(gmet,indkpt1_tmp,dtset%kptns,nkpt,nkpt_rbz,&
      nsym1,option,symrc1,timrev_pert,dtset%wtk,wtk_folded)

     allocate(doccde_rbz(dtset%mband*nkpt_rbz*dtset%nsppol))
     allocate(indkpt1(nkpt_rbz),istwfk_rbz(nkpt_rbz))
     allocate(kpq_rbz(3,nkpt_rbz),kpt_rbz(3,nkpt_rbz))
     allocate(nband_rbz(nkpt_rbz*dtset%nsppol),occ_rbz(dtset%mband*nkpt_rbz*dtset%nsppol))
     allocate(wtk_rbz(nkpt_rbz))
     indkpt1(:)=indkpt1_tmp(1:nkpt_rbz)
     do ikpt=1,nkpt_rbz
      istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt))
      kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt))
      kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt))
      wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt))
     end do
!    Transfer occ to occ_rbz and doccde to doccde_rbz :
!    this is a more delicate issue
!    NOTE : this takes into account that indkpt1 is ordered
     bdtot_index=0
     bdtot1_index=0
     do isppol=1,dtset%nsppol
      ikpt1=1
      do ikpt=1,nkpt
       nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
!      Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1)
       if(ikpt1/=nkpt_rbz+1)then
        if(ikpt==indkpt1(ikpt1))then
         nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k
         occ_rbz(1+bdtot1_index:nband_k+bdtot1_index)=&
&                               occ(1+bdtot_index:nband_k+bdtot_index)
         doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index)=&
&                               doccde(1+bdtot_index:nband_k+bdtot_index)
         ikpt1=ikpt1+1
         bdtot1_index=bdtot1_index+nband_k
        end if
       end if
       bdtot_index=bdtot_index+nband_k
      end do
     end do

     deallocate(indkpt1_tmp,wtk_folded)

     call timab(142,1,tsec)
!    Compute maximum number of planewaves at k
     call status(0,dtfil%filstat,iexit,level,'call getmpw-k ')
     call getmpw(ecut_eff,dtset%exchn2n3,gmet,istwfk_rbz,kpt_rbz,&
&     mpi_enreg,mpw,nkpt_rbz,ucvol)
     call timab(142,2,tsec)

!    Allocate some arrays
     allocate( kg(3,mpw*mkmem) )
     allocate( npwarr(nkpt_rbz),npwtot(nkpt_rbz) )

#if defined MPI
           allocate(mpi_enreg%proc_distrb(nkpt_rbz,dtset%mband,dtset%nsppol))
           mpi_enreg%parareel=0
           call distrb2(dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,mpi_enreg)
#endif
     call status(0,dtfil%filstat,iexit,level,'call kpgio-k  ')

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

     allocate(ylm(mpw*mkmem,mpsang*mpsang*psps%useylm))
     if (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
      allocate(ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm))
     else
      allocate(ylmgr(1,1,1))
     end if

!    Set up the Ylm for each k point
     if (psps%useylm==1) then
      if(mkmem==0) open(dtfil%unylm,file=tmpfil(10),form='unformatted',status='unknown')
      call status(0,dtfil%filstat,iexit,level,'call initylmg ')
      option=0;if (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4) option=1
      call initylmg(gprimd,kg,kpt_rbz,mkmem,mpi_enreg,mpsang,mpw,nband_rbz,nkpt_rbz,&
&                   npwarr,dtset%nsppol,option,rprimd,dtfil%unkg,dtfil%unylm,ylm,ylmgr)
     end if

     bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol))

!    Open and read pseudopotential files

     call status(0,dtfil%filstat,iexit,level,'call pspini-k ')

     write(message, '(a,a,a)' )ch10,&
&     '--------------------------------------------------------------------------------',&
&     ch10
      call wrtout(ab_out,message,'COLL')
!    Initialize band structure datatype
     allocate(eigen0(bantot_rbz))
     eigen0(:)=zero
     call bstruct_init(bantot_rbz,bs_rbz,doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,&
&     nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,occ_rbz,wtk_rbz)
     deallocate(eigen0)

!    Initialize header
     gscase=0 ! A GS WF file is read
     call hdr_init(bs_rbz,codvsn,dtset,hdr0,pawtab,gscase,psps)

!    Update header, with evolving variables, when available
     call hdr_update(bantot_rbz,etotal,fermie,hdr0,dtset%natom,&
&     residm,rprimd,occ_rbz,pawrhoij,psps%usepaw,xred)

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

     call status(0,dtfil%filstat,iexit,level,'call inwffil-k')

!    Initialize wavefunction files and wavefunctions.
     ireadwf0=1 ; formeig=0 ; ask_accurate=1
     allocate( cg(2,mpw*nspinor*dtset%mband*mkmem*dtset%nsppol) )
     allocate( eigen0(dtset%mband*nkpt_rbz*dtset%nsppol) )
     call timab(144,1,tsec)

     call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3,&
&     formeig,gmet,hdr0,ireadwf0,istwfk_rbz,kg,&
&     kpt_rbz,localrdwf,dtset%mband,&
&     mkmem,mpi_enreg,mpw,nband_rbz,dtset%ngfft,nkpt_rbz,npwarr,nspden,&
&     nspinor,dtset%nsppol,nsym,occ_rbz,optorth,psps,prtvol,rprimd,dtset%symafm,&
&     dtset%symrel,dtset%tnons,dtfil%unkg,wffgs,wfftgs,&
&     dtfil%unwffgs,dtfil%unwftgs,dtfil%fnamewffk,wfs,tmpfil(7))
     call timab(144,2,tsec)

!!    Clean hdr
!     call hdr_clean(hdr)

!    Close wffgs%unwff, if it was ever opened (in inwffil)
     if (ireadwf0==1) then
      call WffClose(wffgs,ierr)
     end if


!    Compute maximum number of planewaves at kpq
!    Will be useful for both GS wfs at k+q and RF wavefunctions
     call timab(143,1,tsec)
     call status(0,dtfil%filstat,iexit,level,'call getmpw-kq')
     call getmpw(ecut_eff,dtset%exchn2n3,gmet,istwfk_rbz,kpq_rbz,&
&     mpi_enreg,mpw1,nkpt_rbz,ucvol)
     call timab(143,2,tsec)

!    Allocate some arrays
     allocate( kg1(3,mpw1*mk1mem) )
     allocate( npwar1(nkpt_rbz),npwtot1(nkpt_rbz) )

     call status(0,dtfil%filstat,iexit,level,'call kpgio(2) ')

!    Set up the basis sphere of planewaves at kpq
!    Will be useful for both GS wfs at k+q and RF wavefunctions
     call timab(142,1,tsec)
     call kpgio(ecut_eff,dtset%exchn2n3,gmet,istwfk_rbz,kg1,tmpfil(5),&
&     kpq_rbz,mk1mem,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1,&
&     npwar1,npwtot1,dtset%nsppol,dtfil%unkg1)
     call timab(142,2,tsec)
     allocate(ylm1(mpw1*mk1mem,mpsang*mpsang*psps%useylm))
     if (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
      allocate(ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*psps%useylm))
     else
      allocate(ylmgr1(1,1,1))
     end if

!    Set up the Ylm for each kpq point
     if (psps%useylm==1) then
      if(mkmem==0) open(dtfil%unylm1,file=tmpfil(11),form='unformatted',status='unknown')
      call status(0,dtfil%filstat,iexit,level,'call initylmg ')
      option=0;if (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4) option=1
      call initylmg(gprimd,kg1,kpq_rbz,mk1mem,mpi_enreg,mpsang,mpw1,nband_rbz,nkpt_rbz,&
&                   npwar1,dtset%nsppol,option,rprimd,dtfil%unkg1,dtfil%unylm1,ylm1,ylmgr1)
     end if

     call status(0,dtfil%filstat,iexit,level,'call pspini-kq')

     write(message, '(a,a)' )&
&     '--------------------------------------------------------------------------------',&
&     ch10
     call wrtout(ab_out,message,'COLL')

!    Initialize band structure datatype
     allocate(eigenq(bantot_rbz))
     eigenq(:)=zero
     call bstruct_init(bantot_rbz,bs_rbz,doccde_rbz,eigenq,istwfk_rbz,kpq_rbz,&
&     nband_rbz,nkpt_rbz,npwar1,dtset%nsppol,occ_rbz,wtk_rbz)
     deallocate(eigenq)

!    Initialize header
     call hdr_init(bs_rbz,codvsn,dtset,hdr,pawtab,pertcase,psps)

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

     call status(0,dtfil%filstat,iexit,level,'call inwffilkq')

!    Initialize k+q wavefunction files and wavefunctions.
     ireadwf0=1 ; formeig=0 ; ask_accurate=1
     allocate( cgq(2,mpw1*nspinor*dtset%mband*mkqmem*dtset%nsppol) )
     allocate( eigenq(dtset%mband*nkpt_rbz*dtset%nsppol) )
     call timab(144,1,tsec)
     call inwffil(ask_accurate,cgq,dtset,dtset%ecut,ecut_eff,eigenq,dtset%exchn2n3,&
&     formeig,gmet,hdr,&
&     ireadwf0,istwfk_rbz,kg1,kpq_rbz,localrdwf,dtset%mband,&
&     mkqmem,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,&
&     nspden,nspinor,dtset%nsppol,nsym,occ_rbz,optorth,psps,prtvol,&
&     rprimd,dtset%symafm,dtset%symrel,dtset%tnons,&
&     dtfil%unkg1,wffkq,wfftkq,dtfil%unwffkq,dtfil%unwftkq,dtfil%fnamewffq,wfs,tmpfil(8))
     call timab(144,2,tsec)

!    Close dtfil%unwffkq, if it was ever opened (in inwffil)
     if (ireadwf0==1) then
      call WffClose(wffkq,ierr)
     end if

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

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

     if(prtvol==-level)then
      write(message,'(a,a)') ch10,&
&      ' loper3 : initialisation of q part done. '
      call wrtout(06,message,'COLL')
     end if

!    Initialisation of first-order wavefunctions

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

!    Initialize wavefunction files and wavefunctions.
     formeig=1
     ask_accurate=0
     allocate( cg1(2,mpw1*nspinor*dtset%mband*mk1mem*dtset%nsppol) )
     allocate( eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol) )
     allocate( resid(dtset%mband*nkpt_rbz*dtset%nsppol) )
     call status(pertcase,dtfil%filstat,iexit,level,'call inwffil  ')
     call timab(144,1,tsec)

     call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3,&
&     formeig,gmet,hdr,&
&     dtfil%ireadwf,istwfk_rbz,kg1,kpq_rbz,localrdwf,&
&     dtset%mband,mk1mem,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,&
&     nspden,nspinor,dtset%nsppol,nsym1,occ_rbz,optorth,psps,prtvol,rprimd,&
&     symaf1,symrl1,tnons1,dtfil%unkg1,wff1,wffnow,dtfil%unwff1,dtfil%unwft1,&
&     fiwf1i,wfs,tmpfil(1))
     call timab(144,2,tsec)

!DEBUG
!    write(6,*)' loper3 : cg1(:,1)=',cg1(:,1)
!ENDDEBUG

!    Close wff1 (filename fiwf1i), if it was ever opened (in inwffil)
     if (dtfil%ireadwf==1) then
      call WffClose(wff1,ierr)
     end if

!    Initialize second wavefunction file if needed
     if(mkmem==0 .and. dtset%nstep/=0) then
      write(message, '(a,i4,a,a)' )&
&      ' loper3 about to open unit',dtfil%unwft2,' for file=',trim(tmpfil(2))
      call wrtout(06,  message,'PERS')

#if defined HAVE_NETCDF
     if(dtset%accesswff==2) then
    !  Create empty netCDF file
        ncerr = nf90_create(path=trim(tmpfil(2)), cmode=NF90_CLOBBER, ncid=ncid_hdr)
        call handle_ncerr(ncerr," create netcdf wavefunction file")
        ncerr = nf90_close(ncid_hdr)
        call handle_ncerr(ncerr," close netcdf wavefunction file")
     else if(dtset%accesswff==3) then
        write (std_out,*) "FIXME: ETSF I/O support in loper3"
     end if
#endif

      call WffOpen(dtset%accesswff,spaceComm,tmpfil(2),ierr,wffnew,master,me,dtfil%unwft2)
     end if

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

!    Get first-order local potentials.
!    The correct value of gsqcut should be derived from ecut and not ecutsm,
!    but the present trick works
     gsqcut=gsqcut_eff

!    section for strain perturbation
     if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then

      call status(pertcase,dtfil%filstat,iexit,level,'call vlocalstr')
      call vlocalstr(gmet,gprimd,gsqcut,istr,mgfft,mpi_enreg,&
&      psps%mqgrid_vl,dtset%natom,nattyp,nfft,dtset%ngfft,ntypat,ph1d,psps%qgrid_vl,&
&      ucvol,psps%vlspl,vpsp1)

     else

      call status(pertcase,dtfil%filstat,iexit,level,'call vloca3   ')
      call vloca3(atindx,cplex,gmet,gsqcut,idir,ipert,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
&      nattyp,nfft,dtset%ngfft,ntypat,n1,n2,n3,ph1d,psps%qgrid_vl,&
&      dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
     end if
!    end section for strain perturbation

!    Get first-order core correction density change xccc3d1
!    (but do NOT include it in vpsp1 : this will be done in scfcv3
!     because vpsp1 might become spin-polarized)
     n3xccc=0
     if(psps%n1xccc/=0)n3xccc=cplex*nfft
     allocate(xccc3d1(cplex*nfft))
     if(psps%n1xccc/=0)then
      call status(pertcase,dtfil%filstat,iexit,level,'call mkcor3   ')
      call mkcor3(cplex,idir,ipert,dtset%natom,ntypat,n1,psps%n1xccc,&
&      n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
&      psps%xcccrc,psps%xccc1d,xccc3d1,xred)
!DEBUG
!     write(6,*)' ir    vpsp1     kxc     xccc3d1 '
!     do ir=1,nfft,13
!      write(6, '(i5,3es14.6)' ) ir,vpsp1(ir),kxc(ir,1),xccc3d1(ir)
!     end do
!     stop
!ENDDEBUG
     end if ! psps%n1xccc/=0


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

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

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

      rhor1(:,:)=zero ; rhog1(:,:)=zero

     else

      if(iscf_mod>0)then

!      cplex=2 gets the complex density, =1 only real part
       call status(pertcase,dtfil%filstat,iexit,level,'call mkrho3   ')
       call mkrho3(cg,cg1,cplex,densymop_rf,irrzon1,istwfk_rbz,&
&       kg,kg1,dtset%mband,mgfft,mkmem,mk1mem,mpi_enreg,mpw,mpw1,nband_rbz,&
&       nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,dtset%nsppol,nsym1,&
&       occ_rbz,phnons1,rhog1,rhor1,symaf1,ucvol,dtfil%unkg,&
&       dtfil%unkg1,wffnow,wfftgs,wtk_rbz)

      else

       if (me==0) then
        call status(pertcase,dtfil%filstat,iexit,level,'call ioarr    ')
!       Read rho1(r) from a disk file
        rdwr=1;rdwrpaw=0
        call appdig(pertcase,dtfil%fildens1in,fiden1i)
        call ioarr(accessfil,rhor1,etotal,fformr,fiden1i,hdr,&
&        cplex*nfft,nspden,rdwr,rdwrpaw,dtset%ngfft)
       end if

#if defined MPI
!BEGIN TF_CHANGES
         call MPI_BCAST(rhor1,cplex*nfft*nspden,&
&            MPI_DOUBLE_PRECISION,0,spaceComm,ierr)
!END TF_CHANGES
#endif

!      Compute up+down rho1(G) by fft
       call status(pertcase,dtfil%filstat,iexit,level,'call fourdp   ')
       allocate(work(cplex*nfft))
       work(:)=rhor1(:,1)
       call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfft,dtset%ngfft,0)
       deallocate(work)

      end if

     end if
!    Check whether exiting was required by the user.
!    If found then do not start minimization steps
     openexit=1 ; if(dtset%chkexit==0) openexit=0
     call chkexi(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg,openexit)
!    If immediate exit, and wavefunctions were not read, must zero eigenvalues
     if (iexit/=0)eigen1(:)=0.0_dp
     if (iexit==0) then
      call status(pertcase,dtfil%filstat,iexit,level,'call scfcv3   ')
      call scfcv3(atindx,atindx1,blkflg,cg,cgq,cg1,cplex,cpus,&
&      densymop_rf,doccde_rbz,docckqde,dtfil,dtset,&
&      d2bbb,d2lo,d2nl,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
&      ehart,ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek,ek0,ek1,eloc0,elpsp1,&
&      enl,enl0,enl1,enxc,etotal,exc1,fermie,fform,hdr,idir,indkpt1,&
&      indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
&      kg,kg1,kpt_rbz,kxc,mgfft,mkmem,mkqmem,mk1mem,&
&      mpert,mpi_enreg,mpsang,mpw,mpw1,&
&      nattyp,nband_rbz,nfft,nkpt,nkpt_rbz,nkxc,&
&      npwarr,npwar1,nspden,nspinor,&
&      nsym1,n3xccc,occkq,occ_rbz,pertcase,phnons1,ph1d,prtbbb,psps,&
&      dtset%qptn,resid,residm,rhog,rhog1,&
&      rhor,rhor1,rprimd,symaf1,symrc1,symrl1,timrev,&
&      tnons1,wffddk,wffnew,wffnow,wfftgs,wfftkq,vpsp1,vtrial,&
&      wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,dtset%ngfft)

!    End of the check of hasty exit
     end if

     call timab(146,1,tsec)
     write(message, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,&
&     ' ----iterations are completed or convergence reached----',&
&    ch10
     call wrtout(ab_out,message,'COLL')
     call wrtout(06,  message,'COLL')

!    Update the content of the header (evolving variables)
     call hdr_update(bantot_rbz,etotal,fermie,hdr,dtset%natom,&
&     residm,rprimd,occ_rbz,pawrhoij,psps%usepaw,xred)

     call status(pertcase,dtfil%filstat,iexit,level,'call outwf    ')

     mxfh=0 ; nxfh=0 ; nqpt=1
     allocate(xfhist(3,dtset%natom+4,2,mxfh))
     call outwf(cg1,dtfil,dtset,eigen1,fiwf1o,hdr,kg1,kpt_rbz,&
&     dtset%mband,mk1mem,mpi_enreg,mpw1,mxfh,dtset%natom,nband_rbz,nfft,&
&     dtset%ngfft,nkpt_rbz,npwar1,nqpt,nspinor,dtset%nsppol,dtset%nstep,&
&     nxfh,occ_rbz,resid,response,wffnow,wfs,xfhist)
     deallocate(xfhist)
     if(mkmem==0) then
        call WffDelete(wffnow,ierr)
     end if

!    print _gkk file for this perturbation
     if (dtset%prtgkk == 1) then
      gkkfilnam = trim(dtfil%filnam_ds(4))//"_GKK"
      write (*,*) 'gkkfilnam = ', trim(gkkfilnam)
      call outgkk(bantot_rbz,2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol,gkkfilnam,&
    &      dtset,eigen0,eigen1,hdr0,hdr)
     end if

!    If the perturbation is d/dk, evaluate the f-sum rule.
     if(ipert==dtset%natom+1)then
!     Note : the factor of two is related to the difference
!     between Taylor expansion and perturbation expansion
!     Note : this expression should be modified for ecutsm.
!     Indeed, the present one will NOT tend to 1.0_dp.
      ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect
      fsum=-ek1/ek2
      if(dtset%ecutsm<1.0d-6)then
       write(message, '(a,es20.10,a,a,es20.10)' ) &
&       ' loper3 : ek2=',ek2,ch10,&
&       '          f-sum rule ratio=',fsum
      else
       write(message, '(a,es20.10,a,a,es20.10,a)' ) &
&       ' loper3 : ek2=',ek2,ch10,&
&       '          f-sum rule ratio=',fsum,' (note : ecutsm/=0)'
      end if
      call wrtout(6,message,'COLL')
      call wrtout(ab_out,message,'COLL')


!     Write the diagonal elements of the dH/dk operator
      allocate(eigen1_diag(dtset%mband*nkpt_rbz*dtset%nsppol))
      bdtot_index=0 ; bd2tot_index=0
      do isppol=1,dtset%nsppol
       do ikpt=1,nkpt_rbz
        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
        do iband=1,nband_k
         eigen1_diag(iband+bdtot_index)=&
&         eigen1(2*iband-1 + (iband-1)*2*nband_k + bd2tot_index)
        end do

!       Treat the case of degeneracies : take the mean of degenerate states
        if(nband_k>1)then
         eig0=eigen0(1+bdtot_index)
         ii=1
         do iband=2,nband_k
          if(eigen0(iband+bdtot_index)-eig0<tol8)then
           ii=ii+1
          else
           mean=sum(eigen1_diag(iband-ii+bdtot_index:iband-1+bdtot_index))/ii
           eigen1_diag(iband-ii+bdtot_index:iband-1+bdtot_index)=mean
           ii=1
          end if
          eig0=eigen0(iband+bdtot_index)
          if(iband==nband_k)then
           mean=sum(eigen1_diag(iband-ii+1+bdtot_index:iband+bdtot_index))/ii
           eigen1_diag(iband-ii+1+bdtot_index:iband+bdtot_index)=mean
          end if
         end do
        end if

        bdtot_index=bdtot_index+nband_k
        bd2tot_index=bd2tot_index+2*nband_k**2
       end do
      end do
      option=4
      call prteigrs(eigen1_diag,dtset%enunit,fermie,tmpfil(1),ab_out,iscf_mod,kpt_rbz,dtset%kptopt,&
&      dtset%mband,nband_rbz,nkpt_rbz,dtset%nnsclo,dtset%nsppol,occ_rbz,dtset%occopt,&
&      option,prteig,prtvol,resid,tolwfr,vxcavg,wtk_rbz)
      deallocate(eigen1_diag)

     end if

!    Print the energies
     if (dtset%nline/=0 .or. dtset%nstep/=0)then
      call status(pertcase,dtfil%filstat,iexit,level,'call prtene3  ')
      call prtene3(edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
&      ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,exc1,ab_out,&
&      ipert,dtset%natom)
     end if

     deallocate(cg,cgq,cg1,docckqde,doccde_rbz,eigen0,eigenq,eigen1,kpq)
     deallocate(indkpt1,indsy1,istwfk_rbz)
     deallocate(irrzon1,kg,kg1,kpq_rbz,kpt_rbz)
     deallocate(nband_rbz,npwarr,npwar1,npwtot,npwtot1,occkq,occ_rbz,phnons1,resid,rhog1,rhor1)
     deallocate(symaf1,symrc1,symrl1,tnons1,wtk_rbz,xccc3d1)
     deallocate(ylm,ylm1,ylmgr,ylmgr1)

#if defined MPI
           deallocate(mpi_enreg%proc_distrb)
#endif

!    Clean hdr
     call hdr_clean(hdr)
     call hdr_clean(hdr0)

!    Close the unneeded temporary data files, if any
     if (mkmem==0) then
      close (unit=dtfil%unkg,status='delete')
      close (unit=dtfil%unylm,status='delete')
      call WffDelete(wfftgs,ierr)
     end if
     if (mkqmem==0) then
      call WffDelete(wfftkq,ierr)
     end if
     if (mk1mem==0) then
      close (unit=dtfil%unkg1,status='delete')
      close (unit=dtfil%unylm1,status='delete')
      call WffDelete(wffnew,ierr)
     end if

     call timab(146,2,tsec)
!BEGIN TF_CHANGES
  if(iexit/=0)exit
!End loop over icase
 end do

#if defined MPI
          !test whether _all_ groups are working properly
          !set communicator to old_spaceComm before
           if(mpi_enreg%paral_compil_respfn == 1) then
            mpi_enreg%spaceComm=old_spaceComm
            call leave_test(mpi_enreg)
          ! reset communicator
            mpi_enreg%spaceComm=spaceComm
           end if
#endif

!END TF_CHANGES

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

 call status(0,dtfil%filstat,iexit,level,'end loop      ')

!Deallocate arrays
 deallocate(ph1d,vpsp1)
!BEGIN TF_CHANGES
 deallocate(pert_calc)
!END TF_CHANGES

!BEGIN TF_CHANGES
#if defined MPI
!free groups
 if(mpi_enreg%paral_compil_respfn == 1) then
  do igroup_cnt=1,dtset%ngroup_rf
        if(mpi_enreg%respfn_comm(igroup_cnt) == mpi_enreg%spaceComm) then

                !DEBUG
                !PRINT *,mpi_enreg%me,":Freeing Comm - comm=",mpi_enreg%respfn_comm(igroup_cnt)

                call MPI_COMM_FREE(mpi_enreg%respfn_comm(igroup_cnt),ierr)
                if(ierr /= MPI_SUCCESS) then
                        PRINT *,mpi_enreg%me,': Error on releasing Communicator for group nr ',igroup_cnt
                        call leave_new('COLL')
                end if

        end if
  end do

  deallocate(mpi_enreg%respfn_group)
  deallocate(mpi_enreg%respfn_comm)
  mpi_enreg%spaceComm=old_spaceComm
  mpi_enreg%me_respfn=-1
 end if
#endif
!END TF_CHANGES

 write(message, '(a,a)' ) ch10,' loper3 : exiting '
 call wrtout(06,message,'COLL')

 call status(0,dtfil%filstat,iexit,level,'exit          ')

 call timab(141,2,tsec)
!DEBUG
!write(6,*)' loper3: exit '
!ENDDEBUG
end subroutine loper3
!!***

Generated by  Doxygen 1.6.0   Back to index