TABLE OF CONTENTS


ABINIT/dfpt_dyfro [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_dyfro

FUNCTION

 Compute the different parts of the frozen-wavefunction part of
 the dynamical matrix, except the non-local one, computed previously.
 Also (when installed) symmetrize the different part and their sum.

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
  dyfr_nondiag=1 if dyfrnl and dyfrwf are non diagonal with respect to atoms; 0 otherwise
  gmet(3,3)=reciprocal space metric (bohr^-2)
  gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
  gsqcut=cutoff on G^2 based on ecut
  indsym(4,nsym,natom)=index showing transformation of atom labels
   under symmetry operations (computed in symatm)
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  mqgrid=dimensioned number of q grid points for local psp spline
  natom=number of atoms in unit cell
  nattyp(ntypat)=number of atoms of each type
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT,
     see ~abinit/doc/variables/vargs.htm#ngfft
  nspden=number of spin-density components
  nsym=number of symmetries in space group
  ntypat=number of types of atoms
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  n3xccc=dimension of the xccc3d array (0 or nfft).
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  qgrid(mqgrid)=q point array for local psp spline fits
  qphon(3)=wavevector of the phonon
  rhog(2,nfft)=electron density in G space
  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q
  symrec(3,3,nsym)=symmetries in reciprocal space
  typat(natom)=integer type for each atom in cell
  ucvol=unit cell volume (bohr**3).
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.
  vxc(nfft,nspden)=exchange-correlation potential (hartree) in real
   space--only used when n1xccc/=0
  xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
  xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
   for each type of atom, from psp
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  xred(3,natom)=reduced coordinates for atoms in unit cell

OUTPUT

  dyfrlo(3,3,natom)=frozen wavefunctions part of the dynamical matrix
                    (local only)
  dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
                    frozen wavefunctions part of the dynamical matrix
                    (local + non-local)
                    If NCPP, it depends on one atom
                    If PAW,  it depends on two atoms
  dyfrxc(3,3,natom)=frozen wavefunctions part of the dynamical matrix
                    (non-linear xc core correction)

SIDE EFFECTS

 Input/Output
  dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
                    frozen wavefunctions part of the dynamical matrix
                    (non-local only)
                    If NCPP, it depends on one atom
                    If PAW,  it depends on two atoms

PARENTS

      respfn

CHILDREN

      atm2fft,dfpt_sydy,fourdp,mkcore,mklocl_recipspace,timab,zerosym

SOURCE

3999 subroutine dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrxc,dyfr_cplex,dyfr_nondiag,&
4000 &  gmet,gprimd,gsqcut,indsym,mgfft,mpi_enreg,mqgrid,natom,nattyp,&
4001 &  nfft,ngfft,nspden,nsym,ntypat,n1xccc,n3xccc,paral_kgb,psps,pawtab,ph1d,qgrid,&
4002 &  qphon,rhog,rprimd,symq,symrec,typat,ucvol,usepaw,vlspl,vxc,&
4003 &  xcccrc,xccc1d,xccc3d,xred)
4004 
4005 
4006 !This section has been created automatically by the script Abilint (TD).
4007 !Do not modify the following lines by hand.
4008 #undef ABI_FUNC
4009 #define ABI_FUNC 'dfpt_dyfro'
4010 !End of the abilint section
4011 
4012  implicit none
4013 
4014 !Arguments ------------------------------------
4015 !scalars
4016  integer,intent(in) :: dyfr_cplex,dyfr_nondiag,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden
4017  integer,intent(in) :: nsym,ntypat,paral_kgb,usepaw
4018  real(dp),intent(in) :: gsqcut,ucvol
4019  type(pseudopotential_type),intent(in) :: psps
4020  type(MPI_type),intent(in) :: mpi_enreg
4021 !arrays
4022  integer,intent(in) :: atindx1(natom),indsym(4,nsym,natom),nattyp(ntypat)
4023  integer,intent(in) :: ngfft(18),symq(4,2,nsym),symrec(3,3,nsym),typat(natom)
4024  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
4025  real(dp),intent(in) :: qgrid(mqgrid),qphon(3),rhog(2,nfft),rprimd(3,3)
4026  real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden)
4027  real(dp),intent(in) :: xccc1d(n1xccc,6,ntypat),xcccrc(ntypat),xred(3,natom)
4028  real(dp),intent(inout) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag),xccc3d(n3xccc)
4029  real(dp),intent(out) :: dyfrlo(3,3,natom),dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
4030  real(dp),intent(inout) :: dyfrxc(3,3,natom) !vz_i
4031  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
4032 
4033 !Local variables-------------------------------
4034 !scalars
4035  logical, parameter :: do_final_sym=.true.
4036  integer :: iatom,jatom,n1,n2,n3,optatm,optdyfr,opteltfr,optgr,option
4037  integer :: optn,optn2,optstr,optv
4038  real(dp) :: eei
4039 !arrays
4040  integer  :: qprtrb(3)
4041  real(dp) :: dummy6(6),dum_strn(6),dum_strv(6)
4042  real(dp) :: tsec(2),vprtrb(2)
4043  real(dp) :: dum_atmrho(0),dum_atmvloc(0),dum_gauss(0),dum_grn(0),dum_grv(0),dum_eltfrxc(0)
4044  real(dp),allocatable :: dyfrlo_tmp1(:,:,:),dyfrlo_tmp2(:,:,:,:,:),dyfrsym_tmp(:,:,:,:,:)
4045  real(dp),allocatable :: gr_dum(:,:),v_dum(:),vxctotg(:,:)
4046 
4047 ! *************************************************************************
4048 
4049  if(nspden==4)then
4050    MSG_WARNING('dfpt_dyfro : DFPT with nspden=4 works at the moment just for insulators and norm-conserving psp!')
4051  end if
4052 
4053  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
4054 
4055  if (usepaw==1 .or. psps%nc_xccc_gspace==1) then
4056 
4057 !  PAW or NC with nc_xccc_gspace: compute local psp and core charge contribs together
4058 !  in reciprocal space
4059 !  -----------------------------------------------------------------------
4060    call timab(563,1,tsec)
4061    if (n3xccc>0) then
4062      ABI_ALLOCATE(v_dum,(nfft))
4063      ABI_ALLOCATE(vxctotg,(2,nfft))
4064      v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2))
4065      call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
4066      call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),&
4067 &     comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft)
4068      ABI_DEALLOCATE(v_dum)
4069    else
4070      ABI_ALLOCATE(vxctotg,(0,0))
4071    end if
4072    optatm=0;optdyfr=1;optgr=0;optstr=0;optv=1;optn=n3xccc/nfft;optn2=1;opteltfr=0
4073    call atm2fft(atindx1,dum_atmrho,dum_atmvloc,dyfrxc,dyfrlo,dum_eltfrxc,&
4074 &   dum_gauss,gmet,gprimd,dum_grn,dum_grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,&
4075 &   optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,qprtrb,rhog,&
4076 &   dum_strn,dum_strv,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb,vlspl)
4077    ABI_DEALLOCATE(vxctotg)
4078    if (n3xccc==0) dyfrxc=zero
4079  else
4080 
4081 !  Norm-conserving: compute local psp contribution in reciprocal space
4082 !  and core charge contribution in real space
4083 !  -----------------------------------------------------------------------
4084    option=4
4085    ABI_ALLOCATE(dyfrlo_tmp1,(3,3,natom))
4086    ABI_ALLOCATE(gr_dum,(3,natom))
4087    ABI_ALLOCATE(v_dum,(nfft))
4088    call mklocl_recipspace(dyfrlo_tmp1,eei,gmet,gprimd,&
4089 &   gr_dum,gsqcut,dummy6,mgfft,mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,&
4090 &   ntypat,option,paral_kgb,ph1d,qgrid,qprtrb,rhog,ucvol,vlspl,vprtrb,v_dum)
4091    do iatom=1,natom
4092 !    Reestablish correct order of atoms
4093      dyfrlo(1:3,1:3,atindx1(iatom))=dyfrlo_tmp1(1:3,1:3,iatom)
4094    end do
4095    ABI_DEALLOCATE(dyfrlo_tmp1)
4096    ABI_DEALLOCATE(v_dum)
4097    if(n1xccc/=0)then
4098      call mkcore(dummy6,dyfrxc,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,&
4099 &     n1,n1xccc,n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred)
4100    end if
4101    ABI_DEALLOCATE(gr_dum)
4102  end if
4103 
4104 !Symmetrize dynamical matrix explicitly for given space group:
4105 
4106 !Symmetrize local part of the dynamical matrix dyfrlo:
4107  ABI_ALLOCATE(dyfrsym_tmp,(1,3,3,natom,1))
4108  ABI_ALLOCATE(dyfrlo_tmp2,(1,3,3,natom,1))
4109  dyfrsym_tmp(1,:,:,:,1)=dyfrlo(:,:,:)
4110  call dfpt_sydy(1,dyfrsym_tmp,indsym,natom,0,nsym,qphon,dyfrlo_tmp2,symq,symrec)
4111  dyfrlo(:,:,:)=dyfrlo_tmp2(1,:,:,:,1)
4112  if (do_final_sym) then
4113    dyfrsym_tmp(1,:,:,:,1)=dyfrxc(:,:,:)
4114    call dfpt_sydy(1,dyfrsym_tmp,indsym,natom,0,nsym,qphon,dyfrlo_tmp2,symq,symrec)
4115    dyfrxc(:,:,:)=dyfrlo_tmp2(1,:,:,:,1)
4116  end if
4117  ABI_DEALLOCATE(dyfrsym_tmp)
4118  ABI_DEALLOCATE(dyfrlo_tmp2)
4119 
4120 !Symmetrize nonlocal part of the dynamical matrix dyfrnl:
4121 !atindx1 is used to reestablish the correct order of atoms
4122  if (dyfr_nondiag==0) then
4123    ABI_ALLOCATE(dyfrsym_tmp,(dyfr_cplex,3,3,natom,1))
4124    do iatom=1,natom
4125      dyfrsym_tmp(:,:,:,atindx1(iatom),1)=dyfrnl(:,:,:,iatom,1)
4126    end do
4127  else
4128    ABI_ALLOCATE(dyfrsym_tmp,(dyfr_cplex,3,3,natom,natom))
4129    do jatom=1,natom
4130      do iatom=1,natom
4131        dyfrsym_tmp(:,:,:,atindx1(iatom),atindx1(jatom))=dyfrnl(:,:,:,iatom,jatom)
4132      end do
4133    end do
4134  end if
4135  call dfpt_sydy(dyfr_cplex,dyfrsym_tmp,indsym,natom,dyfr_nondiag,nsym,qphon,dyfrnl,symq,symrec)
4136  ABI_DEALLOCATE(dyfrsym_tmp)
4137 
4138 !Collect local, nl xc core, and non-local part
4139 !of the frozen wf dynamical matrix.
4140  dyfrwf(:,:,:,:,:)=dyfrnl(:,:,:,:,:)
4141  if (dyfr_nondiag==0) then
4142    dyfrwf(1,:,:,:,1)=dyfrwf(1,:,:,:,1)+dyfrlo(:,:,:)+dyfrxc(:,:,:)
4143  else
4144    do iatom=1,natom
4145      dyfrwf(1,:,:,iatom,iatom)=dyfrwf(1,:,:,iatom,iatom)+dyfrlo(:,:,iatom)+dyfrxc(:,:,iatom)
4146    end do
4147  end if
4148 
4149 end subroutine dfpt_dyfro

ABINIT/dfpt_dyout [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_dyout

FUNCTION

 Output of all quantities related to the 2nd-order matrix:
 Ewald part, local and non-local frozen wf part,
 core contributions,
 local and non-local variational part, 2nd-order
 matrix itself, and, for the phonon part,
 eigenfrequencies, in Hartree, meV and cm-1.
 Also output unformatted 2nd-order matrix for later
 use in the Brillouin-zone interpolation

INPUTS

  becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only)
  blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical
  matrix has been calculated ; 0 otherwise )
  carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
  2DTE matrix has been calculated correctly ; 0 otherwise )
  ddboun=unit number for the derivative database output
  ddkfil(3)=components are 1 if corresponding d/dk file exists, otherwise 0
  (in what follows, DYMX means dynamical matrix, and D2MX means 2nd-order matrix)
  dyew(2,3,natom,3,natom)=Ewald part of the DYMX
  dyfrlo(3,3,natom)=frozen wf local part of the DYMX
  dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf nonloc part of the DYMX
  dyfrx1(2,3,natom,3,natom)=frozen wf nonlin. xc core corr.(1)
    part of the DYMX
  dyfrx2(3,3,natom)=frozen wf nonlin. xc core corr.(2) part of the DYMX
  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
  dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix
  d2cart(2,3,mpert,3,mpert)=D2MX in cartesian coordinates
  d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)=
    band by band decomposition of Born effective charges
    (calculated from phonon-type perturbation) in cartesian coordinates
  d2eig0(2,3,mpert,3,mpert)=0-order eigen. station. part of the D2MX
  d2k0(2,3,mpert,3,mpert)=0-order kinet. station. part of the D2MX
  d2lo(2,3,mpert,3,mpert)=nonstation. local part of the D2MX
  d2loc0(2,3,mpert,3,mpert)=0-order loc station. part of the D2MX
  d2matr(2,3,mpert,3,mpert)=D2MX in non-cartesian coordinates
  d2nl(2,3,mpert,3,mpert)=nonstation. nonloc part of the D2MX
  d2nl0(2,3,mpert,3,mpert)=0-order nonloc station. part of the D2MX
  d2nl1(2,3,mpert,3,mpert)=1-order nonloc station. part of the D2MX
  d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs (PAW)
  d2vn(2,3,mpert,3,mpert)=potential*dens station. part of the D2MX and without masses included)
  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
  eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor
  has_full_piezo=the full calculation of the piezoelectric tensor from electric field perturbation
                 is only available if nsym=1 (strain perturbation is not symmetrized)
  has_allddk= True if all ddk file are present on disk, so the effective charge or piezzo
              electric tensor are correctly computed (PAW ONLY)
  iout=unit number for long write-up
  mband=maximum number of bands
  mpert =maximum number of ipert
  natom=number of atoms
  ntypat=number of atom types
  outd2=option for the output of the 2nd-order matrix :
   if outd2=1, non-stationary part
   if outd2=2, stationary part.
  pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only)
  pawpiezo= flag for the computation of frozen part of Piezoelectric tensor (PAW only)
  prtbbb=if 1, print the band-by-band decomposition
  prtvol=print volume
  qphon(3)=phonon wavelength, in reduced coordinates
  qzero=1 if zero phonon wavevector
  rfdir(3)=defines the directions for the perturbations
  rfpert(mpert)=defines the perturbations
  rfphon=if 1, there are phonon perturbations
  rfstrs=if 1,2,3 there are strain perturbations
  typat(natom)=integer label of each type of atom (1,2,...)
  usepaw=1 if PAW, 0 otherwise
  usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
  zion(ntypat)=charge corresponding to the atom type

SIDE EFFECTS

  d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)

NOTES

 This routine is called only by the processor me==0 .
 In consequence, no use of message and wrtout routine.

PARENTS

      respfn

CHILDREN

SOURCE

2148 subroutine dfpt_dyout(becfrnl,berryopt,blkflg,carflg,ddboun,ddkfil,dyew,dyfrlo,dyfrnl,&
2149 & dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,&
2150 & d2eig0,d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,&
2151 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
2152 & has_full_piezo,has_allddk,iout,mband,mpert,natom,ntypat,&
2153 & outd2,pawbec,pawpiezo,piezofrnl,prtbbb,prtvol,qphon,qzero,typat,rfdir,&
2154 & rfpert,rfphon,rfstrs,usepaw,usevdw,zion)
2155 
2156 
2157 !This section has been created automatically by the script Abilint (TD).
2158 !Do not modify the following lines by hand.
2159 #undef ABI_FUNC
2160 #define ABI_FUNC 'dfpt_dyout'
2161 !End of the abilint section
2162 
2163  implicit none
2164 
2165 !Arguments -------------------------------
2166 !scalars
2167  integer,intent(in) :: berryopt,ddboun,dyfr_cplex,dyfr_nondiag,iout,mband,mpert
2168  integer,intent(in) :: natom,ntypat,outd2,pawbec,pawpiezo,prtbbb,prtvol,qzero
2169  integer, intent(in) :: rfphon,rfstrs,usepaw,usevdw
2170 !arrays
2171  integer,intent(in) :: blkflg(3,mpert,3,mpert),carflg(3,mpert,3,mpert)
2172  integer,intent(in) :: ddkfil(3),rfdir(3),rfpert(mpert),typat(natom)
2173  real(dp),intent(in) :: becfrnl(3,natom,3*pawbec)
2174  real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),d2eig0(2,3,mpert,3,mpert)
2175  real(dp),intent(in) :: d2k0(2,3,mpert,3,mpert),d2lo(2,3,mpert,3,mpert)
2176  real(dp),intent(in) :: d2loc0(2,3,mpert,3,mpert),d2matr(2,3,mpert,3,mpert)
2177  real(dp),intent(in) :: d2nl(2,3,mpert,3,mpert),d2nl0(2,3,mpert,3,mpert)
2178  real(dp),intent(in) :: d2nl1(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw)
2179  real(dp),intent(in) :: d2vn(2,3,mpert,3,mpert)
2180  real(dp),intent(in) :: dyew(2,3,natom,3,natom),dyfrlo(3,3,natom)
2181  real(dp),intent(in) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
2182  real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyfrx2(3,3,natom)
2183  real(dp),intent(in) :: dyvdw(2,3,natom,3,natom*usevdw)
2184  real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6)
2185  real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6)
2186  real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6)
2187  real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),piezofrnl(6,3*pawpiezo)
2188  real(dp),intent(in) :: qphon(3),zion(ntypat)
2189  real(dp),intent(inout) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)
2190  logical,intent(in) :: has_allddk,has_full_piezo
2191 
2192 !Local variables -------------------------
2193 !scalars
2194  integer :: iband,idir1,idir2,ii,ipert1,ipert2,jj,nelmts,nline
2195  real(dp) :: qptnrm,zi,zr
2196 !arrays
2197  real(dp) :: delta(3,3)
2198 
2199 ! *********************************************************************
2200 
2201 !DEBUG
2202 !write(std_out,*)' dfpt_dyout : enter '
2203 !write(std_out,*)ddkfil
2204 !ENDDEBUG
2205 
2206 !Long print : includes detail of every part of the 2nd-order energy
2207  if(prtvol>=10)then
2208 
2209 !  In case of phonon
2210    if (rfphon==1)then
2211 
2212 !    write the Ewald part of the dynamical matrix
2213      write(iout,*)' '
2214      write(iout,*)' Ewald part of the dynamical matrix'
2215      write(iout,*)'    j1       j2             matrix element'
2216      write(iout,*)' dir pert dir pert     real part   imaginary part'
2217      do ipert1=1,natom
2218        do idir1=1,3
2219          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2220 &         .or.   outd2==1                           )then
2221            write(iout,*)' '
2222            do ipert2=1,natom
2223              do idir2=1,3
2224                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2225                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2226 &                 dyew(1,idir1,ipert1,idir2,ipert2),&
2227 &                 dyew(2,idir1,ipert1,idir2,ipert2)
2228                end if
2229              end do
2230            end do
2231          end if
2232        end do
2233      end do
2234 
2235 !    Now the local frozen wf part
2236      write(iout,*)' '
2237      write(iout,*)' Frozen wf local part of the dynamical matrix'
2238      write(iout,*)'    j1       j2             matrix element'
2239      write(iout,*)' dir pert dir pert     real part   imaginary part'
2240      do ipert1=1,natom
2241        do idir1=1,3
2242          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2243 &         .or.   outd2==1                           )then
2244            write(iout,*)' '
2245            do ipert2=1,natom
2246              do idir2=1,3
2247                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2248                  if(ipert1==ipert2)then
2249                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2250 &                   dyfrlo(idir1,idir2,ipert2),zero
2251                  else
2252                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2253 &                   zero,zero
2254                  end if
2255                end if
2256              end do
2257            end do
2258          end if
2259        end do
2260      end do
2261 
2262 !    Now the nonlo frozen wf part
2263      write(iout,*)' '
2264      write(iout,*)' Frozen wf non-local part of the dynamical matrix'
2265      write(iout,*)'    j1       j2             matrix element'
2266      write(iout,*)' dir pert dir pert     real part   imaginary part'
2267      do ipert1=1,natom
2268        do idir1=1,3
2269          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2270 &         .or.   outd2==1                           )then
2271            write(iout,*)' '
2272            do ipert2=1,natom
2273              do idir2=1,3
2274                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2275                  if(ipert1==ipert2.or.dyfr_nondiag==1)then
2276                    if (dyfr_cplex==1) then
2277                      write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2278 &                     dyfrnl(1,idir1,idir2,ipert1,1+(ipert2-1)*dyfr_nondiag),zero
2279                    else
2280                      write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2281 &                     dyfrnl(:,idir1,idir2,ipert1,1+(ipert2-1)*dyfr_nondiag)
2282                    end if
2283                  else
2284                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2285 &                   zero,zero
2286                  end if
2287                end if
2288              end do
2289            end do
2290          end if
2291        end do
2292      end do
2293 
2294 !    Now the nonlinear xc core correction(1) part
2295      write(iout,*)' '
2296      write(iout,*)' Frozen wf xc core (1) part',&
2297 &     ' of the dynamical matrix'
2298      write(iout,*)'    j1       j2             matrix element'
2299      write(iout,*)' dir pert dir pert     real part   imaginary part'
2300      do ipert1=1,natom
2301        do idir1=1,3
2302          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2303 &         .or.   outd2==1                           )then
2304            write(iout,*)' '
2305            do ipert2=1,natom
2306              do idir2=1,3
2307                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2308                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2309 &                 dyfrx1(1,idir1,ipert1,idir2,ipert2),&
2310 &                 dyfrx1(2,idir1,ipert1,idir2,ipert2)
2311                end if
2312              end do
2313            end do
2314          end if
2315        end do
2316      end do
2317 
2318 !    Now the nonlinear xc core correction(2) part
2319      write(iout,*)' '
2320      write(iout,*)' Frozen wf xc core (2) part',&
2321 &     ' of the dynamical matrix'
2322      write(iout,*)'    j1       j2             matrix element'
2323      write(iout,*)' dir pert dir pert     real part   imaginary part'
2324      do ipert1=1,natom
2325        do idir1=1,3
2326          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2327 &         .or.   outd2==1                           )then
2328            write(iout,*)' '
2329            do ipert2=1,natom
2330              do idir2=1,3
2331                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2332                  if(ipert1==ipert2)then
2333                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2334 &                   dyfrx2(idir1,idir2,ipert2),zero
2335                  else
2336                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2337 &                   zero,zero
2338                  end if
2339                end if
2340              end do
2341            end do
2342          end if
2343        end do
2344      end do
2345 
2346 !    Now the DFT-D vdw part of the dynamical matrix
2347      if (usevdw==1) then
2348        write(iout,*)' '
2349        write(iout,*)' DFT-D van der Waals part of the dynamical matrix'
2350        write(iout,*)'    j1       j2             matrix element'
2351        write(iout,*)' dir pert dir pert     real part   imaginary part'
2352        do ipert1=1,natom
2353          do idir1=1,3
2354            if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2355 &           .or.   outd2==1                           )then
2356              write(iout,*)' '
2357              do ipert2=1,natom
2358                do idir2=1,3
2359                  if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2360                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2361 &                   dyvdw(1,idir1,ipert1,idir2,ipert2),&
2362 &                   dyvdw(2,idir1,ipert1,idir2,ipert2)
2363                  end if
2364                end do
2365              end do
2366            end if
2367          end do
2368        end do
2369      end if
2370 
2371 !    End of the phonon condition
2372    end if
2373 
2374 !  In case of atom. strain/electric field perturbation (piezoelectric tensor)
2375    if (pawpiezo==1.and.(rfpert(natom+2)==1.or.rfstrs/=0).and.outd2==1)then
2376      write(iout,*)' '
2377      write(iout,*)' Frozen wf part of the piezoelectric tensor'
2378      write(iout,*)'    j1       j2             matrix element'
2379      write(iout,*)' dir pert dir pert     real part   imaginary part'
2380      ipert1=natom+2
2381      do idir1=1,3
2382        write(iout,*)' '
2383        ii=1
2384        do ipert2=natom+3,natom+4
2385          do idir2=1,3
2386            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2387 &           piezofrnl(ii,idir1),zero
2388            ii=ii+1
2389          end do
2390        end do
2391      end do
2392    end if
2393 
2394 !  In case of atom. displ/electric field perturbation (Born Effective Charges)
2395    if (pawbec==1.and.(rfpert(natom+2)==1.or.rfphon==1).and.outd2==1)then
2396      write(iout,*)' '
2397      write(iout,*)' Frozen wf part of the Born Effective Charges'
2398      write(iout,*)'    j1       j2             matrix element'
2399      write(iout,*)' dir pert dir pert     real part   imaginary part'
2400      ipert1 = natom+2
2401      do ipert2=1,natom
2402        do idir1=1,3
2403          write(iout,*)' '
2404          do idir2=1,3
2405            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2406 &           becfrnl(idir2,ipert2,idir1),zero
2407          end do
2408        end do
2409      end do
2410    end if
2411 
2412 !  In case of strain
2413    if (rfstrs/=0)then
2414 
2415 !    Write the Ewald part of the elastic tensor
2416      write(iout,*)' '
2417      write(iout,*)' Ewald part of the elastic tensor in cartesian coordinates'
2418      write(iout,*)'    j1       j2             matrix element'
2419      write(iout,*)' dir pert dir pert     real part   imaginary part'
2420      do ipert1=natom+3,natom+4
2421        do idir1=1,3
2422          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2423 &         .or.   outd2==1                           )then
2424            ii=idir1+3*(ipert1-natom-3)
2425            write(iout,*)' '
2426            do ipert2=natom+3,natom+4
2427              do idir2=1,3
2428                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2429                  jj=idir2+3*(ipert2-natom-3)
2430                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2431 &                 elteew(ii,jj),zero
2432                end if
2433              end do
2434            end do
2435          end if
2436        end do
2437      end do
2438 
2439 !    Write the Ewald part of the internal strain coupling parameters
2440      write(iout,*)' '
2441      write(iout,*)' Ewald part of the internal strain coupling parameters'
2442      write(iout,*)'  (cartesian strain, reduced atomic coordinates)'
2443      write(iout,*)'    j1       j2             matrix element'
2444      write(iout,*)' dir pert dir pert     real part   imaginary part'
2445      do ipert1=1,natom
2446        do idir1=1,3
2447          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2448 &         .or.   outd2==1                           )then
2449            ii=idir1+6+3*(ipert1-1)
2450            write(iout,*)' '
2451            do ipert2=natom+3,natom+4
2452              do idir2=1,3
2453                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2454                  jj=idir2+3*(ipert2-natom-3)
2455                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2456 &                 elteew(ii,jj),zero
2457                end if
2458              end do
2459            end do
2460          end if
2461        end do
2462      end do
2463 
2464 !    Now the local frozen wf part
2465      write(iout,*)' '
2466      write(iout,*)' Frozen wf local part of the elastic tensor in cartesian coordinates'
2467      write(iout,*)'    j1       j2             matrix element'
2468      write(iout,*)' dir pert dir pert     real part   imaginary part'
2469      do ipert1=natom+3,natom+4
2470        do idir1=1,3
2471          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2472 &         .or.   outd2==1                           )then
2473            ii=idir1+3*(ipert1-natom-3)
2474            write(iout,*)' '
2475            do ipert2=natom+3,natom+4
2476              do idir2=1,3
2477                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2478                  jj=idir2+3*(ipert2-natom-3)
2479                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2480 &                 eltfrloc(ii,jj),zero
2481                end if
2482              end do
2483            end do
2484          end if
2485        end do
2486      end do
2487 
2488      write(iout,*)' '
2489      write(iout,*)' Frozen wf local part of the internal strain coupling parameters '
2490      write(iout,*)'  (cartesian strain, reduced atomic coordinates)'
2491      write(iout,*)'    j1       j2             matrix element'
2492      write(iout,*)' dir pert dir pert     real part   imaginary part'
2493      do ipert1=1,natom
2494        do idir1=1,3
2495          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2496 &         .or.   outd2==1                           )then
2497            ii=idir1+6+3*(ipert1-1)
2498            write(iout,*)' '
2499            do ipert2=natom+3,natom+4
2500              do idir2=1,3
2501                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2502                  jj=idir2+3*(ipert2-natom-3)
2503                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2504 &                 eltfrloc(ii,jj),zero
2505                end if
2506              end do
2507            end do
2508          end if
2509        end do
2510      end do
2511 
2512 !    Now the nonlo frozen wf part
2513      write(iout,*)' '
2514      write(iout,*)' Frozen wf nonlocal part of the elastic tensor in cartesian coordinates'
2515      write(iout,*)'    j1       j2             matrix element'
2516      write(iout,*)' dir pert dir pert     real part   imaginary part'
2517      do ipert1=natom+3,natom+4
2518        do idir1=1,3
2519          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2520 &         .or.   outd2==1                           )then
2521            ii=idir1+3*(ipert1-natom-3)
2522            write(iout,*)' '
2523            do ipert2=natom+3,natom+4
2524              do idir2=1,3
2525                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2526                  jj=idir2+3*(ipert2-natom-3)
2527                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2528 &                 eltfrnl(ii,jj),zero
2529                end if
2530              end do
2531            end do
2532          end if
2533        end do
2534      end do
2535 
2536      write(iout,*)' '
2537      write(iout,*)' Frozen wf nonlocal part of the internal strain coupling parameters '
2538      write(iout,*)'  (cartesian strain, reduced atomic coordinates)'
2539      write(iout,*)'    j1       j2             matrix element'
2540      write(iout,*)' dir pert dir pert     real part   imaginary part'
2541      do ipert1=1,natom
2542        do idir1=1,3
2543          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2544 &         .or.   outd2==1                           )then
2545            ii=idir1+6+3*(ipert1-1)
2546            write(iout,*)' '
2547            do ipert2=natom+3,natom+4
2548              do idir2=1,3
2549                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2550                  jj=idir2+3*(ipert2-natom-3)
2551                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2552 &                 eltfrnl(ii,jj),zero
2553                end if
2554              end do
2555            end do
2556          end if
2557        end do
2558      end do
2559 
2560 !    Now the xc part
2561      write(iout,*)' '
2562      write(iout,*)' Frozen wf xc part of the elastic tensor in cartesian coordinates'
2563      write(iout,*)'    j1       j2             matrix element'
2564      write(iout,*)' dir pert dir pert     real part   imaginary part'
2565      do ipert1=natom+3,natom+4
2566        do idir1=1,3
2567          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2568 &         .or.   outd2==1                           )then
2569            ii=idir1+3*(ipert1-natom-3)
2570            write(iout,*)' '
2571            do ipert2=natom+3,natom+4
2572              do idir2=1,3
2573                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2574                  jj=idir2+3*(ipert2-natom-3)
2575                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2576 &                 eltfrxc(ii,jj),zero
2577                end if
2578              end do
2579            end do
2580          end if
2581        end do
2582      end do
2583 
2584      write(iout,*)' '
2585      write(iout,*)' Frozen wf xc part of the internal strain coupling parameters '
2586      write(iout,*)'  (cartesian strain, reduced atomic coordinates)'
2587      write(iout,*)'    j1       j2             matrix element'
2588      write(iout,*)' dir pert dir pert     real part   imaginary part'
2589      do ipert1=1,natom
2590        do idir1=1,3
2591          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2592 &         .or.   outd2==1                           )then
2593            ii=idir1+6+3*(ipert1-1)
2594            write(iout,*)' '
2595            do ipert2=natom+3,natom+4
2596              do idir2=1,3
2597                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2598                  jj=idir2+3*(ipert2-natom-3)
2599                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2600 &                 eltfrxc(ii,jj),zero
2601                end if
2602              end do
2603            end do
2604          end if
2605        end do
2606      end do
2607 
2608 !    Now the kinetic frozen wf part
2609      write(iout,*)' '
2610      write(iout,*)' Frozen wf kinetic part of the elastic tensor in cartesian coordinates'
2611      write(iout,*)'    j1       j2             matrix element'
2612      write(iout,*)' dir pert dir pert     real part   imaginary part'
2613      do ipert1=natom+3,natom+4
2614        do idir1=1,3
2615          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2616 &         .or.   outd2==1                           )then
2617            ii=idir1+3*(ipert1-natom-3)
2618            write(iout,*)' '
2619            do ipert2=natom+3,natom+4
2620              do idir2=1,3
2621                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2622                  jj=idir2+3*(ipert2-natom-3)
2623                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2624 &                 eltfrkin(ii,jj),zero
2625                end if
2626              end do
2627            end do
2628          end if
2629        end do
2630      end do
2631 
2632 !    Now the hartree frozen wf part
2633      write(iout,*)' '
2634      write(iout,*)' Frozen wf hartree part of the elastic tensor in cartesian coordinates'
2635      write(iout,*)'    j1       j2             matrix element'
2636      write(iout,*)' dir pert dir pert     real part   imaginary part'
2637      do ipert1=natom+3,natom+4
2638        do idir1=1,3
2639          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2640 &         .or.   outd2==1                           )then
2641            ii=idir1+3*(ipert1-natom-3)
2642            write(iout,*)' '
2643            do ipert2=natom+3,natom+4
2644              do idir2=1,3
2645                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2646                  jj=idir2+3*(ipert2-natom-3)
2647                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2648 &                 eltfrhar(ii,jj),zero
2649                end if
2650              end do
2651            end do
2652          end if
2653        end do
2654      end do
2655 
2656 !    Now the psp core part
2657      write(iout,*)' '
2658      write(iout,*)' Psp core part of the elastic tensor in cartesian coordinates'
2659      write(iout,*)'    j1       j2             matrix element'
2660      write(iout,*)' dir pert dir pert     real part   imaginary part'
2661      do ipert1=natom+3,natom+4
2662        do idir1=1,3
2663          if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2664 &         .or.   outd2==1                           )then
2665            ii=idir1+3*(ipert1-natom-3)
2666            write(iout,*)' '
2667            do ipert2=natom+3,natom+4
2668              do idir2=1,3
2669                if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2670                  jj=idir2+3*(ipert2-natom-3)
2671                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2672 &                 eltcore(ii,jj),zero
2673                end if
2674              end do
2675            end do
2676          end if
2677        end do
2678      end do
2679 
2680 !    Now the DFT-D vdw part
2681      if (usevdw==1) then
2682        write(iout,*)' '
2683        write(iout,*)' DFT-D van der Waals part of the elastic tensor in cartesian coordinates'
2684        write(iout,*)'    j1       j2             matrix element'
2685        write(iout,*)' dir pert dir pert     real part   imaginary part'
2686        do ipert1=natom+3,natom+4
2687          do idir1=1,3
2688            if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2689 &           .or.   outd2==1                           )then
2690              ii=idir1+3*(ipert1-natom-3)
2691              write(iout,*)' '
2692              do ipert2=natom+3,natom+4
2693                do idir2=1,3
2694                  if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2695                    jj=idir2+3*(ipert2-natom-3)
2696                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2697 &                   eltvdw(ii,jj),zero
2698                  end if
2699                end do
2700              end do
2701            end if
2702          end do
2703        end do
2704 
2705        write(iout,*)' '
2706        write(iout,*)' DFT-D2 van der Waals part of the internal strain coupling parameters'
2707        write(iout,*)'  (cartesian strain, reduced atomic coordinates)'
2708        write(iout,*)'    j1       j2             matrix element'
2709        write(iout,*)' dir pert dir pert     real part   imaginary part'
2710        do ipert1=1,natom
2711          do idir1=1,3
2712            if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)&
2713 &           .or.   outd2==1                           )then
2714              ii=idir1+6+3*(ipert1-1)
2715              write(iout,*)' '
2716              do ipert2=natom+3,natom+4
2717                do idir2=1,3
2718                  if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2719                    jj=idir2+3*(ipert2-natom-3)
2720                    write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2721 &                   eltvdw(ii,jj),zero
2722                  end if
2723                end do
2724              end do
2725            end if
2726          end do
2727        end do
2728      end if ! usevdw
2729 
2730 !    End of the strain condition
2731    end if
2732 
2733 !  Now the local nonstationary nonfrozenwf part
2734    if (outd2==1)then
2735      write(iout,*)' '
2736      write(iout,*)' Non-stationary local part of the 2-order matrix'
2737      write(iout,*)'    j1       j2             matrix element'
2738      write(iout,*)' dir pert dir pert     real part   imaginary part'
2739      do ipert1=1,mpert
2740        do idir1=1,3
2741          if ((ipert1<=natom .or.&
2742 &         (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.&
2743 &         ((ipert1==natom+3.or.ipert1==natom+4).and.&
2744 &         (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then
2745            write(iout,*)' '
2746            do ipert2=1,mpert
2747              do idir2=1,3
2748                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2749                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2750 &                 d2lo(1,idir1,ipert1,idir2,ipert2),&
2751 &                 d2lo(2,idir1,ipert1,idir2,ipert2)
2752                end if
2753              end do
2754            end do
2755          end if
2756        end do
2757      end do
2758    end if
2759 
2760 !  Now the nonlocal nonstationary nonfrozenwf part
2761    if (outd2==1)then
2762      write(iout,*)' '
2763      write(iout,*)' Non-stationary non-local part of the 2nd-order matrix'
2764      write(iout,*)'    j1       j2             matrix element'
2765      write(iout,*)' dir pert dir pert     real part   imaginary part'
2766      do ipert1=1,mpert
2767        do idir1=1,3
2768          if ((ipert1<=natom .or.&
2769 &         (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.&
2770 &         ((ipert1==natom+3.or.ipert1==natom+4).and.&
2771 &         (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then
2772            write(iout,*)' '
2773            do ipert2=1,mpert
2774              do idir2=1,3
2775                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2776                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2777 &                 d2nl(1,idir1,ipert1,idir2,ipert2),&
2778 &                 d2nl(2,idir1,ipert1,idir2,ipert2)
2779                end if
2780              end do
2781            end do
2782          end if
2783        end do
2784      end do
2785    end if
2786 
2787 !  Now the overlap change nonstationnary nonfrozenwf part (PAW only)
2788    if (outd2==1.and.usepaw==1)then
2789      write(iout,*)' '
2790      write(iout,*)' PAW: Non-stationary WF-overlap part of the 2nd-order matrix'
2791      write(iout,*)'    j1       j2             matrix element'
2792      write(iout,*)' dir pert dir pert     real part   imaginary part'
2793      do ipert1=1,mpert
2794        do idir1=1,3
2795          if ((ipert1<=natom .or.&
2796 &         (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.&
2797 &         ((ipert1==natom+3.or.ipert1==natom+4).and.&
2798 &         (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then
2799            write(iout,*)' '
2800            do ipert2=1,mpert
2801              do idir2=1,3
2802                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2803                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2804 &                 d2ovl(1,idir1,ipert1,idir2,ipert2),&
2805 &                 d2ovl(2,idir1,ipert1,idir2,ipert2)
2806                end if
2807              end do
2808            end do
2809          end if
2810        end do
2811      end do
2812    end if
2813 
2814 !  Now the 0-order local stationary nonfrozenwf part
2815    if (outd2==2)then
2816      write(iout,*)' '
2817      write(iout,*)' Stationary 0-order local part of the 2nd-order matrix'
2818      write(iout,*)'    j1       j2             matrix element'
2819      write(iout,*)' dir pert dir pert     real part   imaginary part'
2820      do ipert1=1,mpert
2821        do idir1=1,3
2822          if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2823            write(iout,*)' '
2824            do ipert2=1,mpert
2825              do idir2=1,3
2826                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2827                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2828 &                 d2loc0(1,idir1,ipert1,idir2,ipert2),&
2829 &                 d2loc0(2,idir1,ipert1,idir2,ipert2)
2830                end if
2831              end do
2832            end do
2833          end if
2834        end do
2835      end do
2836    end if
2837 
2838 !  Now the stationary 0-order kinetic nonfrozenwf part
2839    if (outd2==2)then
2840      write(iout,*)' '
2841      write(iout,*)' Stationary 0-order kinetic part of the 2nd-order matrix'
2842      write(iout,*)'    j1       j2             matrix element'
2843      write(iout,*)' dir pert dir pert     real part   imaginary part'
2844      do ipert1=1,mpert
2845        do idir1=1,3
2846          if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2847            write(iout,*)' '
2848            do ipert2=1,mpert
2849              do idir2=1,3
2850                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2851                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2852 &                 d2k0(1,idir1,ipert1,idir2,ipert2),&
2853 &                 d2k0(2,idir1,ipert1,idir2,ipert2)
2854                end if
2855              end do
2856            end do
2857          end if
2858        end do
2859      end do
2860    end if
2861 
2862 !  Now the stationary 0-order eigenvalue nonfrozenwf part
2863    if (outd2==2)then
2864      write(iout,*)' '
2865      write(iout,*)' Stationary 0-order eigenvalue part of the'&
2866 &     ,' 2nd-order matrix'
2867      write(iout,*)'    j1       j2             matrix element'
2868      write(iout,*)' dir pert dir pert     real part   imaginary part'
2869      do ipert1=1,mpert
2870        do idir1=1,3
2871          if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2872            write(iout,*)' '
2873            do ipert2=1,mpert
2874              do idir2=1,3
2875                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2876                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2877 &                 d2eig0(1,idir1,ipert1,idir2,ipert2),&
2878 &                 d2eig0(2,idir1,ipert1,idir2,ipert2)
2879                end if
2880              end do
2881            end do
2882          end if
2883        end do
2884      end do
2885    end if
2886 
2887 !  Now the stationary potential-density nonfrozenwf part
2888    if (outd2==2)then
2889      write(iout,*)' '
2890      write(iout,*)' Station. potential-density part of the ',&
2891 &     ' 2nd-order matrix'
2892      write(iout,*)'  (Note : include some xc core-correction) '
2893      write(iout,*)'    j1       j2             matrix element'
2894      write(iout,*)' dir pert dir pert     real part   imaginary part'
2895      do ipert1=1,mpert
2896        do idir1=1,3
2897          if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2898            write(iout,*)' '
2899            do ipert2=1,mpert
2900              do idir2=1,3
2901                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2902                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2903 &                 d2vn(1,idir1,ipert1,idir2,ipert2),&
2904 &                 d2vn(2,idir1,ipert1,idir2,ipert2)
2905                end if
2906              end do
2907            end do
2908          end if
2909        end do
2910      end do
2911    end if
2912 
2913 !  Now the stationary 0-order nonloc nonfrozenwf part
2914    if (outd2==2)then
2915      write(iout,*)' '
2916      write(iout,*)' Stationary 0-order nonlocal part of the 2-order'&
2917 &     ,' matrix'
2918      write(iout,*)'    j1       j2             matrix element'
2919      write(iout,*)' dir pert dir pert     real part   imaginary part'
2920      do ipert1=1,mpert
2921        do idir1=1,3
2922          if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2923            write(iout,*)' '
2924            do ipert2=1,mpert
2925              do idir2=1,3
2926                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2927                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2928 &                 d2nl0(1,idir1,ipert1,idir2,ipert2),&
2929 &                 d2nl0(2,idir1,ipert1,idir2,ipert2)
2930                end if
2931              end do
2932            end do
2933          end if
2934        end do
2935      end do
2936    end if
2937 
2938 !  Now the stationary 1-order nonloc nonfrozenwf part
2939    if (outd2==2)then
2940      write(iout,*)' '
2941      write(iout,*)' Stationary 1-order nonlocal part of the'&
2942 &     ,' 2nd-order matrix'
2943      write(iout,*)' (or the ddk wf part of it, in case of',&
2944 &     ' an electric field perturbation )'
2945      write(iout,*)'    j1       j2             matrix element'
2946      write(iout,*)' dir pert dir pert     real part   imaginary part'
2947      do ipert1=1,mpert
2948        do idir1=1,3
2949          if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
2950            write(iout,*)' '
2951            do ipert2=1,mpert
2952              do idir2=1,3
2953                if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then
2954                  write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
2955 &                 d2nl1(1,idir1,ipert1,idir2,ipert2),&
2956 &                 d2nl1(2,idir1,ipert1,idir2,ipert2)
2957                end if
2958              end do
2959            end do
2960          end if
2961        end do
2962      end do
2963    end if
2964 
2965 !  End of the long print out condition
2966  end if
2967 
2968 
2969 !Derivative database initialisation
2970 
2971 !Calculation of the number of elements
2972  nelmts=0
2973  do ipert1=1,mpert
2974    do idir1=1,3
2975      do ipert2=1,mpert
2976        do idir2=1,3
2977          nelmts=nelmts+blkflg(idir1,ipert1,idir2,ipert2)
2978        end do
2979      end do
2980    end do
2981  end do
2982 
2983  if(outd2==2)then
2984    write(ddboun, '(/,a,i8)' ) ' 2nd derivatives (stationary) - # elements :',nelmts
2985  else if(outd2==1)then
2986    write(ddboun, '(/,a,i8)' ) ' 2nd derivatives (non-stat.)  - # elements :',nelmts
2987  end if
2988 
2989 !Phonon wavevector
2990  qptnrm=1.0_dp
2991 
2992 !Note : if qptnrm should assume another value, it should
2993 !be checked if the f6.1 format is OK.
2994  write(ddboun, '(a,3es16.8,f6.1)' ) ' qpt',(qphon(ii),ii=1,3),qptnrm
2995 
2996 !Now the whole 2nd-order matrix, but not in cartesian coordinates,
2997 !and masses not included
2998  write(iout,*)' '
2999  write(iout,*)' 2nd-order matrix (non-cartesian coordinates,',' masses not included,'
3000  write(iout,*)'  asr not included )'
3001  if(rfstrs/=0) then
3002    write(iout,*)' cartesian coordinates for strain terms (1/ucvol factor '
3003    write(iout,*)'  for elastic tensor components not included) '
3004  end if
3005  write(iout,*)'    j1       j2             matrix element'
3006  write(iout,*)' dir pert dir pert     real part     imaginary part'
3007  nline=1
3008  do ipert1=1,mpert
3009    do idir1=1,3
3010      if(nline/=0)write(iout,*)' '
3011      nline=0
3012      do ipert2=1,mpert
3013        do idir2=1,3
3014          if(blkflg(idir1,ipert1,idir2,ipert2)==1)then
3015            nline=nline+1
3016            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3017 &           d2matr(1,idir1,ipert1,idir2,ipert2),&
3018 &           d2matr(2,idir1,ipert1,idir2,ipert2)
3019            write(ddboun, '(4i4,2d22.14)' )idir1,ipert1,idir2,ipert2,&
3020 &           d2matr(1,idir1,ipert1,idir2,ipert2),&
3021 &           d2matr(2,idir1,ipert1,idir2,ipert2)
3022          end if
3023        end do
3024      end do
3025    end do
3026  end do
3027 
3028 !Now the dynamical matrix
3029  if(rfphon==1)then
3030    write(iout,*)' '
3031    write(iout,*)' Dynamical matrix, in cartesian coordinates,'
3032    write(iout,*)'  if specified in the inputs, asr has been imposed'
3033    write(iout,*)'    j1       j2             matrix element'
3034    write(iout,*)' dir pert dir pert     real part    imaginary part'
3035    nline=1
3036    do ipert1=1,natom
3037      do idir1=1,3
3038        if(nline/=0)write(iout,*)' '
3039        nline=0
3040        do ipert2=1,natom
3041          do idir2=1,3
3042            if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3043              nline=nline+1
3044              write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3045 &             d2cart(1,idir1,ipert1,idir2,ipert2),&
3046 &             d2cart(2,idir1,ipert1,idir2,ipert2)
3047            end if
3048          end do
3049        end do
3050      end do
3051    end do
3052  end if
3053 
3054 !Now the dielectric tensor ! normal case
3055  if(rfpert(natom+2)==1)then
3056 
3057    write(iout,*)' '
3058    write(iout,*)' Dielectric tensor, in cartesian coordinates,'
3059    write(iout,*)'    j1       j2             matrix element'
3060    write(iout,*)' dir pert dir pert     real part    imaginary part'
3061    ipert1=natom+2
3062    ipert2=natom+2
3063    nline=1
3064    do idir1=1,3
3065      if(nline/=0)write(iout,*)' '
3066      nline=0
3067      do idir2=1,3
3068        if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3069          nline=nline+1
3070          write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3071 &         d2cart(1,idir1,ipert1,idir2,ipert2),&
3072 &         d2cart(2,idir1,ipert1,idir2,ipert2)
3073        end if
3074      end do
3075    end do
3076 
3077    if (prtbbb == 1) then
3078 
3079      delta(:,:) = zero
3080      delta(1,1) = one ; delta(2,2) = one ; delta(3,3) = one
3081 
3082      write(iout,*)
3083      write(iout,*)'Band by band decomposition of the dielectric tensor'
3084      write(iout,*)' '
3085 
3086      write(iout,*)' Vacuum polarization'
3087      write(iout,*)'    j1       j2             matrix element'
3088      write(iout,*)' dir pert dir pert     real part    imaginary part'
3089      nline=1
3090      do idir1=1,3
3091        if(nline/=0)write(iout,*)' '
3092        nline=0
3093        do idir2=1,3
3094          nline=nline+1
3095          write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3096 &         delta(idir2,idir1),zero
3097        end do
3098      end do
3099 
3100      do iband = 1,mband
3101        write(iout,*)' '
3102        write(iout,*)' Dielectric tensor, in cartesian coordinates, for band',iband
3103        write(iout,*)'    j1       j2             matrix element'
3104        write(iout,*)' dir pert dir pert     real part    imaginary part'
3105        ipert1 = natom + 2
3106        ipert2 = natom + 2
3107        nline=1
3108        do idir1=1,3
3109          if(nline/=0)write(iout,*)' '
3110          nline=0
3111          do idir2=1,3
3112            if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3113 !            substract vacuum polarization
3114              if (idir1 == idir2) then
3115                d2cart_bbb(1,idir1,idir2,ipert2,iband,iband) = &
3116 &               d2cart_bbb(1,idir1,idir2,ipert2,iband,iband) - 1
3117              end if
3118              nline=nline+1
3119              write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3120 &             d2cart_bbb(1,idir1,idir2,ipert2,iband,iband),&
3121 &             d2cart_bbb(2,idir1,idir2,ipert2,iband,iband)
3122            end if
3123          end do
3124        end do
3125      end do !iband
3126 
3127    end if !prtbbb
3128 
3129  end if ! end natom+2 dielectric output
3130 
3131 !Now the effective charges
3132 !In case of the stationary calculation
3133  if(outd2==2 .and. rfpert(natom+2)==1 .and.rfphon==1)then
3134    write(iout,*)' '
3135    write(iout,*)' Effective charges, in cartesian coordinates,'
3136    write(iout,*)'  if specified in the inputs, asr has been imposed'
3137    write(iout,*)'    j1       j2             matrix element'
3138    write(iout,*)' dir pert dir pert     real part    imaginary part'
3139    ipert1=natom+2
3140    nline=1
3141    do idir1=1,3
3142      if(nline/=0)write(iout,*)' '
3143      nline=0
3144      do ipert2=1,natom
3145        do idir2=1,3
3146          if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3147            nline=nline+1
3148            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3149 &           d2cart(1,idir1,ipert1,idir2,ipert2),&
3150 &           d2cart(2,idir1,ipert1,idir2,ipert2)
3151          end if
3152        end do
3153      end do
3154    end do
3155  end if
3156 
3157 !Now in case of the non-stationary calculation
3158  if(outd2==1 .and. rfpert(natom+2)==1)then
3159    write(iout,*)' '
3160    if(usepaw==1.and..not.(has_allddk))then
3161      write(iout,*)' Warning: Born effectives charges are not correctly computed'
3162      write(iout,*)' you need all ddk perturbations!'
3163    end if
3164    write(iout,*)' Effective charges, in cartesian coordinates,'
3165    write(iout,*)' (from electric field response) '
3166    write(iout,*)'  if specified in the inputs, asr has been imposed'
3167    write(iout,*)'    j1       j2             matrix element'
3168    write(iout,*)' dir pert dir pert     real part    imaginary part'
3169    ipert2=natom+2
3170    nline=1
3171    do idir2=1,3
3172      if(nline/=0)write(iout,*)' '
3173      nline=0
3174      do ipert1=1,natom
3175        do idir1=1,3
3176          if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3177            nline=nline+1
3178            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3179 &           d2cart(1,idir1,ipert1,idir2,ipert2),&
3180 &           d2cart(2,idir1,ipert1,idir2,ipert2)
3181          end if
3182        end do
3183      end do
3184    end do
3185  end if
3186 
3187  if(outd2==1 .and. rfphon==1 .and. qzero==1&
3188 & .and. ( (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0) .or.   &
3189 & berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 ) )then  !!HONG  need to test for fixed E and D
3190    write(iout,*)' '
3191    if(usepaw==1.and..not.(has_allddk))then
3192      write(iout,*)' Warning: Born effectives charges are not correctly computed'
3193      write(iout,*)' you need all ddk perturbations!'
3194    end if
3195    write(iout,*)' Effective charges, in cartesian coordinates,'
3196    write(iout,*)' (from phonon response) '
3197    write(iout,*)'  if specified in the inputs, asr has been imposed'
3198    write(iout,*)'    j1       j2             matrix element'
3199    write(iout,*)' dir pert dir pert     real part    imaginary part'
3200    nline=1
3201    do ipert2=1,natom
3202      do idir2=1,3
3203        if(nline/=0)write(iout,*)' '
3204        nline=0
3205        ipert1=natom+2
3206        do idir1=1,3
3207          if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3208            nline=nline+1
3209            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3210 &           d2cart(1,idir1,ipert1,idir2,ipert2),&
3211 &           d2cart(2,idir1,ipert1,idir2,ipert2)
3212          end if
3213        end do
3214      end do
3215    end do
3216    write(iout,*)' '
3217    write(iout,*)' '
3218    write(iout,*)' '
3219 
3220    if (prtbbb == 1) then
3221 
3222      write(iout,*)'Band by band decomposition of the Born effective charges'
3223      write(iout,*)' '
3224      write(iout,*)'Ionic charges in cartesian coordinates'
3225      write(iout,*)'    j1       j2             matrix element'
3226      write(iout,*)' dir pert dir pert     real part    imaginary part'
3227      zr = zero
3228      zi = zero
3229      do ipert2=1,natom
3230        do idir2=1,3
3231          if(nline/=0)write(iout,*)' '
3232          nline=0
3233          ipert1=natom+2
3234          do idir1=1,3
3235            zr = zero
3236            if (idir1 == idir2) then
3237              zr = zion(typat(ipert2))
3238            end if
3239            nline=nline+1
3240            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3241 &           zr,zi
3242          end do
3243        end do
3244      end do
3245 
3246      do iband = 1,mband
3247        write(iout,*)' '
3248        write(iout,*)' Effective charges, in cartesian coordinates, for band',iband
3249        write(iout,*)' (from phonon response) '
3250        write(iout,*)'  if specified in the inputs, asr has been imposed'
3251        write(iout,*)'    j1       j2             matrix element'
3252        write(iout,*)' dir pert dir pert     real part    imaginary part'
3253        nline=1
3254        do ipert2=1,natom
3255          do idir2=1,3
3256            if(nline/=0)write(iout,*)' '
3257            nline=0
3258            ipert1=natom+2
3259            do idir1=1,3
3260              if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3261                nline=nline+1
3262                write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3263 &               d2cart_bbb(1,idir1,idir2,ipert2,iband,iband),&
3264 &               d2cart_bbb(2,idir1,idir2,ipert2,iband,iband)
3265              end if
3266            end do
3267          end do
3268        end do
3269      end do !iband
3270    end if !prtbbb
3271  end if ! end of print effective charges
3272 
3273 !Now the elastic tensor
3274  if(rfstrs/=0) then
3275    write(iout,*)' '
3276    write(iout,*)' Rigid-atom elastic tensor , in cartesian coordinates,'
3277    write(iout,*)'    j1       j2             matrix element'
3278    write(iout,*)' dir pert dir pert     real part    imaginary part'
3279    nline=1
3280    do ipert1=natom+3,natom+4
3281      do idir1=1,3
3282        if(nline/=0)write(iout,*)' '
3283        nline=0
3284        do ipert2=natom+3,natom+4
3285          do idir2=1,3
3286            if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3287              nline=nline+1
3288              write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3289 &             d2cart(1,idir1,ipert1,idir2,ipert2),&
3290 &             d2cart(2,idir1,ipert1,idir2,ipert2)
3291            end if
3292          end do
3293        end do
3294      end do
3295    end do
3296  end if
3297 
3298 !Now the internal strain coupling parameters
3299  if(rfstrs/=0) then
3300    write(iout,*)' '
3301    write(iout,*)' Internal strain coupling parameters, in cartesian coordinates,'
3302    write(iout,*)'  zero average net force deriv. has been imposed  '
3303    write(iout,*)'    j1       j2             matrix element'
3304    write(iout,*)' dir pert dir pert     real part    imaginary part'
3305    nline=1
3306    do ipert1=1,natom
3307      do idir1=1,3
3308        if(nline/=0)write(iout,*)' '
3309        nline=0
3310        do ipert2=natom+3,natom+4
3311          do idir2=1,3
3312            if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3313              nline=nline+1
3314              write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3315 &             d2cart(1,idir1,ipert1,idir2,ipert2),&
3316 &             d2cart(2,idir1,ipert1,idir2,ipert2)
3317            end if
3318          end do
3319        end do
3320      end do
3321    end do
3322  end if
3323 
3324 !Now the piezoelectric tensor
3325  if(rfstrs/=0 .and. (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0))then
3326    write(iout,*)' '
3327    if(usepaw==1.and..not.(has_allddk))then
3328      write(iout,*)' Warning: Rigid-atom proper piezoelectric tensor is not correctly computed'
3329      write(iout,*)' you need all ddk perturbations!'
3330    end if
3331    write(iout,*)' Rigid-atom proper piezoelectric tensor, in cartesian coordinates,'
3332    write(iout,*)' (from strain response)'
3333    write(iout,*)'    j1       j2             matrix element'
3334    write(iout,*)' dir pert dir pert     real part    imaginary part'
3335    nline=1
3336    ipert1=natom+2
3337    do idir1=1,3
3338      if(nline/=0)write(iout,*)' '
3339      nline=0
3340      do ipert2=natom+3,natom+4
3341        do idir2=1,3
3342          if(carflg(idir1,ipert1,idir2,ipert2)==1)then
3343            nline=nline+1
3344            write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,&
3345 &           d2cart(1,idir1,ipert1,idir2,ipert2),&
3346 &           d2cart(2,idir1,ipert1,idir2,ipert2)
3347          end if
3348        end do
3349      end do
3350    end do
3351  end if
3352 
3353 !Now the piezoelectric tensor
3354  if(outd2==1 .and. (pawpiezo==1.and.rfpert(natom+2)==1)&
3355 & .and. (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0)) then
3356    write(iout,*)' '
3357    if(usepaw==1.and..not.(has_allddk))then
3358      write(iout,*)' Warning: Rigid-atom proper piezoelectric tensor is not correctly computed'
3359      write(iout,*)' you need all ddk perturbations!'
3360    end if
3361    if(usepaw==1.and..not.has_full_piezo)then
3362      write(iout,*)' Warning: The rigid-atom proper piezoelectric tensor'
3363      write(iout,*)' from  electric field response requires nsym=1'
3364    end if
3365    if (has_full_piezo) then
3366      write(iout,*)' Rigid-atom proper piezoelectric tensor, in cartesian coordinates,'
3367      write(iout,*)' (from electric field response)'
3368      write(iout,*)'    j1       j2             matrix element'
3369      write(iout,*)' dir pert dir pert     real part    imaginary part'
3370      nline=1
3371      ipert1=natom+2
3372      do idir1=1,3
3373        if(nline/=0)write(iout,*)' '
3374        nline=0
3375        do ipert2=natom+3,natom+4
3376          do idir2=1,3
3377            if(carflg(idir2,ipert2,idir1,ipert1)==1)then
3378              nline=nline+1
3379              write(iout,'(2(i4,i5),2(1x,f20.10))')idir2,ipert2,idir1,ipert1,&
3380 &             d2cart(1,idir2,ipert2,idir1,ipert1),&
3381 &             d2cart(2,idir2,ipert2,idir1,ipert1)
3382            end if
3383          end do
3384        end do
3385      end do
3386    end if
3387  end if
3388 
3389 end subroutine dfpt_dyout

ABINIT/dfpt_dyxc1 [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_dyxc1

FUNCTION

 Compute 2nd-order non-linear xc core-correction (part1) to the dynamical matrix.
 In case of derivative with respect to k or electric field perturbation,
 the 1st-order local potential vanishes.

INPUTS

  atindx(natom)=index table for atoms
  gmet(3,3)=metrix tensor in G space in Bohr**-2.
  gsqcut=cutoff value on G**2 for sphere inside fft box.
  ixc= choice of exchange-correlation scheme
  kxc(nfft,nkxc)=first-order derivative of the xc potential
    if (nkxc=1) LDA kxc(:,1)= d2Exc/drho2
    if (nkxc=2) LDA kxc(:,1)= d2Exc/drho_up drho_up
                    kxc(:,2)= d2Exc/drho_up drho_dn
                    kxc(:,3)= d2Exc/drho_dn drho_dn
    if (nkxc=7) GGA kxc(:,1)= d2Exc/drho2
                    kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
                    kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
                    kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
                    kxc(:,5)= gradx(rho)
                    kxc(:,6)= grady(rho)
                    kxc(:,7)= gradz(rho)
    if (nkxc=19) spin-polarized GGA case (same as nkxc=7 with up and down components)
  mgfft=maximum size of 1D FFTs
  mpert=maximum number of ipert
  mpi_enreg=information about MPI parallelization
  mqgrid=number of grid pts in q array for f(q) spline.
  natom=number of atoms in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(3)=fft grid dimensions.
  nkxc=second dimension of the kxc array
   (=1 for non-spin-polarized case, =3 for spin-polarized case)
  nspden=number of spin-density components
  ntypat=number of types of atoms in cell.
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
  qgrid(mqgrid)=q grid for spline from 0 to qmax.
  qphon(3)=wavevector of the phonon
  rfdir(3)=array that define the directions of perturbations
  rfpert(mpert)=array defining the type of perturbations that have to be computed
  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  timrev=1 if time-reversal preserves the q wavevector; 0 otherwise.
  typat(natom)=integer type for each atom in cell
  ucvol=unit cell volume (bohr**3).
  usepaw= 0 for non paw calculation; =1 for paw calculation
  xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
  xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
   for each type of atom, from psp
  xred(3,natom)=fractional coordinates for atoms in unit cell

OUTPUT

  blkflgfrx1(3,natom,3,natom)=flag to indicate whether an element has been computed or not
  dyfrx1(2,3,natom,3,natom)=2nd-order non-linear xc
    core-correction (part1) part of the dynamical matrix

PARENTS

      respfn

CHILDREN

      dfpt_atm2fft,dfpt_mkcore,dfpt_mkvxc,dfpt_mkvxc_noncoll,dotprod_vn,timab
      xmpi_sum

SOURCE

4222 subroutine dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,ixc,kxc,mgfft,mpert,mpi_enreg,mqgrid,&
4223 &          natom,nfft,ngfft,nkxc,nspden,ntypat,n1xccc,paral_kgb,psps,pawtab,&
4224 &          ph1d,qgrid,qphon,rfdir,rfpert,rprimd,timrev,typat,ucvol,usepaw,xcccrc,xccc1d,xred,rhor,vxc)
4225 
4226  use m_cgtools,       only : dotprod_vn
4227  use m_atm2fft,       only : dfpt_atm2fft
4228  use m_dfpt_mkvxc,    only : dfpt_mkvxc, dfpt_mkvxc_noncoll
4229 
4230 !This section has been created automatically by the script Abilint (TD).
4231 !Do not modify the following lines by hand.
4232 #undef ABI_FUNC
4233 #define ABI_FUNC 'dfpt_dyxc1'
4234 !End of the abilint section
4235 
4236  implicit none
4237 
4238 !Arguments ------------------------------------
4239 !scalars
4240  integer,intent(in) :: ixc,mgfft,mpert,mqgrid,n1xccc,natom,nfft,nkxc,nspden,ntypat
4241  integer,intent(in) :: paral_kgb,timrev,usepaw
4242  real(dp),intent(in) :: gsqcut,ucvol
4243  type(pseudopotential_type),intent(in) :: psps
4244  type(MPI_type),intent(in) :: mpi_enreg
4245 !arrays
4246  integer,intent(in) :: atindx(natom),ngfft(18),rfdir(3),rfpert(mpert),typat(natom)
4247  real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc)
4248  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3)
4249  real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat),xcccrc(ntypat)
4250  real(dp),intent(in) :: xred(3,natom)
4251  integer,intent(out) :: blkflgfrx1(3,natom,3,natom)
4252  real(dp),intent(out) :: dyfrx1(2,3,natom,3,natom)
4253  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
4254 !optional
4255  real(dp),optional,intent(in) :: rhor(nfft,nspden)
4256  real(dp),optional,intent(in) :: vxc(nfft,nspden)
4257 
4258 !Local variables-------------------------------
4259 !scalars
4260  integer :: cplex,iat1,iatom1,iatom2,idir1,idir2,ierr,ifft,my_natom,comm_atom
4261  integer :: n1,n2,n3,n3xccc,nfftot,option,upperdir,optnc
4262  logical :: paral_atom
4263  real(dp) :: valuei,valuer
4264 !arrays
4265  integer,pointer :: my_atmtab(:)
4266  real(dp) :: tsec(2),gprimd_dummy(3,3)
4267  real(dp) :: dum_nhat(0)
4268  real(dp),allocatable :: rhor1(:,:),vxc10(:,:),xcccwk1(:),xcccwk2(:)
4269 ! *********************************************************************
4270 
4271  call timab(182,1,tsec)
4272 
4273  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
4274  nfftot=n1*n2*n3
4275 
4276 !Set up parallelism over atoms
4277  my_natom=mpi_enreg%my_natom
4278  my_atmtab=>mpi_enreg%my_atmtab
4279  comm_atom=mpi_enreg%comm_atom
4280  paral_atom=(my_natom/=natom)
4281 
4282 !Zero out the output arrays :
4283  blkflgfrx1(:,:,:,:)=0
4284  dyfrx1(:,:,:,:,:)=zero
4285 
4286  cplex=2-timrev ; n3xccc=nfft
4287  ABI_ALLOCATE(vxc10,(cplex*nfft,nspden))
4288 
4289 
4290 !Loop on the perturbation j1
4291  do iat1=1,my_natom
4292    iatom1=iat1; if(paral_atom)iatom1=my_atmtab(iat1)
4293    do idir1=1,3
4294 
4295 !    Compute the derivative of the core charge with respect to j1
4296      ABI_ALLOCATE(xcccwk1,(cplex*n3xccc))
4297 
4298 !    PAW or NC with nc_xccc_gspace: 1st-order core charge in reciprocal space
4299      if (usepaw==1 .or. psps%nc_xccc_gspace==1) then
4300        call dfpt_atm2fft(atindx,cplex,gmet,gprimd_dummy,gsqcut,idir1,iatom1,&
4301 &       mgfft,mqgrid,natom,1,nfft,ngfft,ntypat,ph1d,qgrid,&
4302 &       qphon,typat,ucvol,usepaw,xred,psps,pawtab,atmrhor1=xcccwk1,optn2_in=1)
4303 
4304 !      Norm-conserving psp: 1st-order core charge in real space
4305      else
4306        call dfpt_mkcore(cplex,idir1,iatom1,natom,ntypat,n1,n1xccc,&
4307 &       n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xcccwk1,xred)
4308      end if
4309 
4310 !    Compute the corresponding potential
4311      option=0
4312      ABI_ALLOCATE(rhor1,(cplex*nfft,nspden))
4313      rhor1=zero
4314 !FR SPr EB Non-collinear magnetism
4315      if (nspden==4.and.present(rhor).and.present(vxc)) then
4316        optnc=1
4317        call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,dum_nhat,0,dum_nhat,0,dum_nhat,0,nkxc,&
4318 &       nspden,n3xccc,optnc,option,paral_kgb,qphon,rhor,rhor1,rprimd,0,vxc,vxc10,xcccwk1)
4319      else
4320        call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,dum_nhat,0,dum_nhat,0,nkxc,&
4321 &       nspden,n3xccc,option,paral_kgb,qphon,rhor1,rprimd,0,vxc10,xcccwk1)
4322      end if
4323      ABI_DEALLOCATE(rhor1)
4324      ABI_DEALLOCATE(xcccwk1)
4325 
4326 !    vxc10 will couple with xcccwk2, that behaves like
4327 !    a total density (ispden=1). Only the spin-up + spin-down
4328 !    average of vxc10 is needed.
4329      if (nspden/=1)then
4330        do ifft=1,cplex*nfft
4331          vxc10(ifft,1)=(vxc10(ifft,1)+vxc10(ifft,2))*half
4332        end do
4333      end if
4334 
4335 !    Loop on the perturbation j2
4336      do iatom2=1,iatom1
4337        upperdir=3
4338        if(iatom1==iatom2)upperdir=idir1
4339        do idir2=1,upperdir
4340          if( (rfpert(iatom1)==1 .and. rfdir(idir1) == 1) .or. &
4341 &         (rfpert(iatom2)==1 .and. rfdir(idir2) == 1)    )then
4342 
4343 !          Compute the derivative of the core charge with respect to j2
4344            ABI_ALLOCATE(xcccwk2,(cplex*n3xccc))
4345 
4346 !          PAW or NC with nc_xccc_gspace: 1st-order core charge in reciprocal space
4347            if (usepaw==1 .or. psps%nc_xccc_gspace==1) then
4348              call dfpt_atm2fft(atindx,cplex,gmet,gprimd_dummy,gsqcut,idir2,iatom2,&
4349 &             mgfft,mqgrid,natom,1,nfft,ngfft,ntypat,ph1d,qgrid,&
4350 &             qphon,typat,ucvol,usepaw,xred,psps,pawtab,atmrhor1=xcccwk2,optn2_in=1)
4351 
4352 !            Norm-conserving psp: 1st-order core charge in real space
4353            else
4354              call dfpt_mkcore(cplex,idir2,iatom2,natom,ntypat,n1,n1xccc,&
4355 &             n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xcccwk2,xred)
4356            end if
4357 
4358 !          Get the matrix element j1,j2
4359 
4360            call dotprod_vn(cplex,xcccwk2,valuer,valuei,nfft,nfftot,1,2,vxc10,ucvol)
4361 
4362            ABI_DEALLOCATE(xcccwk2)
4363 
4364            dyfrx1(1,idir1,iatom1,idir2,iatom2)= valuer
4365            dyfrx1(2,idir1,iatom1,idir2,iatom2)= valuei
4366            dyfrx1(1,idir2,iatom2,idir1,iatom1)= valuer
4367            dyfrx1(2,idir2,iatom2,idir1,iatom1)=-valuei
4368            blkflgfrx1(idir1,iatom1,idir2,iatom2)=1
4369            blkflgfrx1(idir2,iatom2,idir1,iatom1)=1
4370          end if
4371        end do
4372      end do
4373    end do
4374  end do
4375 
4376  if (paral_atom) then
4377    call timab(48,1,tsec)
4378    call xmpi_sum(dyfrx1,comm_atom,ierr)
4379    call xmpi_sum(blkflgfrx1,comm_atom,ierr)
4380    call timab(48,2,tsec)
4381  end if
4382 
4383  ABI_DEALLOCATE(vxc10)
4384 
4385  call timab(182,2,tsec)
4386 
4387 end subroutine dfpt_dyxc1

ABINIT/dfpt_gatherdy [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_gatherdy

FUNCTION

 Sum (gather) the different parts of the 2nd-order matrix,
 to get the matrix of second-order derivatives (d2matr)
 Then, generates the dynamical matrix, not including the masses,
 but the correct non-cartesian coordinates ( => d2cart)

INPUTS

 becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only)
 berryopt=option for berry phase treatment
 blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical
  matrix has been calculated ; 0 otherwise )
 dyew(2,3,natom,3,natom)=Ewald part of the dyn.matrix
 dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf part of the dyn.matrix (except xc1)
 dyfrx1(2,3,natom,3,natom)=xc core correction (1) part of the frozen-wf
  part of the dynamical matrix.
 dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
 dyfr_nondiag=1 if dyfrwf is non diagonal with respect to atoms; 0 otherwise
 dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix
 d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
       second order derivatives
 d2nfr(2,3,mpert,3,mpert)=non-frozen wf part of the 2nd-order matr
 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
 eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor
 gprimd(3,3)=basis vector in the reciprocal space
 mband=maximum number of bands
 mpert =maximum number of ipert
 natom=number of atoms in unit cell
 ntypat=number of atom types
 outd2=option for the output of the 2nd-order matrix :
  if outd2=1, non-stationary part
  if outd2=2, stationary part.
 pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only)
 prtbbb=if 1, print the band-by-band decomposition, otherwise, prtbbb=0
 rfasr= (0=> no acoustic sum rule [asr] imposed), (1=> asr is imposed,
  in the democratic way for the effective charges),
 (2=> asr is imposed, in the aristocratic way for the effective
  charges)
 rfpert(mpert)=define the perturbations
 rprimd(3,3)=dimensional primitive translations (bohr)
 typat(natom)=integer label of each type of atom (1,2,...)
 ucvol=unit cell volume
 usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
 zion(ntypat)=charge corresponding to the atom type

OUTPUT

 carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
  2DTE matrix has been calculated correctly ; 0 otherwise )
 d2cart(2,3,mpert,3,mpert)=
  dynamical matrix, effective charges, dielectric tensor,....
  all in cartesian coordinates
 d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)=
  band by band decomposition of Born effective charges
  (calculated from phonon-type perturbation) in cartesian coordinates
 d2matr(2,3,mpert,3,mpert)=2nd-order matrix (masses non included,
  no cartesian coordinates : simply second derivatives)

PARENTS

      respfn

CHILDREN

      asria_calc,asria_corr,cart29,cart39,chneu9

SOURCE

3467 subroutine dfpt_gatherdy(becfrnl,berryopt,blkflg,carflg,dyew,dyfrwf,dyfrx1,&
3468 & dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,&
3469 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
3470 & gprimd,mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,prtbbb,&
3471 & rfasr,rfpert,rprimd,typat,ucvol,usevdw,zion)
3472 
3473 
3474 !This section has been created automatically by the script Abilint (TD).
3475 !Do not modify the following lines by hand.
3476 #undef ABI_FUNC
3477 #define ABI_FUNC 'dfpt_gatherdy'
3478 !End of the abilint section
3479 
3480  implicit none
3481 
3482 !Arguments -------------------------------
3483 !scalars
3484  integer,intent(in) :: berryopt,dyfr_cplex,dyfr_nondiag,mband,mpert,natom,ntypat,outd2
3485  integer,intent(in) :: pawbec,pawpiezo,prtbbb,rfasr,usevdw
3486  real(dp),intent(in) :: ucvol
3487 !arrays
3488  integer,intent(in) :: rfpert(mpert),typat(natom)
3489  integer,intent(inout) :: blkflg(3,mpert,3,mpert)
3490  integer,intent(out) :: carflg(3,mpert,3,mpert)
3491  real(dp),intent(in) :: becfrnl(3,natom,3*pawbec)
3492  real(dp),intent(in) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb)
3493  real(dp),intent(in) :: d2nfr(2,3,mpert,3,mpert),dyew(2,3,natom,3,natom)
3494  real(dp),intent(in) :: dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
3495  real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyvdw(2,3,natom,3,natom*usevdw)
3496  real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6)
3497  real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6)
3498  real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6)
3499  real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),gprimd(3,3)
3500  real(dp),intent(in) :: piezofrnl(6,3*pawpiezo),rprimd(3,3),zion(ntypat)
3501  real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
3502  real(dp),intent(out) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)
3503  real(dp),intent(out) :: d2matr(2,3,mpert,3,mpert)
3504 
3505 !Local variables -------------------------
3506 !scalars
3507  integer :: chneut,iband,iblok,idir,idir1,idir2,ii,ipert,ipert1,ipert2
3508  integer :: jj,nblok,selectz
3509  character(len=500) :: message
3510 !arrays
3511  integer :: flg1(3),flg2(3)
3512  real(dp) :: vec1(3),vec2(3)
3513 ! real(dp) :: ter(3,3) ! this variable appears commented out below
3514  real(dp),allocatable :: d2tmp(:,:,:,:,:),d2work(:,:,:,:,:),elfrtot(:,:)
3515 
3516 ! *********************************************************************
3517 
3518 
3519 !DEBUG
3520 !write(std_out,*)' dfpt_gatherdy : enter '
3521 !write(std_out,*)' outd2,mpert =',outd2,mpert
3522 !write(std_out,*)' blkflg(:,natom+2,:,natom+2)=',blkflg(:,natom+2,:,natom+2)
3523 !ENDDEBUG
3524 
3525  if(outd2/=3)then
3526 
3527 !  Initialise the 2nd-derivative matrix
3528    d2matr(:,:,:,:,:)=0.0_dp
3529 
3530 !  Add the non-frozen-part, the
3531 !  Ewald part and the xc1 part of the frozen-wf part
3532 !  Add the vdw part (if any)
3533    do ipert2=1,mpert
3534      do idir2=1,3
3535        do ipert1=1,mpert
3536          do idir1=1,3
3537            if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3538              do ii=1,2
3539                d2matr(ii,idir1,ipert1,idir2,ipert2)=&
3540 &               d2nfr(ii,idir1,ipert1,idir2,ipert2)
3541                if(ipert1<=natom .and. ipert2<=natom) then
3542                  d2matr(ii,idir1,ipert1,idir2,ipert2)=&
3543 &                 d2matr(ii,idir1,ipert1,idir2,ipert2)+&
3544 &                 dyew(ii,idir1,ipert1,idir2,ipert2)  +&
3545 &                 dyfrx1(ii,idir1,ipert1,idir2,ipert2)
3546                  if (usevdw==1) then
3547                    d2matr(ii,idir1,ipert1,idir2,ipert2)=&
3548 &                   d2matr(ii,idir1,ipert1,idir2,ipert2)+&
3549 &                   dyvdw(ii,idir1,ipert1,idir2,ipert2)
3550                  end if
3551                end if
3552              end do
3553            end if
3554          end do
3555        end do
3556      end do
3557    end do
3558 
3559 !  Add the frozen-wavefunction part
3560    if (dyfr_nondiag==0) then
3561      do ipert2=1,natom
3562        do idir2=1,3
3563          do idir1=1,3
3564            if( blkflg(idir1,ipert2,idir2,ipert2)==1 ) then
3565              d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)=&
3566 &             d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)&
3567 &             +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert2,1)
3568            end if
3569          end do
3570        end do
3571      end do
3572    else
3573      do ipert2=1,natom
3574        do ipert1=1,natom
3575          do idir2=1,3
3576            do idir1=1,3
3577              if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3578                d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)=&
3579 &               d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)&
3580 &               +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert1,ipert2)
3581              end if
3582            end do
3583          end do
3584        end do
3585      end do
3586    end if
3587 
3588 !  Add the frozen-wavefunction part of Born Effective Charges
3589    if (pawbec==1) then
3590      ipert2=natom+2
3591      do idir2=1,3            ! Direction of electric field
3592        do ipert1=1,natom     ! Atom
3593          do idir1=1,3        ! Direction of atom
3594            if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3595              d2matr(1,idir1,ipert1,idir2,ipert2)=&
3596 &             d2matr(1,idir1,ipert1,idir2,ipert2)+becfrnl(idir1,ipert1,idir2)
3597            end if
3598            if(blkflg(idir2,ipert2,idir1,ipert1)==1 ) then
3599              d2matr(1,idir2,ipert2,idir1,ipert1)=&
3600 &             d2matr(1,idir2,ipert2,idir1,ipert1)+becfrnl(idir1,ipert1,idir2)
3601            end if
3602          end do
3603        end do
3604      end do
3605    end if
3606 
3607 !  Section for piezoelectric tensor (from electric field response only for PAW)
3608    if(rfpert(natom+2)==1.and.pawpiezo==1) then
3609      ipert2=natom+2
3610      do idir2=1,3            ! Direction of electric field
3611        do ipert1=natom+3,natom+4     ! Strain
3612          do idir1=1,3
3613            ii=idir1+3*(ipert1-natom-3)
3614            if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3615              d2matr(1,idir1,ipert1,idir2,ipert2)=&
3616 &             d2nfr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir2)
3617            end if
3618          end do
3619        end do
3620      end do
3621    end if
3622 
3623 !  Section for strain perturbation
3624    if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then
3625 !    Make sure relevant columns of output are nulled
3626      d2matr(:,:,:,:,natom+3:natom+4)=0.0_dp
3627 !    Accumulate all frozen parts of the elastic tensor
3628      ABI_ALLOCATE(elfrtot,(6+3*natom,6))
3629      elfrtot(:,:)=elteew(:,:)+eltfrloc(:,:)+eltfrnl(:,:)+eltfrxc(:,:)
3630      elfrtot(1:6,1:6)=elfrtot(1:6,1:6)+eltcore(:,:)+eltfrhar(:,:)+eltfrkin(:,:)
3631      if (usevdw==1) elfrtot(:,:)=elfrtot(:,:)+eltvdw(:,:)
3632 
3633      do ipert2=natom+3,natom+4
3634        do idir2=1,3
3635 !        Internal strain components first
3636          do ipert1=1,natom
3637            do idir1=1,3
3638              if( blkflg(1,ipert1,idir2,ipert2)==1 ) then
3639                ii=idir1+6+3*(ipert1-1)
3640                jj=idir2+3*(ipert2-natom-3)
3641                d2matr(1,idir1,ipert1,idir2,ipert2)=&
3642 &               d2nfr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj)
3643              end if
3644            end do
3645          end do
3646 !        Now, electric field - strain mixed derivative (piezoelectric tensor)
3647          ipert1=natom+2
3648          do idir1=1,3
3649            if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3650              d2matr(1,idir1,ipert1,idir2,ipert2)=&
3651 &             d2nfr(1,idir1,ipert1,idir2,ipert2)
3652              if (pawpiezo==1) then
3653                ii=idir2+3*(ipert2-natom-3)
3654                d2matr(1,idir1,ipert1,idir2,ipert2)=&
3655 &               d2matr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir1)
3656              end if
3657            end if
3658          end do
3659 !        Now, strain-strain 2nd derivatives
3660          do ipert1=natom+3,natom+4
3661            do idir1=1,3
3662              if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3663                ii=idir1+3*(ipert1-natom-3)
3664                jj=idir2+3*(ipert2-natom-3)
3665                d2matr(1,idir1,ipert1,idir2,ipert2)=&
3666 &               d2nfr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj)
3667              end if
3668            end do
3669          end do
3670        end do
3671      end do
3672      ABI_DEALLOCATE(elfrtot)
3673    end if
3674 !  End section for strain perturbation
3675 
3676 !  The second-order matrix has been computed.
3677 
3678 !  Filter now components smaller in absolute value than 1.0d-20,
3679 !  for automatic testing reasons
3680    do ipert2=1,mpert
3681      do idir2=1,3
3682        do ipert1=1,mpert
3683          do idir1=1,3
3684            if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then
3685              do ii=1,2
3686                if( d2matr(ii,idir1,ipert1,idir2,ipert2)**2 < 1.0d-40)then
3687                  d2matr(ii,idir1,ipert1,idir2,ipert2)=zero
3688                end if
3689              end do
3690            end if
3691          end do
3692        end do
3693      end do
3694    end do
3695 
3696 !  DEBUG
3697 !  write(std_out,*)' d2matr '
3698 !  ipert2=natom+2
3699 !  do idir2=1,3
3700 !  ipert1=natom+2
3701 !  do idir1=1,3
3702 !  write(std_out,'(4i4,2es16.6)' )idir1,ipert1,idir2,ipert2,&
3703 !  &       d2matr(1,idir1,ipert1,idir2,ipert2),&
3704 !  &       d2matr(2,idir1,ipert1,idir2,ipert2)
3705 !  end do
3706 !  end do
3707 !  ENDDEBUG
3708 
3709 !  Cartesian coordinates transformation
3710    iblok=1 ; nblok=1
3711 
3712 !  In the case of finite electric field, the convention for the
3713 !  direction of the electric field perturbation was NOT the usual convention ...
3714 !  So, there is a transformation to the usual conventions
3715 !  to be done first ...
3716    if((berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 )  &
3717 &   .and. minval(abs(blkflg(:,natom+2,:,natom+2)))/=0)then   !!HONG  need to check for fixed D and E calculation
3718      if(minval(abs(blkflg(:,natom+2,:,natom+2)-1))/=0)then
3719        write(message,'(5a)')&
3720 &       '  In case of finite electric field, and electric field perturbation,',ch10,&
3721 &       '  the three directions for the perturbations must be treated.',ch10,&
3722 &       '  Action : set idir to 1 1 1, or forget about finite electric field.'
3723        MSG_ERROR(message)
3724      end if
3725      do ipert=1,mpert
3726        do idir=1,3
3727          do ii=1,2
3728            vec1(:)=d2matr(ii,idir,ipert,:,natom+2)
3729            flg1(:)=blkflg(idir,ipert,:,natom+2)
3730            call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2)
3731            d2matr(ii,idir,ipert,:,natom+2)=vec2(:)*two_pi
3732            blkflg(idir,ipert,:,natom+2)=flg2(:)
3733          end do
3734        end do
3735      end do
3736      do ipert=1,mpert
3737        do idir=1,3
3738          do ii=1,2
3739            vec1(:)=d2matr(ii,:,natom+2,idir,ipert)
3740            flg1(:)=blkflg(:,natom+2,idir,ipert)
3741            call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2)
3742            d2matr(ii,:,natom+2,idir,ipert)=vec2(:)*two_pi
3743            blkflg(:,natom+2,idir,ipert)=flg2(:)
3744          end do
3745        end do
3746      end do
3747 !    Also to be done, a change of sign, that I do not understand (XG071110)
3748 !    Perhaps due to d/dk replacing id/dk ? !
3749      d2matr(:,:,natom+2,:,natom+2)=-d2matr(:,:,natom+2,:,natom+2)
3750    end if
3751 
3752    call cart29(blkflg,d2matr,carflg,d2cart,&
3753 &   gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)
3754 
3755 !  Band by band decomposition of the Born effective charges
3756    if(prtbbb==1)then
3757      ABI_ALLOCATE(d2work,(2,3,mpert,3,mpert))
3758      ABI_ALLOCATE(d2tmp,(2,3,mpert,3,mpert))
3759      do iband=1,mband
3760        d2work(:,:,:,:,:)=0.0_dp
3761        d2tmp(:,:,:,:,:)=0.0_dp
3762        d2work(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband)
3763        call cart29(blkflg,d2work,carflg,d2tmp,&
3764 &       gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion)
3765 
3766 !      Remove the ionic part
3767        do ipert1=1,natom
3768          do idir1=1,3
3769            d2tmp(1,idir1,natom+2,idir1,ipert1) = &
3770 &           d2tmp(1,idir1,natom+2,idir1,ipert1) - zion(typat(ipert1))
3771          end do
3772        end do
3773 
3774        d2cart_bbb(:,:,:,:,iband,iband)=d2tmp(:,:,natom+2,:,:)
3775 
3776      end do
3777      ABI_DEALLOCATE(d2tmp)
3778      ABI_DEALLOCATE(d2work)
3779    end if ! prtbbb==1
3780 
3781 !
3782 !  Now, the cartesian elements are ready for output
3783 !  carflg give the information on their correctness
3784  end if
3785 
3786 !  Imposition of the ASR on the analytical part of the DynMat
3787 !  Assume that if rfasr/=0, the whole cartesian matrix is correct
3788  if(rfasr/=0)then
3789 
3790    ABI_ALLOCATE(d2work,(2,3,mpert,3,mpert))
3791    call asria_calc(rfasr,d2work,d2cart,mpert,natom)
3792 !  The following line imposes ASR:
3793    call asria_corr(rfasr,d2work,d2cart,mpert,natom)
3794 
3795    ABI_DEALLOCATE(d2work)
3796 
3797 !  Imposition of the ASR on the effective charges.
3798    if(rfpert(natom+2)==1)then
3799      chneut=rfasr
3800      selectz=0
3801      call chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion)
3802    end if
3803 
3804  end if
3805 
3806 !Additional operations on cartesian strain derivatives
3807  if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then
3808 !  Impose zero-net-force condition on internal strain tensor
3809    do ipert2=natom+3,natom+4
3810      do idir2=1,3
3811        vec1(:)=0.0_dp
3812        do ipert1=1,natom
3813          do idir1=1,3
3814            if(carflg(idir1,ipert1,idir2,ipert2)==1) then
3815              vec1(idir1)=vec1(idir1)+d2cart(1,idir1,ipert1,idir2,ipert2)
3816            end if
3817          end do
3818        end do
3819        vec1(:)=vec1(:)/dble(natom)
3820        do ipert1=1,natom
3821          do idir1=1,3
3822            if(carflg(idir1,ipert1,idir2,ipert2)==1) then
3823 !            Note minus sign to convert gradients to forces
3824              d2cart(1,idir1,ipert1,idir2,ipert2)=&
3825 &             -(d2cart(1,idir1,ipert1,idir2,ipert2)-vec1(idir1))
3826            end if
3827          end do
3828        end do
3829      end do
3830    end do
3831 !  Divide strain 2nd deriviative by ucvol to give elastic tensor
3832    do ipert2=natom+3,natom+4
3833      do idir2=1,3
3834        do ipert1=natom+3,natom+4
3835          do idir1=1,3
3836            if(carflg(idir1,ipert1,idir2,ipert2)==1) then
3837              d2cart(1,idir1,ipert1,idir2,ipert2)=&
3838 &             d2cart(1,idir1,ipert1,idir2,ipert2)/ucvol
3839            end if
3840          end do
3841        end do
3842      end do
3843    end do
3844  end if
3845 
3846 !calculate Born effective charges from electric field perturbation
3847 !do ipert1=1,natom
3848 !ter(:,:)=zero
3849 !do idir1=1,3
3850 !do ii=1,3
3851 !do jj=1,3
3852 !if(abs(gprimd(idir1,ii))>1.0d-10)then
3853 !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,idir1,natom+2,jj,ipert1)*gprimd(jj,ii)
3854 !endif
3855 !enddo
3856 !enddo
3857 !add zion to bec
3858 !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1))
3859 !enddo
3860 !d2cart(1,:,ipert1,:,natom+2)=ter(:,:)
3861 !enddo
3862 !carflg(:,1:natom,:,natom+2)=1
3863 
3864 !Born effective charges from phonon perturbation
3865 !do ipert1=1,natom
3866 !ter(:,:)=zero
3867 !do idir1=1,3
3868 !do ii=1,3
3869 !do jj=1,3
3870 !if(abs(gprimd(idir1,ii))>1.0d-10)then
3871 !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,jj,ipert1,idir1,natom+2)*gprimd(jj,ii)
3872 !endif
3873 !enddo
3874 !enddo
3875 !! add zion to bec
3876 !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1))
3877 !enddo
3878 !d2cart(1,:,natom+2,:,ipert1)=ter(:,:)
3879 !enddo
3880 !carflg(:,natom+2,:,1:natom)=1
3881 
3882 
3883 !!Dielectric constant
3884 !do ii=1,3
3885 !do jj=1,3
3886 !ter(ii,jj)=d2nfr(1,ii,natom+2,jj,natom+2)
3887 !end do
3888 !end do
3889 !ter(:,:)=pi*four*ter(:,:)/ucvol
3890 !
3891 !do ii=1,3
3892 !ter(ii,ii)=ter(ii,ii)+one
3893 !end do
3894 !d2cart(1,:,natom+2,:,natom+2)=ter(:,:)
3895 !carflg(:,natom+2,:,1:natom+2)=1
3896 
3897 !DEBUG
3898 !Debugging, but needed to make the stuff work on the IBM Dirac ? !
3899 !write(std_out,*)' d2cart '
3900 !ipert2=natom+2
3901 !do idir2=1,3
3902 !ipert1=natom+2
3903 !do idir1=1,3
3904 !write(std_out,'(5i4,2d20.10)' )idir1,ipert1,idir2,ipert2,&
3905 !&      carflg(idir1,ipert1,idir2,ipert2),&
3906 !&      d2cart(1,idir1,ipert1,idir2,ipert2),&
3907 !&      d2cart(2,idir1,ipert1,idir2,ipert2)
3908 !end do
3909 !end do
3910 !ENDDEBUG
3911 
3912 !DEBUG
3913 !write(std_out,*)' dfpt_gatherdy : exit '
3914 !ENDDEBUG
3915 
3916 end subroutine dfpt_gatherdy

ABINIT/m_respfn_driver [ Modules ]

[ Top ] [ Modules ]

NAME

  m_respfn_driver

FUNCTION

  Subdriver for DFPT calculations.

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XG, DRH, MT, MKV)
  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_respfn_driver
 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_exit
 37  use m_wffile
 38  use m_errors
 39  use m_ebands
 40  use m_results_respfn
 41  use m_hdr
 42  use m_crystal
 43  use m_xcdata
 44 
 45  use m_dtfil,       only : status
 46  use m_time,        only : timab
 47  use m_fstrings,    only : strcat
 48  use m_symtk,       only : matr3inv, littlegroup_q, symmetrize_xred
 49  use m_fft,         only : zerosym, fourdp
 50  use m_kpts,        only : symkchk
 51  use m_geometry,    only : irreducible_set_pert
 52  use m_dynmat,      only : chkph3, d2sym3, q0dy3_apply, q0dy3_calc, wings3, dfpt_phfrq, sytens, dfpt_prtph, &
 53                            asria_calc, asria_corr, cart29, cart39, chneu9, dfpt_sydy
 54  use m_ddb,         only : DDB_VERSION
 55  use m_ddb_hdr,     only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write
 56  use m_ddb_interpolate, only : outddbnc
 57  use m_occ,         only : newocc
 58  use m_efmas,       only : efmasdeg_free_array, efmasval_free_array
 59  use m_wfk,         only : wfk_read_eigenvalues
 60  use m_ioarr,       only : read_rhor
 61  use m_pawang,      only : pawang_type
 62  use m_pawrad,      only : pawrad_type
 63  use m_pawtab,      only : pawtab_type, pawtab_get_lsize
 64  use m_paw_an,      only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify
 65  use m_paw_ij,      only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify
 66  use m_pawfgrtab,   only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
 67  use m_pawrhoij,    only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, &
 68                            pawrhoij_bcast, pawrhoij_nullify, pawrhoij_get_nspden
 69  use m_pawdij,      only : pawdij, symdij
 70  use m_pawfgr,      only : pawfgr_type, pawfgr_init, pawfgr_destroy
 71  use m_paw_finegrid,only : pawexpiqr
 72  use m_pawxc,       only : pawxc_get_nkxc
 73  use m_paw_dmft,    only : paw_dmft_type
 74  use m_paw_sphharm, only : setsym_ylm
 75  use m_paw_nhat,    only : nhatgrid,pawmknhat
 76  use m_paw_tools,   only : chkpawovlp
 77  use m_paw_denpot,  only : pawdenpot
 78  use m_paw_init,    only : pawinit,paw_gencond
 79  use m_kg,          only : getcut, getph, kpgio
 80  use m_eig2d,       only : eig2tot, elph2_fanddw
 81  use m_inwffil,     only : inwffil
 82  use m_spacepar,    only : hartre, setsym
 83  use m_mkrho,       only : mkrho
 84  use m_vdw_dftd2,   only : vdw_dftd2
 85  use m_vdw_dftd3,   only : vdw_dftd3
 86  use m_initylmg,    only : initylmg
 87  use m_pspini,      only : pspini
 88  use m_atm2fft,     only : atm2fft
 89  use m_dfpt_loopert,only : dfpt_looppert, eigen_meandege
 90  use m_rhotoxc,     only : rhotoxc
 91  use m_drivexc,     only : check_kxc
 92  use m_mklocl,      only : mklocl, mklocl_recipspace
 93  use m_common,      only : setup1, prteigrs
 94  use m_fourier_interpol, only : transgrid
 95  use m_paral_atom,     only : get_my_atmtab, free_my_atmtab
 96  use m_paw_occupancies, only : initrhoij
 97  use m_paw_correlations,only : pawpuxinit
 98  use m_mkcore,     only : mkcore, dfpt_mkcore
 99  use m_dfpt_elt,   only : dfpt_eltfrxc, dfpt_eltfrloc, dfpt_eltfrkin, dfpt_eltfrhar, elt_ewald, dfpt_ewald
100  use m_d2frnl,     only : d2frnl
101 
102 #if defined HAVE_GPU_CUDA
103  use m_alloc_hamilt_gpu, only : alloc_hamilt_gpu, dealloc_hamilt_gpu
104 #endif
105 
106  implicit none
107 
108  private

m_respfn_driver/respfn [ Functions ]

[ Top ] [ m_respfn_driver ] [ Functions ]

NAME

 respfn

FUNCTION

 Primary routine for conducting DFT calculations of Response functions.

INPUTS

  codvsn=code version
  cpui=initial cpu time
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
   | mband=maximum number of bands
   | mgfft=maximum single fft dimension
   | mkmem=Number of k points treated by this node
   | mpw=maximum number of planewaves in basis sphere (large number)
   | natom=number of atoms in unit cell
   | nfft=(effective) number of FFT grid points (for this processor)
   | nkpt=number of k points
   | nspden=number of spin-density components
   | nsppol=number of channels for spin-polarization (1 or 2)
   | nsym=number of symmetry elements in space group
  mkmems(3)=array containing the tree values of mkmem (see above) (k-GS, k+q-GS and RF)
  mpi_enreg=information about MPI parallelization
  npwtot(nkpt)=number of planewaves in basis and boundary at each k point
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  etotal=total energy (sum of 7 or 8 contributions) (hartree)

SIDE EFFECTS

  iexit=index of "exit" on first line of file (0 if not found)
  occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
    Occupations number may have been read from a previous dataset...
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
    Some dimensions in pawrad have been set in driver.f
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
    Some dimensions in pawtab have been set in driver.f
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
    Before entering the first time in respfn, a significant part of psps
    has been initialized: the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,
    mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
    All the remaining components of psps are to be initialized in the call
    to pspini.  The next time the code enters respfn, psps might be identical
    to the one of the previous dtset, in which case, no reinitialisation
    is scheduled in pspini.f .
  results_respfn <type(results_respfn_type)>=stores some results of respfn calls

NOTES

 USE OF FFT GRIDS:
 =================
 In case of PAW:
 ---------------
    Two FFT grids are used:
    - A "coarse" FFT grid (defined by ecut)
      for the application of the Hamiltonian on the plane waves basis.
      It is defined by nfft, ngfft, mgfft, ...
      Hamiltonian, wave-functions, density related to WFs (rhor here), ...
      are expressed on this grid.
    - A "fine" FFT grid (defined) by ecutdg)
      for the computation of the density inside PAW spheres.
      It is defined by nfftf, ngfftf, mgfftf, ...
      Total density, potentials, ...
      are expressed on this grid.
 In case of norm-conserving:
 ---------------------------
    - Only the usual FFT grid (defined by ecut) is used.
      It is defined by nfft, ngfft, mgfft, ...
      For compatibility reasons, (nfftf,ngfftf,mgfftf)
      are set equal to (nfft,ngfft,mgfft) in that case.

PARENTS

      driver

CHILDREN

      alloc_hamilt_gpu,atm2fft,check_kxc,chkpawovlp,chkph3,crystal_free
      crystal_init,d2frnl,d2sym3,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write
      dealloc_hamilt_gpu,dfpt_dyfro,dfpt_dyout,dfpt_dyxc1,dfpt_eltfrhar
      dfpt_eltfrkin,dfpt_eltfrloc,dfpt_eltfrxc,dfpt_ewald,dfpt_gatherdy
      dfpt_looppert,dfpt_phfrq,dfpt_prtph,ebands_free,efmasdeg_free_array
      efmasval_free_array,eig2tot,eigen_meandege,elph2_fanddw,elt_ewald
      exit_check,fourdp,getcut,getph,hartre,hdr_free,hdr_init,hdr_update
      initrhoij,initylmg,inwffil,irreducible_set_pert,kpgio,littlegroup_q
      matr3inv,mkcore,mklocl,mkrho,newocc,nhatgrid,outddbnc,paw_an_free
      paw_an_init,paw_an_nullify,paw_gencond,paw_ij_free,paw_ij_init
      paw_ij_nullify,pawdenpot,pawdij,pawexpiqr,pawfgr_destroy,pawfgr_init
      pawfgrtab_free,pawfgrtab_init,pawinit,pawmknhat,pawpuxinit
      pawrhoij_alloc,pawrhoij_bcast,pawrhoij_copy,pawrhoij_free
      pawrhoij_nullify,pawtab_get_lsize,prteigrs,pspini,q0dy3_apply
      q0dy3_calc,read_rhor,rhotoxc,setsym,setsym_ylm,setup1,status,symdij
      symmetrize_xred,sytens,timab,transgrid,vdw_dftd2,vdw_dftd3,wffclose
      wings3,wrtloctens,wrtout,xcdata_init,xmpi_bcast

SOURCE

 213 subroutine respfn(codvsn,cpui,dtfil,dtset,etotal,iexit,&
 214 &  mkmems,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,results_respfn,xred)
 215 
 216 
 217 !This section has been created automatically by the script Abilint (TD).
 218 !Do not modify the following lines by hand.
 219 #undef ABI_FUNC
 220 #define ABI_FUNC 'respfn'
 221 !End of the abilint section
 222 
 223  implicit none
 224 
 225 !Arguments ------------------------------------
 226  integer,intent(inout) :: iexit
 227  real(dp),intent(in) :: cpui
 228  real(dp),intent(inout) :: etotal !vz_i
 229  character(len=6),intent(in) :: codvsn
 230  type(MPI_type),intent(inout) :: mpi_enreg
 231  type(datafiles_type),intent(in) :: dtfil
 232  type(dataset_type),intent(in) :: dtset
 233  type(pawang_type),intent(inout) :: pawang
 234  type(pseudopotential_type),intent(inout) :: psps
 235  integer,intent(in) :: mkmems(3)
 236  integer,intent(inout) :: npwtot(dtset%nkpt)
 237  real(dp),intent(inout) :: xred(3,dtset%natom)
 238  real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 239  type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
 240  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 241  type(results_respfn_type),intent(inout) :: results_respfn
 242 
 243 !Local variables-------------------------------
 244  integer,parameter :: formeig=0,level=10
 245  integer,parameter :: response=1,syuse=0,master=0,cplex1=1
 246  integer :: nk3xc
 247  integer :: analyt,ask_accurate,band_index,bantot,bdeigrf,coredens_method,cplex
 248  integer :: dim_eig2nkq,dim_eigbrd,dyfr_cplex,dyfr_nondiag,gnt_option
 249  integer :: gscase,has_dijnd,has_diju,has_kxc,iatom,iatom_tot,iband,idir,ider,ierr,ifft,ii,ikpt,indx
 250  integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
 251  integer :: initialized,ipert,ipert2,ireadwf0,iscf,iscf_eff,ispden,isppol
 252  integer :: itypat,izero,mcg,me,mgfftf,mk1mem,mkqmem,mpert,mu
 253  integer :: my_natom,n1,natom,n3xccc,nband_k,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim
 254  integer :: nkpt_eff,nkpt_max,nkpt_rbz,nkxc,nkxc1,nspden_rhoij,ntypat,nzlmopt,openexit
 255  integer :: optcut,option,optgr0,optgr1,optgr2,optorth,optrad
 256  integer :: optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv
 257  integer :: outd2,pawbec,pawpiezo,prtbbb,psp_gencond,qzero,rdwr,rdwrpaw
 258  integer :: rfasr,rfddk,rfelfd,rfphon,rfstrs,rfuser,rf2_dkdk,rf2_dkde,rfmagn
 259  integer :: spaceworld,sumg0,sz1,sz2,tim_mkrho,timrev,usecprj,usevdw
 260  integer :: usexcnhat,use_sym,vloc_method,zero_by_symm
 261  logical :: has_full_piezo,has_allddk,paral_atom,qeq0,use_nhat_gga,call_pawinit,non_magnetic_xc
 262  real(dp) :: boxcut,compch_fft,compch_sph,cpus,ecore,ecut_eff,ecutdg_eff,ecutf
 263  real(dp) :: eei,eew,ehart,eii,ek,enl,entropy,enxc
 264  real(dp) :: epaw,epawdc,etot,evdw,fermie,gsqcut,gsqcut_eff,gsqcutc_eff,qphnrm,residm
 265  real(dp) :: ucvol,vxcavg
 266  character(len=fnlen) :: dscrpt
 267  character(len=fnlen) :: filename
 268  character(len=500) :: message
 269  type(ebands_t) :: bstruct
 270  type(hdr_type) :: hdr,hdr_fine,hdr0,hdr_den
 271  type(ddb_hdr_type) :: ddb_hdr
 272  type(paw_dmft_type) :: paw_dmft
 273  type(pawfgr_type) :: pawfgr
 274  type(wffile_type) :: wffgs,wfftgs
 275  type(wvl_data) :: wvl
 276  type(crystal_t) :: Crystal
 277  type(xcdata_type) :: xcdata
 278  integer :: ddkfil(3),ngfft(18),ngfftf(18),rfdir(3),rf2_dirs_from_rfpert_nl(3,3)
 279  integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:),blkflgfrx1(:,:,:,:),blkflg1(:,:,:,:)
 280  integer,allocatable :: blkflg2(:,:,:,:),carflg(:,:,:,:),clflg(:,:),indsym(:,:,:)
 281  integer,allocatable :: irrzon(:,:,:),kg(:,:),l_size_atm(:),nattyp(:),npwarr(:)
 282  integer,allocatable :: pertsy(:,:),rfpert(:),rfpert_nl(:,:,:,:,:,:),symq(:,:,:),symrec(:,:,:)
 283  real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0)
 284  real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0)
 285  real(dp) :: dummy6(6),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3),qphon(3)
 286  real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),strn_dummy6(6),strv_dummy6(6),strsxc(6),tsec(2)
 287  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
 288  real(dp),allocatable :: becfrnl(:,:,:),cg(:,:),d2bbb(:,:,:,:,:,:),d2cart(:,:,:,:,:)
 289  real(dp),allocatable :: d2cart_bbb(:,:,:,:,:,:),d2eig0(:,:,:,:,:)
 290  real(dp),allocatable :: d2k0(:,:,:,:,:),d2lo(:,:,:,:,:),d2loc0(:,:,:,:,:)
 291  real(dp),allocatable :: d2matr(:,:,:,:,:),d2nfr(:,:,:,:,:),d2nl(:,:,:,:,:),d2ovl(:,:,:,:,:)
 292  real(dp),allocatable :: d2nl0(:,:,:,:,:),d2nl1(:,:,:,:,:),d2tmp(:,:,:,:,:),d2vn(:,:,:,:,:)
 293  real(dp),allocatable :: displ(:),doccde(:)
 294  real(dp),allocatable :: dyew(:,:,:,:,:),dyewq0(:,:,:),dyfrlo(:,:,:),dyfrlo_indx(:,:,:)
 295  real(dp),allocatable :: dyfrnl(:,:,:,:,:),dyfrwf(:,:,:,:,:),dyfrx1(:,:,:,:,:),dyvdw(:,:,:,:,:)
 296  real(dp),allocatable :: dyfrx2(:,:,:),eigen0(:),eigval(:),eigvec(:)
 297  real(dp),allocatable :: eig2nkq(:,:,:,:,:,:,:),eigbrd(:,:,:,:,:,:,:)
 298  real(dp),allocatable :: eigen_fan(:),eigen_ddw(:),eigen_fanddw(:)
 299  real(dp),allocatable :: eigen_fan_mean(:),eigen_ddw_mean(:)
 300  real(dp),allocatable :: eltcore(:,:),elteew(:,:),eltfrhar(:,:),eltfrkin(:,:)
 301  real(dp),allocatable :: eltfrloc(:,:),eltfrnl(:,:),eltfrxc(:,:),eltvdw(:,:),grtn_indx(:,:)
 302  real(dp),allocatable :: grxc(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:)
 303  real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phfrq(:),phnons(:,:,:),piezofrnl(:,:)
 304  real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:)
 305  real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:)
 306  real(dp),allocatable :: vxc(:,:),work(:),xccc3d(:),ylm(:,:),ylmgr(:,:,:)
 307  real(dp),pointer :: eigenq_fine(:,:,:),eigen1_pert(:,:,:)
 308  real(dp),allocatable :: eigen0_pert(:),eigenq_pert(:),occ_rbz_pert(:)
 309  type(efmasdeg_type),allocatable :: efmasdeg(:)
 310  type(efmasval_type),allocatable :: efmasval(:,:)
 311  type(paw_an_type),allocatable :: paw_an(:)
 312  type(paw_ij_type),allocatable :: paw_ij(:)
 313  type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
 314  type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:)
 315 ! ***********************************************************************
 316 
 317  DBG_ENTER("COLL")
 318 
 319  call timab(132,1,tsec)
 320  call timab(133,1,tsec)
 321 
 322  call status(0,dtfil%filstat,iexit,level,'init          ')
 323 
 324 ! Initialise non_magnetic_xc for rhohxc
 325  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
 326 
 327 !Some data for parallelism
 328  nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1
 329  my_natom=mpi_enreg%my_natom
 330  paral_atom=(my_natom/=dtset%natom)
 331 !Define FFT grid(s) sizes (be careful !)
 332 !See NOTES in the comments at the beginning of this file.
 333  call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
 334 
 335 !Structured debugging if dtset%prtvol==-level
 336  if(dtset%prtvol==-level)then
 337    write(message,'(80a,a,a)')  ('=',ii=1,80),ch10,' respfn : enter , debug mode '
 338    call wrtout(std_out,message,'COLL')
 339  end if
 340 
 341 !Option input variables
 342  iscf=dtset%iscf
 343 
 344 !Respfn input variables
 345  rfasr=dtset%rfasr   ; rfdir(1:3)=dtset%rfdir(1:3)
 346  rfddk=dtset%rfddk   ; rfelfd=dtset%rfelfd ; rfmagn=dtset%rfmagn
 347  rfphon=dtset%rfphon ; rfstrs=dtset%rfstrs
 348  rfuser=dtset%rfuser ; rf2_dkdk=dtset%rf2_dkdk ; rf2_dkde=dtset%rf2_dkde
 349 
 350  pawbec=0  ; if(psps%usepaw==1.and.(rfphon==1.or.(rfelfd==1.or.rfelfd==3))) pawbec=1
 351  pawpiezo=0; if(psps%usepaw==1.and.(rfstrs/=0.or.(rfelfd==1.or.rfelfd==3))) pawpiezo=1
 352 !AM 10152015 -- WARNING --- the full calculation of the piezoelectric tensor
 353 !from electric field perturbation is only available
 354 !if nsym/=1 (strain perturbation is not symmetrized):
 355  has_full_piezo=.False. ; if(pawpiezo==1.and.dtset%nsym==1)  has_full_piezo=.True.
 356  usevdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) usevdw=1
 357 !mkmem variables (mkmem is already argument)
 358  mkqmem=mkmems(2) ; mk1mem=mkmems(3)
 359 
 360  ntypat=psps%ntypat
 361  natom=dtset%natom
 362  nfftot=product(ngfft(1:3))
 363  nfftotf=product(ngfftf(1:3))
 364 
 365  call status(0,dtfil%filstat,iexit,level,'call setup1   ')
 366 
 367 !LIKELY TO BE TAKEN AWAY
 368  initialized=0
 369  ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero
 370  eii=zero ; eew=zero ; ecore=zero
 371 
 372 !Set up for iterations
 373  call setup1(dtset%acell_orig(1:3,1),bantot,dtset,&
 374 & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
 375 & natom,ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
 376 & response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw)
 377 
 378 !In some cases (e.g. getcell/=0), the plane wave vectors have
 379 ! to be generated from the original simulation cell
 380  rprimd_for_kg=rprimd
 381  if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=dtset%rprimd_orig(:,:,1)
 382  call matr3inv(rprimd_for_kg,gprimd_for_kg)
 383  gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg)
 384 
 385 !Define the set of admitted perturbations
 386  mpert=natom+7
 387  if (rf2_dkdk>0.or.rf2_dkde>0) mpert=natom+11
 388 
 389 !Initialize the list of perturbations rfpert
 390  ABI_ALLOCATE(rfpert,(mpert))
 391  rfpert(:)=0
 392  if(rfphon==1)rfpert(dtset%rfatpol(1):dtset%rfatpol(2))=1
 393 
 394  if(rfddk==1)rfpert(natom+1)=1
 395  if(rfddk==2)rfpert(natom+6)=1
 396 
 397  if(rf2_dkdk/=0)rfpert(natom+10)=1
 398  if(rf2_dkde/=0)rfpert(natom+11)=1
 399 
 400  if(rfelfd==1.or.rfelfd==2)rfpert(natom+1)=1
 401  if(rfelfd==1.or.rfelfd==3)rfpert(natom+2)=1
 402 
 403  if(rfstrs==1.or.rfstrs==3)rfpert(natom+3)=1
 404  if(rfstrs==2.or.rfstrs==3)rfpert(natom+4)=1
 405 
 406  if(rfuser==1.or.rfuser==3)rfpert(natom+6)=1
 407  if(rfuser==2.or.rfuser==3)rfpert(natom+7)=1
 408 
 409  if(rfmagn==1) rfpert(natom+5)=1
 410 
 411  qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14)
 412 
 413 !Init spaceworld
 414  spaceworld=mpi_enreg%comm_cell
 415  me = xmpi_comm_rank(spaceworld)
 416 
 417 !Set up the basis sphere of planewaves
 418  ABI_ALLOCATE(kg,(3,dtset%mpw*dtset%mkmem))
 419  ABI_ALLOCATE(npwarr,(dtset%nkpt))
 420  call status(0,dtfil%filstat,iexit,level,'call kpgio(1) ')
 421  call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg,&
 422 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,&
 423 & dtset%nsppol)
 424 
 425 !Set up the Ylm for each k point
 426  ABI_ALLOCATE(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
 427  if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) then
 428    ABI_ALLOCATE(ylmgr,(dtset%mpw*dtset%mkmem,9,psps%mpsang*psps%mpsang*psps%useylm))
 429  else
 430    ABI_ALLOCATE(ylmgr,(0,0,psps%useylm))
 431  end if
 432  if (psps%useylm==1) then
 433    call status(0,dtfil%filstat,iexit,level,'call initylmg ')
 434    option=0
 435    if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) option=2
 436    call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,&
 437 &   npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
 438  end if
 439 
 440  call timab(133,2,tsec)
 441  call timab(134,1,tsec)
 442 
 443 !Open and read pseudopotential files
 444  call status(0,dtfil%filstat,iexit,level,'call pspini(1)')
 445  call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,&
 446 & psps,rprimd,comm_mpi=mpi_enreg%comm_cell)
 447 
 448  call timab(134,2,tsec)
 449  call timab(135,1,tsec)
 450 
 451 !Initialize band structure datatype
 452  bstruct = ebands_from_dtset(dtset, npwarr)
 453 
 454 !Initialize PAW atomic occupancies
 455  if (psps%usepaw==1) then
 456    ABI_DATATYPE_ALLOCATE(pawrhoij,(my_natom))
 457    call pawrhoij_nullify(pawrhoij)
 458    call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, &
 459 &   my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,&
 460 &   pawrhoij,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,&
 461 &   comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
 462  else
 463    ABI_DATATYPE_ALLOCATE(pawrhoij,(0))
 464  end if
 465 
 466 !Initialize header
 467  gscase=0
 468  call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr, &
 469 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
 470 
 471 !Update header, with evolving variables, when available
 472 !Here, rprimd, xred and occ are available
 473  etot=hdr%etot ; fermie=hdr%fermie ; residm=hdr%residm
 474 !If parallelism over atom, hdr is distributed
 475  call hdr_update(hdr,bantot,etot,fermie,&
 476 & residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), &
 477 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
 478 
 479 !Clean band structure datatype (should use it more in the future !)
 480  call ebands_free(bstruct)
 481 
 482  call status(0,dtfil%filstat,iexit,level,'call inwffil(1')
 483 
 484 !Initialize wavefunction files and wavefunctions.
 485  ireadwf0=1
 486 
 487  mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
 488  ABI_STAT_ALLOCATE(cg,(2,mcg), ierr)
 489  ABI_CHECK(ierr==0, "out-of-memory in cg")
 490 
 491  ABI_ALLOCATE(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol))
 492  eigen0(:)=zero ; ask_accurate=1
 493  optorth=0
 494 
 495  hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors
 496  call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,&
 497 & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,&
 498 & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,&
 499 & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,&
 500 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
 501 & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl)
 502  hdr%rprimd=rprimd
 503 
 504 !Close wffgs, if it was ever opened (in inwffil)
 505  if (ireadwf0==1) then
 506    call WffClose(wffgs,ierr)
 507  end if
 508 
 509  if (psps%usepaw==1.and.ireadwf0==1) then
 510 !  if parallelism, pawrhoij is distributed, hdr%pawrhoij is not
 511    call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,&
 512 &   mpi_atmtab=mpi_enreg%my_atmtab)
 513  end if
 514 
 515  call timab(135,2,tsec)
 516  call timab(136,1,tsec)
 517 
 518 !Report on eigen0 values   ! Should use prteigrs.F90
 519  write(message, '(a,a)' )
 520  call wrtout(std_out,ch10//' respfn : eigen0 array','COLL')
 521  nkpt_eff=dtset%nkpt
 522  if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max
 523  band_index=0
 524  do isppol=1,dtset%nsppol
 525    do ikpt=1,dtset%nkpt
 526      nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 527      if(ikpt<=nkpt_eff)then
 528        write(message, '(a,i2,a,i5)' )'  isppol=',isppol,', k point number',ikpt
 529        call wrtout(std_out,message,'COLL')
 530        do iband=1,nband_k,4
 531          write(message, '(a,4es16.6)')'  ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index)
 532          call wrtout(std_out,message,'COLL')
 533        end do
 534      else if(ikpt==nkpt_eff+1)then
 535        write(message,'(a,a)' )'  respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10
 536        call wrtout(std_out,message,'COLL')
 537      end if
 538      band_index=band_index+nband_k
 539    end do
 540  end do
 541 
 542 !Allocation for forces and atomic positions (should be taken away, also argument ... )
 543  ABI_ALLOCATE(grxc,(3,natom))
 544 
 545 !Do symmetry stuff
 546  ABI_ALLOCATE(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 547  ABI_ALLOCATE(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 548  ABI_ALLOCATE(indsym,(4,dtset%nsym,natom))
 549  ABI_ALLOCATE(symrec,(3,3,dtset%nsym))
 550  irrzon=0;indsym=0;symrec=0;phnons=zero
 551 !If the density is to be computed by mkrho, need irrzon and phnons
 552  iscf_eff=0;if(dtset%getden==0)iscf_eff=1
 553  call setsym(indsym,irrzon,iscf_eff,natom,&
 554 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
 555 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
 556 
 557 !Symmetrize atomic coordinates over space group elements:
 558  call symmetrize_xred(indsym,natom,dtset%nsym,dtset%symrel,dtset%tnons,xred)
 559 
 560 !Examine the symmetries of the q wavevector
 561  ABI_ALLOCATE(symq,(4,2,dtset%nsym))
 562  timrev=1
 563 
 564 ! By default use symmetries.
 565  use_sym = 1
 566  if (dtset%prtgkk == 1)then
 567    use_sym = 0
 568    call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym)
 569  else
 570    call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol)
 571  end if
 572 
 573 !Generate an index table of atoms, in order for them to be used
 574 !type after type.
 575  ABI_ALLOCATE(atindx,(natom))
 576  ABI_ALLOCATE(atindx1,(natom))
 577  ABI_ALLOCATE(nattyp,(ntypat))
 578  indx=1
 579  do itypat=1,ntypat
 580    nattyp(itypat)=0
 581    do iatom=1,natom
 582      if(dtset%typat(iatom)==itypat)then
 583        atindx(iatom)=indx
 584        atindx1(indx)=iatom
 585        indx=indx+1
 586        nattyp(itypat)=nattyp(itypat)+1
 587      end if
 588    end do
 589  end do
 590 
 591 !Here allocation of GPU for fft calculations
 592 #if defined HAVE_GPU_CUDA
 593  if (dtset%use_gpu_cuda==1) then
 594    call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,0,psps,dtset%use_gpu_cuda)
 595  end if
 596 #endif
 597 
 598 !Compute structure factor phases for current atomic pos:
 599  ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*natom))
 600  ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*natom))
 601  call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
 602 
 603  if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
 604    call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
 605  else
 606    ph1df(:,:)=ph1d(:,:)
 607  end if
 608 
 609 !Compute occupation numbers and fermi energy, in case occupation scheme is metallic.
 610  ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
 611  if( dtset%occopt>=3.and.dtset%occopt<=8 ) then
 612 
 613    call newocc(doccde,eigen0,entropy,fermie,dtset%spinmagntarget,dtset%mband,dtset%nband,&
 614 &   dtset%nelect,dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,dtset%occopt,dtset%prtvol,dtset%stmbias,&
 615 &   dtset%tphysel,dtset%tsmear,dtset%wtk)
 616 
 617 !  Update fermie and occ
 618    etot=hdr%etot ; residm=hdr%residm
 619    call hdr_update(hdr,bantot,etot,fermie,residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1))
 620 
 621  else
 622 !  doccde is irrelevant in this case
 623    doccde(:)=zero
 624  end if
 625 
 626 !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx
 627  ecutf=dtset%ecut
 628  if (psps%usepaw==1) then
 629    ecutf=dtset%pawecutdg
 630    call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:','COLL')
 631  end if
 632 
 633  call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
 634 
 635 !PAW: 1- Initialize values for several arrays depending only on atomic data
 636 !2- Check overlap
 637 !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed)
 638 !4- Allocate PAW specific arrays
 639 !5- Compute perturbed local potential inside spheres
 640 !6- Eventually open temporary storage files
 641  if(psps%usepaw==1) then
 642 !  1-Initialize values for several arrays depending only on atomic data
 643    gnt_option=1
 644    if (dtset%pawxcdev==2.or.dtset%rfphon/=0.or.dtset%rfstrs/=0.or.dtset%rfelfd==1.or.&
 645    dtset%rfelfd==3.or.dtset%rf2_dkde==1) gnt_option=2
 646 
 647    ! Test if we have to call pawinit
 648    call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
 649 
 650    if (psp_gencond==1.or.call_pawinit) then
 651 !    Some gen-cond have to be added...
 652      call timab(553,1,tsec)
 653      call pawinit(gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,&
 654 &     psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,&
 655 &     pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%xclevel,dtset%usepotzero)
 656      call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,&
 657 &     rprimd,symrec,pawang%zarot)
 658 
 659      ! Update internal values
 660      call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
 661 
 662      call timab(553,2,tsec)
 663    else
 664      if (pawtab(1)%has_kij  ==1) pawtab(1:psps%ntypat)%has_kij  =2
 665      if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2
 666    end if
 667    psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
 668    call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot)
 669    pawtab(:)%usepawu=0
 670    pawtab(:)%useexexch=0
 671    pawtab(:)%exchmix=zero
 672 !  if (dtset%usepawu>0.or.dtset%useexexch>0) then
 673    call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
 674 &   dtset%jpawu,dtset%lexexch,dtset%lpawu,ntypat,pawang,dtset%pawprtvol,pawrad,&
 675 &   pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
 676 !  end if
 677    compch_fft=-1.d5;compch_sph=-1.d5
 678    usexcnhat=maxval(pawtab(:)%usexcnhat)
 679    usecprj=dtset%pawusecp
 680 !  2-Check overlap
 681    call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred)
 682 !  3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r)
 683    ABI_DATATYPE_ALLOCATE(pawfgrtab,(my_natom))
 684    if (my_natom>0) then
 685      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
 686 &     mpi_atmtab=mpi_enreg%my_atmtab)
 687      call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,&
 688 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 689      ABI_DEALLOCATE(l_size_atm)
 690    end if
 691    use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0)
 692    optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
 693    if (use_nhat_gga) then
 694      optgr1=dtset%pawstgylm
 695      if (rfphon==1) optgr2=1
 696    end if
 697    if (rfphon==1.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1) then ! LB2016-11-28 : Why not rfelfd==1?
 698      if (optgr1==0) optgr1=dtset%pawstgylm
 699      if (optgr2==0) optgr2=dtset%pawstgylm
 700      if (optrad==0.and.(.not.qeq0.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1)) optrad=1  ! LB2016-11-28 : Why not rfelfd==1?
 701    end if
 702    if (rfelfd==1.or.rfelfd==3.or.rf2_dkde==1) then
 703      if (optgr1==0) optgr1=dtset%pawstgylm
 704    end if
 705    call status(0,dtfil%filstat,iexit,level,'call nhatgrid ')
 706    call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,&
 707 &   optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
 708 &   comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab )
 709 !  Compute exp(iq.r) factors around the atoms
 710    if (.not.qeq0) then
 711      do iatom=1,my_natom
 712        iatom_tot=iatom; if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom)
 713        if (allocated(pawfgrtab(iatom)%expiqr)) then
 714          ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
 715        end if
 716        ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd))
 717        call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,dtset%qptn,&
 718 &       pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
 719        pawfgrtab(iatom)%expiqr_allocated=1
 720      end do
 721    end if
 722 !  4-Allocate PAW specific arrays
 723    ABI_DATATYPE_ALLOCATE(paw_an,(my_natom))
 724    ABI_DATATYPE_ALLOCATE(paw_ij,(my_natom))
 725    call paw_an_nullify(paw_an)
 726    call paw_ij_nullify(paw_ij)
 727    has_kxc=0;nkxc1=0;cplex=1
 728    has_dijnd=0; if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
 729    has_diju=0; if(dtset%usepawu==5.or.dtset%usepaw==6) has_diju=1
 730    if (rfphon/=0.or.rfelfd==1.or.rfelfd==3.or.rfstrs/=0.or.rf2_dkde/=0) then
 731      has_kxc=1
 732      call pawxc_get_nkxc(nkxc1,dtset%nspden,dtset%xclevel)
 733    end if
 734    call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,0,dtset%nspden,&
 735 &   cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,has_vxc_ex=1,has_kxc=has_kxc,&
 736 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 737    call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,&
 738 &   natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,&
 739 &   has_dijso=1,has_dijU=has_diju,has_pawu_occ=1,has_exexch_pot=1,nucdipmom=dtset%nucdipmom,&
 740 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 741 
 742  else ! PAW vs NCPP
 743    usexcnhat=0;usecprj=0
 744    use_nhat_gga=.false.
 745    ABI_DATATYPE_ALLOCATE(paw_an,(0))
 746    ABI_DATATYPE_ALLOCATE(paw_ij,(0))
 747    ABI_DATATYPE_ALLOCATE(pawfgrtab,(0))
 748  end if
 749 
 750  ABI_ALLOCATE(rhog,(2,nfftf))
 751  ABI_ALLOCATE(rhor,(nfftf,dtset%nspden))
 752 
 753 !Read ground-state charge density from diskfile in case getden /= 0
 754 !or compute it from wfs that were read previously : rhor as well as rhog
 755 
 756  if (dtset%getden /= 0 .or. dtset%irdden /= 0) then
 757    ! Read rho1(r) from a disk file and broadcast data.
 758    ! This part is not compatible with MPI-FFT (note single_proc=.True. below)
 759 
 760    rdwr=1;rdwrpaw=psps%usepaw;if(ireadwf0/=0) rdwrpaw=0
 761    if (rdwrpaw/=0) then
 762      ABI_DATATYPE_ALLOCATE(pawrhoij_read,(natom))
 763      call pawrhoij_nullify(pawrhoij_read)
 764      nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
 765      call pawrhoij_alloc(pawrhoij_read,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,&
 766 &     dtset%nsppol,dtset%typat,pawtab=pawtab)
 767    else
 768      ABI_DATATYPE_ALLOCATE(pawrhoij_read,(0))
 769    end if
 770 
 771 !    MT july 2013: Should we read rhoij from the density file ?
 772    call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, &
 773    hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr)
 774    etotal = hdr_den%etot; call hdr_free(hdr_den)
 775 
 776    if (rdwrpaw/=0) then
 777      call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld)
 778      call pawrhoij_free(pawrhoij_read)
 779    end if
 780    ABI_DATATYPE_DEALLOCATE(pawrhoij_read)
 781 
 782 !  Compute up+down rho(G) by fft
 783    ABI_ALLOCATE(work,(nfftf))
 784    work(:)=rhor(:,1)
 785    call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 786    ABI_DEALLOCATE(work)
 787 
 788  else
 789    izero=0
 790 !  Obtain the charge density from read wfs
 791 !  Be careful: in PAW, compensation density has to be added !
 792    tim_mkrho=4
 793    paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented
 794    paw_dmft%use_dmft=0 ! respfn with dmft not implemented
 795    if (psps%usepaw==1) then
 796      ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
 797      ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
 798      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
 799 &     mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
 800      call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
 801      ABI_DEALLOCATE(rhowfg)
 802      ABI_DEALLOCATE(rhowfr)
 803    else
 804      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
 805 &     mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
 806    end if
 807  end if ! getden
 808 
 809 !In PAW, compensation density has eventually to be added
 810  nhatgrdim=0;nhatdim=0
 811  ABI_ALLOCATE(nhatgr,(0,0,0))
 812  if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then
 813    nhatdim=1
 814    ABI_ALLOCATE(nhat,(nfftf,dtset%nspden))
 815    call timab(558,1,tsec)
 816    nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat
 817    ider=2*nhatgrdim
 818    if (nhatgrdim>0)  then
 819      ABI_DEALLOCATE(nhatgr)
 820      ABI_ALLOCATE(nhatgr,(nfftf,dtset%nspden,3))
 821    end if
 822    izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero
 823    call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,&
 824 &   nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,&
 825 &   nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, &
 826 &   mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom)
 827    if (dtset%getden==0) then
 828      rhor(:,:)=rhor(:,:)+nhat(:,:)
 829      call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 830    end if
 831    call timab(558,2,tsec)
 832  else
 833    ABI_ALLOCATE(nhat,(0,0))
 834  end if
 835 
 836 !The GS irrzon and phnons were only needed to symmetrize the GS density
 837  ABI_DEALLOCATE(irrzon)
 838  ABI_DEALLOCATE(phnons)
 839 
 840 !jmb 2012 write(std_out,'(a)')' ' ! needed to make ibm6_xlf12 pass tests. No idea why this works. JWZ 5 Sept 2011
 841 !Will compute now the total potential
 842 
 843 !Compute local ionic pseudopotential vpsp and core electron density xccc3d:
 844  n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
 845  ABI_ALLOCATE(xccc3d,(n3xccc))
 846  ABI_ALLOCATE(vpsp,(nfftf))
 847 
 848 !Determine by which method the local ionic potential and/or
 849 ! the pseudo core charge density have to be computed
 850 !Local ionic potential:
 851 ! Method 1: PAW ; Method 2: Norm-conserving PP
 852  vloc_method=1;if (psps%usepaw==0) vloc_method=2
 853 !Pseudo core charge density:
 854 ! Method 1: PAW, nc_xccc_gspace ; Method 2: Norm-conserving PP
 855  coredens_method=1;if (psps%usepaw==0) coredens_method=2
 856  if (psps%nc_xccc_gspace==1) coredens_method=1
 857  if (psps%nc_xccc_gspace==0) coredens_method=2
 858 
 859 !Local ionic potential and/or pseudo core charge by method 1
 860  if (vloc_method==1.or.coredens_method==1) then
 861    call timab(562,1,tsec)
 862    optv=0;if (vloc_method==1) optv=1
 863    optn=0;if (coredens_method==1) optn=n3xccc/nfftf
 864    optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
 865    call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,&
 866 &   dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,&
 867 &   ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,&
 868 &   dtset%qprtrb,dum_rhog,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl)
 869    call timab(562,2,tsec)
 870  end if
 871 
 872 !Local ionic potential by method 2
 873  if (vloc_method==2) then
 874    option=1
 875    ABI_ALLOCATE(dyfrlo_indx,(3,3,natom))
 876    ABI_ALLOCATE(grtn_indx,(3,natom))
 877    call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,&
 878 &   grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,&
 879 &   nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,&
 880 &   dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
 881    ABI_DEALLOCATE(dyfrlo_indx)
 882    ABI_DEALLOCATE(grtn_indx)
 883  end if
 884 
 885 !Pseudo core electron density by method 2
 886  if (coredens_method==2.and.psps%n1xccc/=0) then
 887    option=1
 888    ABI_ALLOCATE(dyfrx2,(3,3,natom))
 889    ABI_ALLOCATE(vxc,(0,0)) ! dummy
 890    call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,&
 891 &   ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,&
 892 &   psps%xcccrc,psps%xccc1d,xccc3d,xred)
 893    ABI_DEALLOCATE(dyfrx2)
 894    ABI_DEALLOCATE(vxc) ! dummy
 895  end if
 896 
 897 !Set up hartree and xc potential. Compute kxc here.
 898  ABI_ALLOCATE(vhartr,(nfftf))
 899 
 900  call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr)
 901 
 902  option=2 ; nk3xc=1
 903  nkxc=2*min(dtset%nspden,2)-1;if(dtset%xclevel==2)nkxc=12*min(dtset%nspden,2)-5
 904  call check_kxc(dtset%ixc,dtset%optdriver)
 905  ABI_ALLOCATE(kxc,(nfftf,nkxc))
 906  ABI_ALLOCATE(vxc,(nfftf,dtset%nspden))
 907 
 908  _IBM6("Before rhotoxc")
 909 
 910  call xcdata_init(xcdata,dtset=dtset)
 911  call rhotoxc(enxc,kxc,mpi_enreg,nfftf,ngfftf,&
 912 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,&
 913 & rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartr)
 914 
 915 !Compute local + Hxc potential, and subtract mean potential.
 916  ABI_ALLOCATE(vtrial,(nfftf,dtset%nspden))
 917  do ispden=1,min(dtset%nspden,2)
 918    do ifft=1,nfftf
 919      vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
 920    end do
 921  end do
 922  if (dtset%nspden==4) then
 923    do ispden=3,4
 924      do ifft=1,nfftf
 925        vtrial(ifft,ispden)=vxc(ifft,ispden)
 926      end do
 927    end do
 928  end if
 929  ABI_DEALLOCATE(vhartr)
 930 
 931 
 932  if(dtset%prtvol==-level)then
 933    call wrtout(std_out,' respfn: ground-state density and potential set up.','COLL')
 934  end if
 935 
 936 !PAW: compute Dij quantities (psp strengths)
 937  if (psps%usepaw==1)then
 938    cplex=1;ipert=0;option=1
 939    nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1
 940 
 941    call pawdenpot(compch_sph,epaw,epawdc,ipert,dtset%ixc,my_natom,natom,dtset%nspden,&
 942 &   ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,&
 943 &   pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,&
 944 &   dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, &
 945 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 946 
 947    call timab(561,1,tsec)
 948    call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,&
 949 &   dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,&
 950 &   pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,&
 951 &   dtset%spnorbscl,ucvol,dtset%charge,vtrial,vxc,xred,nucdipmom=dtset%nucdipmom,&
 952 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 953    call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,&
 954 &   paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
 955 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 956    call timab(561,2,tsec)
 957 
 958  end if
 959 
 960 !-----2. Frozen-wavefunctions and Ewald(q=0) parts of 2DTE
 961 
 962  dyfr_nondiag=0;if (psps%usepaw==1.and.rfphon==1) dyfr_nondiag=1
 963  dyfr_cplex=1;if (psps%usepaw==1.and.rfphon==1.and.(.not.qeq0)) dyfr_cplex=2
 964  ABI_ALLOCATE(dyew,(2,3,natom,3,natom))
 965  ABI_ALLOCATE(dyewq0,(3,3,natom))
 966  ABI_ALLOCATE(dyfrlo,(3,3,natom))
 967  ABI_ALLOCATE(dyfrx2,(3,3,natom))
 968  ABI_ALLOCATE(dyfrnl,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag))
 969  ABI_ALLOCATE(dyfrwf,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag))
 970  ABI_ALLOCATE(becfrnl,(3,natom,3*pawbec))
 971  ABI_ALLOCATE(piezofrnl,(6,3*pawpiezo))
 972  ABI_ALLOCATE(dyvdw,(2,3,natom,3,natom*usevdw))
 973  dyew(:,:,:,:,:)=zero
 974  dyewq0(:,:,:)=zero
 975  dyfrnl(:,:,:,:,:)=zero
 976  dyfrwf(:,:,:,:,:)=zero
 977  dyfrlo(:,:,:)=zero
 978  dyfrx2(:,:,:)=zero
 979  if (usevdw==1) dyvdw(:,:,:,:,:)=zero
 980  if (pawbec==1) becfrnl(:,:,:)=zero
 981  if (pawpiezo==1) piezofrnl(:,:)=zero
 982 
 983  ABI_ALLOCATE(eltcore,(6,6))
 984  ABI_ALLOCATE(elteew,(6+3*natom,6))
 985  ABI_ALLOCATE(eltfrhar,(6,6))
 986  ABI_ALLOCATE(eltfrnl,(6+3*natom,6))
 987  ABI_ALLOCATE(eltfrloc,(6+3*natom,6))
 988  ABI_ALLOCATE(eltfrkin,(6,6))
 989  ABI_ALLOCATE(eltfrxc,(6+3*natom,6))
 990  ABI_ALLOCATE(eltvdw,(6+3*natom,6*usevdw))
 991  eltcore(:,:)=zero
 992  elteew(:,:)=zero
 993  eltfrnl(:,:)=zero
 994  eltfrloc(:,:)=zero
 995  eltfrkin(:,:)=zero
 996  eltfrhar(:,:)=zero
 997  eltfrxc(:,:)=zero
 998  if (usevdw==1) eltvdw(:,:)=zero
 999 
1000 !Section common to all perturbations
1001 !Compute the nonlocal part of the elastic tensor and/or dynamical matrix
1002  if (rfstrs/=0.or.rfphon==1.or.dtset%efmas>0.or.pawbec==1.or.pawpiezo==1)then
1003    call d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,&
1004 &   dyfr_nondiag,efmasdeg,efmasval,eigen0,eltfrnl,gsqcut,has_allddk,indsym,kg,mgfftf,&
1005 &   mpi_enreg,psps%mpsang,my_natom,natom,nfftf,ngfft,ngfftf,&
1006 &   npwarr,occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,&
1007 &   pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,rprimd,rfphon,&
1008 &   rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr)
1009  end if
1010 
1011 !No more need of these local derivatives
1012  if (rfphon==1.and.psps%usepaw==1.and.(.not.use_nhat_gga)) then
1013    do iatom=1,my_natom
1014      if (allocated(pawfgrtab(iatom)%gylmgr2)) then
1015        ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
1016      end if
1017      pawfgrtab(iatom)%gylmgr2_allocated=0
1018    end do
1019  end if
1020 
1021 !Section for the atomic displacement/electric field perturbations
1022  if (rfphon==1) then
1023 
1024 !  Compute the local of the dynamical matrix
1025 !  dyfrnl has not yet been symmetrized, but will be in the next routine
1026    call dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrx2,dyfr_cplex,dyfr_nondiag,&
1027 &   gmet,gprimd,gsqcut,indsym,mgfftf,mpi_enreg,psps%mqgrid_vl,&
1028 &   natom,nattyp, nfftf,ngfftf,dtset%nspden,dtset%nsym,ntypat,&
1029 &   psps%n1xccc,n3xccc,dtset%paral_kgb,psps,pawtab,ph1df,psps%qgrid_vl,&
1030 &   dtset%qptn,rhog,rprimd,symq,symrec,dtset%typat,ucvol,&
1031 &   psps%usepaw,psps%vlspl,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred)
1032 
1033    _IBM6("Before dfpt_ewald")
1034 
1035 !  Compute Ewald (q=0) contribution
1036    sumg0=0;qphon(:)=zero
1037    call status(0,dtfil%filstat,iexit,level,'call dfpt_ewald(1)')
1038    call dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,psps%ziontypat,&
1039 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1040    option=1
1041    call q0dy3_calc(natom,dyewq0,dyew,option)
1042 !  Calculate the DFT-D2 vdw part of the dynamical matrix
1043    if (usevdw==1.and.dtset%vdw_xc==5) then
1044      call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,dtset%vdw_tol,&
1045 &     xred,psps%znuclpsp,dyn_vdw_dftd2=dyvdw,qphon=dtset%qptn)
1046    end if
1047 ! Calculate the DFT-D3/D3(BJ) vdw part of the dynamical matrix
1048    if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then
1049      call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,&
1050 &     dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,&
1051 &     dyn_vdw_dftd3=dyvdw,qphon=dtset%qptn)
1052    end if
1053 !The frozen-wavefunction part of the dynamical matrix is now:
1054 !  d2frnl:  non-local contribution
1055 !  dyfrlo:  local contribution
1056 !  dyfrx2:  2nd order xc core correction contribution
1057 !  dyew  : Ewald contribution
1058 !  dyvdw : vdw DFT-D contribution
1059 !  dyfrwf:  all contributions
1060 !  In case of PAW, it misses a term coming from the perturbed overlap operator
1061  end if
1062 
1063 !Section for the strain perturbation
1064  if(rfstrs/=0) then
1065 
1066 !  Verify that k-point set has full space-group symmetry; otherwise exit
1067    call status(0,dtfil%filstat,iexit,level,'call symkchk ')
1068    timrev=1
1069    if (symkchk(dtset%kptns,dtset%nkpt,dtset%nsym,symrec,timrev,message) /= 0) then
1070      MSG_ERROR(message)
1071    end if
1072 
1073 !  Calculate the kinetic part of the elastic tensor
1074    call dfpt_eltfrkin(cg,eltfrkin,dtset%ecut,dtset%ecutsm,dtset%effmass_free,&
1075 &   dtset%istwfk,kg,dtset%kptns,dtset%mband,dtset%mgfft,dtset%mkmem,mpi_enreg,&
1076 &   dtset%mpw,dtset%nband,dtset%nkpt,ngfft,npwarr,&
1077 &   dtset%nspinor,dtset%nsppol,occ,rprimd,dtset%wtk)
1078 
1079 !  Calculate the hartree part of the elastic tensor
1080    call dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfftf,ngfftf,rhog)
1081 
1082 !  Calculate the xc part of the elastic tensor
1083    call dfpt_eltfrxc(atindx,dtset,eltfrxc,enxc,gsqcut,kxc,mpi_enreg,mgfftf,&
1084 &   nattyp,nfftf,ngfftf,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1df,psps,rhor,rprimd,&
1085 &   usexcnhat,vxc,xccc3d,xred)
1086 
1087 !  Calculate the local potential part of the elastic tensor
1088    call dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfftf,mpi_enreg,psps%mqgrid_vl,&
1089 &   natom,nattyp,nfftf,ngfftf,ntypat,ph1df,psps%qgrid_vl,rhog,psps%vlspl)
1090 
1091 !  Calculate the Ewald part of the elastic tensor
1092    call elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,&
1093 &   dtset%typat,ucvol,xred,psps%ziontypat,&
1094 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1095 
1096 !  Calculate the DFT-D2 vdw part of the elastic tensor
1097    if (usevdw==1.and.dtset%vdw_xc==5) then
1098      option=1; if (rfphon==1) option=0
1099      call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,dtset%vdw_tol,&
1100 &     xred,psps%znuclpsp,elt_vdw_dftd2=eltvdw)
1101    end if
1102 !  Calculate the DFT-D3/D3(BJ) vdw part of the elastic tensor
1103    if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then
1104      option=1; if (rfphon==1) option=0
1105      call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,&
1106 &     dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,&
1107 &     elt_vdw_dftd3=eltvdw)
1108    end if
1109 !  Calculate the psp core energy part of elastic tensor (trivial)
1110    eltcore(1:3,1:3)=ecore/ucvol
1111 
1112 !The frozen-wavefunction part of the elastic tensor is now:
1113 !  eltfrnl:  non-local contribution
1114 !  eltfrloc: local contribution
1115 !  eltfrkin: kinetic contribution
1116 !  eltfrhar: Hartree contribution
1117 !  eltfrx:   XC contribution
1118 !  eltcore:  psps core contribution
1119 !  elteew:   Ewald contribution
1120 !  eltvdw:   vdw DFT-D contribution
1121 !  In case of PAW, it misses a term coming from the perturbed overlap operator
1122  end if
1123 
1124  ABI_DEALLOCATE(vpsp)
1125  ABI_DEALLOCATE(xccc3d)
1126 
1127  if(dtset%prtvol==-level)then
1128    call wrtout(std_out,' respfn: frozen wavef. and Ewald(q=0) part of 2DTE done.','COLL')
1129  end if
1130 
1131  call timab(136,2,tsec)
1132 
1133 !-----3. Initialisation of 1st response, taking into account the q vector.
1134 
1135  call timab(137,1,tsec)
1136 
1137  write(message,'(3a)')ch10,' ==>  initialize data related to q vector <== ',ch10
1138  call wrtout(std_out,message,'COLL')
1139  call wrtout(ab_out,message,'COLL')
1140 
1141  qphon(:)=dtset%qptn(:)
1142  sumg0=1
1143 
1144 !Treat the case of q=0 or q too close to 0
1145  qzero=0
1146  if(qeq0)then
1147    qphon(:)=zero
1148    write(message,'(3a)')&
1149 &   ' respfn : the norm of the phonon wavelength (as input) was small (<1.d-7).',ch10,&
1150 &   '  q has been set exactly to (0 0 0)'
1151    call wrtout(std_out,message,'COLL')
1152    sumg0=0
1153    qzero=1
1154  else
1155    if(rfelfd/=0 .or. rfstrs/=0 .or. rfddk /= 0  .or. rf2_dkdk /= 0 .or. rf2_dkde /= 0) then
1156 !    Temporarily, ...
1157      write(message, '(a,a,a,3es16.6,a,a,5(a,i2),a,a,a)' )ch10,&
1158 &     'The treatment of non-zero wavevector q is restricted to phonons.',&
1159 &     'However, the input normalized qpt is',qphon(:),',',ch10,&
1160 &     'while rfelfd=',rfelfd,', rfddk=',rfddk,', rf2_dkdk=',rf2_dkdk,', rf2_dkde=',rf2_dkde,&
1161 &     ' and rfstrs=',rfstrs,'.',ch10,&
1162 &     'Action: change qpt, or rfelfd, or rfstrs in the input file.'
1163      MSG_ERROR(message)
1164    else if(rfasr.eq.2)then
1165      write(message,'(2a)')ch10,' rfasr=2 not allowed with q/=0 => rfasr was reset to 0.'
1166      MSG_WARNING(message)
1167      rfasr=0
1168    end if
1169  end if
1170 
1171  _IBM6("Before irreducible_set_pert")
1172 
1173 !Determine the symmetrical perturbations
1174  ABI_ALLOCATE(pertsy,(3,natom+6))
1175  call irreducible_set_pert(indsym,natom+6,natom,dtset%nsym,pertsy,rfdir,rfpert,symq,symrec,dtset%symrel)
1176  write(message,'(a)') ' The list of irreducible perturbations for this q vector is:'
1177  call wrtout(ab_out,message,'COLL')
1178  call wrtout(std_out,message,'COLL')
1179  ii=1
1180  do ipert=1,natom+6
1181    do idir=1,3
1182      if(rfpert(ipert)==1.and.rfdir(idir)==1)then
1183        if( pertsy(idir,ipert)==1 )then
1184          write(message, '(i5,a,i2,a,i4)' )ii,')    idir=',idir,'    ipert=',ipert
1185          call wrtout(ab_out,message,'COLL')
1186          call wrtout(std_out,message,'COLL')
1187          ii=ii+1
1188        end if
1189      end if
1190    end do
1191  end do
1192 
1193 !test if the user left default rfdir 0 0 0
1194  if (ii==1 .and. rf2_dkdk==0 .and. rf2_dkde==0) then
1195    write(message,'(5a)')ch10,&
1196 &   ' WARNING: no perturbations to be done at this q-point.',ch10,&
1197 &   ' You may have forgotten to set the rfdir or rfatpol variables. Continuing normally.',ch10
1198    call wrtout(ab_out,message,'COLL')
1199    MSG_WARNING(message)
1200  end if
1201 
1202  if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then
1203    ABI_ALLOCATE(rfpert_nl,(3,natom+2,3,natom+2,3,natom+2))
1204    rfpert_nl = 0
1205    rfpert_nl(:,natom+2,:,natom+2,:,natom+2) = 1
1206    rfpert_nl(:,1:natom,:,natom+2,:,natom+2) = 1
1207    rfpert_nl(:,natom+2,:,1:natom,:,natom+2) = 1
1208    rfpert_nl(:,natom+2,:,natom+2,:,1:natom) = 1
1209    call sytens(indsym,natom+2,natom,dtset%nsym,rfpert_nl,symrec,dtset%symrel)
1210    write(message, '(a,a,a,a,a)' ) ch10, &
1211 &   ' The list of irreducible elements of the Raman and non-linear',&
1212 &   ch10,' optical susceptibility tensors is:',ch10
1213    call wrtout(std_out,message,'COLL')
1214 
1215    write(message,'(12x,a)')&
1216 &   'i1pert  i1dir   i2pert  i2dir   i3pert  i3dir'
1217    call wrtout(std_out,message,'COLL')
1218    n1 = 0
1219    rf2_dirs_from_rfpert_nl(:,:) = 0
1220    do i1pert = 1, natom + 2
1221      do i1dir = 1, 3
1222        do i2pert = 1, natom + 2
1223          do i2dir = 1, 3
1224            do i3pert = 1, natom + 2
1225              do i3dir = 1,3
1226                if (rfpert_nl(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
1227                  n1 = n1 + 1
1228                  write(message,'(2x,i4,a,6(5x,i3))') n1,')', &
1229 &                 i1pert,i1dir,i2pert,i2dir,i3pert,i3dir
1230                  call wrtout(std_out,message,'COLL')
1231                  if (i2pert==natom+2) then
1232                    if (i3pert==natom+2) then
1233                      rf2_dirs_from_rfpert_nl(i3dir,i2dir) = 1
1234                    else if (i1pert==natom+2) then
1235                      rf2_dirs_from_rfpert_nl(i1dir,i2dir) = 1
1236                    end if
1237                  end if
1238                end if
1239              end do
1240            end do
1241          end do
1242        end do
1243      end do
1244    end do
1245    write(message,'(a,a)') ch10,ch10
1246    call wrtout(std_out,message,'COLL')
1247 
1248    write(message,'(a)') 'rf2_dirs_from_rfpert_nl :'
1249    call wrtout(std_out,message,'COLL')
1250    do i1dir = 1, 3
1251      do i2dir = 1, 3
1252        write(message,'(3(a,i1))') ' ',i1dir,' ',i2dir,' : ',rf2_dirs_from_rfpert_nl(i1dir,i2dir)
1253        call wrtout(std_out,message,'COLL')
1254      end do
1255    end do
1256  end if
1257 
1258 !Contribution to the dynamical matrix from ion-ion energy
1259  if(rfphon==1)then
1260    call status(0,dtfil%filstat,iexit,level,'call dfpt_ewald(2)')
1261    call dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,psps%ziontypat, &
1262 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1263    call q0dy3_apply(natom,dyewq0,dyew)
1264  end if
1265 
1266 !1-order contribution of the xc core correction to the dynamical matrix
1267  ABI_ALLOCATE(dyfrx1,(2,3,natom,3,natom))
1268  dyfrx1(:,:,:,:,:)=zero
1269  if(rfphon==1.and.psps%n1xccc/=0)then
1270    ABI_ALLOCATE(blkflgfrx1,(3,natom,3,natom))
1271 !FR non-collinear magnetism
1272    if (dtset%nspden==4) then
1273      call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,&
1274 &     psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,dtset%nspden,&
1275 &     ntypat,psps%n1xccc,dtset%paral_kgb,psps,pawtab,ph1df,psps%qgrid_vl,qphon,&
1276 &     rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred,rhor=rhor,vxc=vxc)
1277    else
1278      call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,&
1279 &     psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,dtset%nspden,&
1280 &     ntypat,psps%n1xccc,dtset%paral_kgb,psps,pawtab,ph1df,psps%qgrid_vl,qphon,&
1281 &     rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred)
1282    end if
1283  end if
1284 
1285 !Deallocate the arrays that were needed only for the frozen wavefunction part
1286  ABI_DEALLOCATE(ph1d)
1287  ABI_DEALLOCATE(ph1df)
1288  ABI_DEALLOCATE(cg)
1289  ABI_DEALLOCATE(kg)
1290  ABI_DEALLOCATE(npwarr)
1291  if(xmpi_paral==1) then
1292    ABI_DEALLOCATE(mpi_enreg%proc_distrb)
1293  end if
1294 
1295  ABI_ALLOCATE(blkflg,(3,mpert,3,mpert))
1296  ABI_ALLOCATE(d2eig0,(2,3,mpert,3,mpert))
1297  ABI_ALLOCATE(d2k0,(2,3,mpert,3,mpert))
1298  ABI_ALLOCATE(d2lo,(2,3,mpert,3,mpert))
1299  ABI_ALLOCATE(d2loc0,(2,3,mpert,3,mpert))
1300  ABI_ALLOCATE(d2nfr,(2,3,mpert,3,mpert))
1301  ABI_ALLOCATE(d2nl,(2,3,mpert,3,mpert))
1302  ABI_ALLOCATE(d2nl0,(2,3,mpert,3,mpert))
1303  ABI_ALLOCATE(d2nl1,(2,3,mpert,3,mpert))
1304  ABI_ALLOCATE(d2vn,(2,3,mpert,3,mpert))
1305  ABI_ALLOCATE(d2ovl,(2,3,mpert,3,mpert*psps%usepaw))
1306  blkflg(:,:,:,:)=0
1307  d2eig0(:,:,:,:,:)=zero ; d2k0(:,:,:,:,:)=zero
1308  d2lo(:,:,:,:,:)=zero   ; d2loc0(:,:,:,:,:)=zero
1309  d2nfr(:,:,:,:,:)=zero  ; d2nl(:,:,:,:,:)=zero
1310  d2nl0(:,:,:,:,:)=zero  ; d2nl1(:,:,:,:,:)=zero
1311  d2vn(:,:,:,:,:)=zero
1312  if (psps%usepaw==1) d2ovl(:,:,:,:,:)=zero
1313 
1314  prtbbb=dtset%prtbbb
1315  ABI_ALLOCATE(d2bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb))
1316  ABI_ALLOCATE(d2cart_bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb))
1317  if(prtbbb==1)then
1318    d2cart_bbb(:,:,:,:,:,:)=zero
1319    d2bbb(:,:,:,:,:,:)=zero
1320  end if
1321 
1322  dim_eig2nkq = 0
1323  if(dtset%ieig2rf /= 0) dim_eig2nkq = 1
1324  ABI_ALLOCATE(eig2nkq,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq))
1325  dim_eigbrd=0
1326  if(dtset%ieig2rf /= 0 .and. dtset%smdelta>0 ) dim_eigbrd = 1
1327  ABI_ALLOCATE(eigbrd,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eigbrd))
1328 
1329  call timab(137,2,tsec)
1330 
1331 
1332 !Check whether exiting was required by the user.
1333 !If found then do not start minimization steps
1334 !At this first call to exit_check, initialize cpus
1335  cpus=dtset%cpus
1336  if(abs(cpus)>1.0d-5)cpus=cpus+cpui
1337  openexit=1; if(dtset%chkexit==0) openexit=0
1338  call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
1339 
1340 !TEMPORARY: for testing purpose only
1341 ! if (rfstrs/=0.and.dtset%usepaw==1) iexit=1
1342 
1343  _IBM6("Before dfpt_looppert")
1344 
1345  if (iexit==0) then
1346 !  #######################################################################
1347    write(message,'(a,80a)')ch10,('=',mu=1,80)
1348    call wrtout(ab_out,message,'COLL')
1349    call wrtout(std_out,message,'COLL')
1350 
1351    ddkfil(:)=0
1352 
1353    ! MGNAG: WHY THIS?  all these variables are declared as optional pointers in dfpt_looppert!
1354    ! but they are allocated here so why using pointers! Moreover OPTIONAL arguments MUST
1355    ! be passed by keyword for better clarity and robustness!
1356    ! People should learn how to program Fort90 before being allowed to change the code!
1357    ! v5[26] crashes in dfpt_looppert
1358    ! The best solution would be using a datatype to gather the results!
1359    ABI_ALLOCATE(clflg,(3,mpert))
1360    if(dtset%ieig2rf > 0) then
1361      ABI_ALLOCATE(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1362      ABI_ALLOCATE(eigenq_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1363      ABI_ALLOCATE(occ_rbz_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1364    end if
1365    if(dtset%efmas > 0) then
1366      ABI_ALLOCATE(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol))
1367    end if
1368 !  Note that kg, cg, eigen0, mpw and npwarr are NOT passed to dfpt_looppert :
1369 !  they will be reinitialized for each perturbation, with an eventually
1370 !  reduced set of k point, thanks to the use of symmetry operations.
1371    call dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,&
1372 &   ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,dyfr_cplex,dyfr_nondiag,&
1373 &   d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,&
1374 &   eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
1375 &   etotal,fermie,iexit,indsym,kxc,&
1376 &   dtset%mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,&
1377 &   nfftf,nhat,dtset%nkpt,nkxc,dtset%nspden,dtset%nsym,occ,&
1378 &   paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
1379 &   pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,&
1380 &   usecprj,usevdw,vtrial,vxc,vxcavg,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,&
1381 &   eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0)
1382 
1383 !  #####################################################################
1384  end if
1385 
1386  call timab(138,1,tsec)
1387 
1388  write(message, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,&
1389 & ' ---- first-order wavefunction calculations are completed ----',ch10
1390  call wrtout(ab_out,message,'COLL')
1391  call wrtout(std_out,message,'COLL')
1392 
1393  ABI_DEALLOCATE(vxc)
1394 
1395  if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then
1396    ABI_DEALLOCATE(rfpert_nl)
1397  end if
1398 
1399 !Output of the localization tensor
1400  if ( rfpert(natom+1) /= 0 .and. (me == 0) .and. dtset%occopt<=2) then
1401    call wrtloctens(blkflg,d2bbb,d2nl,dtset%mband,mpert,natom,dtset%prtbbb,rprimd,psps%usepaw)
1402  end if
1403 
1404 !The perturbation  natom+1 was only an auxiliary perturbation,
1405 !needed to construct the electric field response, so its flag is now set to 0.
1406 !rfpert(natom+1)=0
1407 
1408 !Were 2DTE computed ?
1409  if(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfuser==0 .and. rfmagn==0)then
1410 
1411    write(message,'(a,a)' )ch10,' respfn : d/dk was computed, but no 2DTE, so no DDB output.'
1412    call wrtout(std_out,message,'COLL')
1413    call wrtout(ab_out,message,'COLL')
1414 
1415 !  If 2DTE were computed, only one processor must output them and compute
1416 !  frequencies.
1417  else if(me==0)then
1418 
1419    write(message,'(a,a)' )ch10,' ==> Compute Derivative Database <== '
1420    call wrtout(std_out,message,'COLL')
1421    call wrtout(ab_out,message,'COLL')
1422 
1423 !  In the RESPFN code, dfpt_nstdy and stady3 were called here
1424    d2nfr(:,:,:,:,:)=d2lo(:,:,:,:,:)+d2nl(:,:,:,:,:)
1425    if (psps%usepaw==1) d2nfr(:,:,:,:,:)=d2nfr(:,:,:,:,:)+d2ovl(:,:,:,:,:)
1426 
1427    zero_by_symm=1
1428    if(dtset%rfmeth<0)zero_by_symm=0
1429 
1430 !  In case of bbb decomposition
1431    if(prtbbb==1)then
1432      ABI_ALLOCATE(blkflg1,(3,mpert,3,mpert))
1433      ABI_ALLOCATE(blkflg2,(3,mpert,3,mpert))
1434      blkflg2(:,:,:,:) = blkflg(:,:,:,:)
1435      do ipert = 1, mpert
1436        do ipert2 = 1, mpert
1437          if ((ipert /= natom + 2).and.(ipert>natom).and.(ipert2/=natom+2)) then
1438            blkflg2(:,ipert2,:,ipert) = 0
1439          end if
1440        end do
1441      end do
1442      ABI_ALLOCATE(d2tmp,(2,3,mpert,3,mpert))
1443      do iband = 1,dtset%mband
1444        d2tmp(:,:,:,:,:)=zero
1445        blkflg1(:,:,:,:) = blkflg2(:,:,:,:)
1446        d2tmp(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband)
1447        call d2sym3(blkflg1,d2tmp,indsym,mpert,natom,dtset%nsym,qphon,symq,&
1448 &       symrec,dtset%symrel,timrev,zero_by_symm)
1449        d2bbb(:,:,:,:,iband,iband) = d2tmp(:,:,natom+2,:,:)
1450      end do
1451      ABI_DEALLOCATE(blkflg1)
1452      ABI_DEALLOCATE(blkflg2)
1453      ABI_DEALLOCATE(d2tmp)
1454    end if
1455 
1456 !  Complete the d2nfr matrix by symmetrization of the existing elements
1457    !write(std_out,*)"blkflg before d2sym3: ", blkflg
1458    call d2sym3(blkflg,d2nfr,indsym,mpert,natom,dtset%nsym,qphon,symq,symrec,dtset%symrel,timrev,zero_by_symm)
1459    !write(std_out,*)"blkflg after d2sym3: ", blkflg
1460 
1461    if(rfphon==1.and.psps%n1xccc/=0)then
1462 !    Complete the dyfrx1 matrix by symmetrization of the existing elements
1463      call d2sym3(blkflgfrx1,dyfrx1,indsym,natom,natom,dtset%nsym,qphon,symq,symrec,dtset%symrel,timrev,zero_by_symm)
1464    end if
1465 
1466 !  Note that d2sym3 usually complete the 2nd-order matrix
1467 !  with elements that are zero by symmetry, automatically,
1468 !  unless it has been explicitly asked not to do so.
1469 !  blkflg is then set to 1 for these matrix elements, even if there has be no calculation.
1470 
1471 !  Add the frozen-wf (dyfrwf) part to the ewald part (dyew),
1472 !  the part 1 of the frozen wf part of the xc core correction
1473 !  (dyfrx1) and the non-frozen part (dynfr) to get the second-order
1474 !  derivative matrix (d2matr), then
1475 !  take account of the non-cartesian coordinates (d2cart).
1476    ABI_ALLOCATE(d2cart,(2,3,mpert,3,mpert))
1477    ABI_ALLOCATE(carflg,(3,mpert,3,mpert))
1478    ABI_ALLOCATE(d2matr,(2,3,mpert,3,mpert))
1479    outd2=1
1480    call status(0,dtfil%filstat,iexit,level,'call dfpt_gatherdy    ')
1481    call dfpt_gatherdy(becfrnl,dtset%berryopt,blkflg,carflg,&
1482 &   dyew,dyfrwf,dyfrx1,dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,&
1483 &   eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
1484 &   gprimd,dtset%mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,&
1485 &   rfasr,rfpert,rprimd,dtset%typat,ucvol,usevdw,psps%ziontypat)
1486 
1487    dscrpt=' Note : temporary (transfer) database '
1488 
1489 !  Initialize the header of the DDB file
1490    call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,&
1491 &   1,xred=xred,occ=occ,ngfft=ngfft)
1492 
1493 !  Open the formatted derivative database file, and write the header
1494    call status(0,dtfil%filstat,iexit,level,'call ddb_hdr_open_write')
1495    call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_ddb, dtfil%unddb)
1496 
1497    call ddb_hdr_free(ddb_hdr)
1498 
1499 !  Output of the dynamical matrix (master only)
1500    call status(0,dtfil%filstat,iexit,level,'call dfpt_dyout   ')
1501    call dfpt_dyout(becfrnl,dtset%berryopt,blkflg,carflg,dtfil%unddb,ddkfil,dyew,dyfrlo,&
1502 &   dyfrnl,dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,d2eig0,&
1503 &   d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,&
1504 &   eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
1505 &   has_full_piezo,has_allddk,ab_out,dtset%mband,mpert,natom,ntypat,&
1506 &   outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,dtset%prtvol,qphon,qzero,&
1507 &   dtset%typat,rfdir,rfpert,rfphon,rfstrs,psps%usepaw,usevdw,psps%ziontypat)
1508 
1509    close(dtfil%unddb)
1510 
1511 #ifdef HAVE_NETCDF
1512    ! Output dynamical matrix in NetCDF format.
1513    call crystal_init(dtset%amu_orig(:,1), Crystal, &
1514 &   dtset%spgroup, dtset%natom, dtset%npsp, psps%ntypat, &
1515 &   dtset%nsym, rprimd, dtset%typat, xred, dtset%ziontypat, dtset%znucl, 1, &
1516 &   dtset%nspden==2.and.dtset%nsppol==1, .false., hdr%title, &
1517 &   dtset%symrel, dtset%tnons, dtset%symafm)
1518 
1519    filename = strcat(dtfil%filnam_ds(4),"_DDB.nc")
1520 
1521    call outddbnc(filename, mpert, d2matr, blkflg, dtset%qptn, Crystal)
1522 
1523    call crystal_free(Crystal)
1524 
1525 #endif
1526 
1527 
1528 !  In case of phonons, diagonalize the dynamical matrix
1529    if(rfphon==1)then
1530 
1531 !    First, suppress the 'wings' elements,
1532 !    for which the diagonal element is not known
1533      call wings3(carflg,d2cart,mpert)
1534 
1535 !    Check the analyticity of the dynamical matrix
1536      analyt=0
1537      if (rfpert(natom+2)==0 .or. rfpert(natom+2)==2 .or. sumg0==1 ) analyt=1
1538 
1539 !    Diagonalize the analytic part
1540      ABI_ALLOCATE(displ,(2*3*natom*3*natom))
1541      ABI_ALLOCATE(eigval,(3*natom))
1542      ABI_ALLOCATE(eigvec,(2*3*natom*3*natom))
1543      ABI_ALLOCATE(phfrq,(3*natom))
1544      qphnrm=one
1545      call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,&
1546 &     dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,&
1547 &     dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol)
1548 
1549 !    Print the phonon frequencies
1550      call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon)
1551 
1552 !    Check the completeness of the dynamical matrix and eventually send a warning
1553      call chkph3(carflg,0,mpert,natom)
1554    end if ! end case of phonons
1555  end if !end me == 0
1556 
1557 !Compute the other terms for AHC dynamic and AHC full
1558  if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0.or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfuser==0 &
1559 & .and. rfmagn==0)) then
1560    if(rfphon==1) then ! AHC can only be computed in case of phonons
1561 
1562 !    Stuff for parallelism
1563      if(master /= me) then
1564        ABI_ALLOCATE(phfrq,(3*natom))
1565        ABI_ALLOCATE(displ,(2*3*natom*3*natom))
1566      end if
1567      call xmpi_bcast (phfrq,master,mpi_enreg%comm_cell,ierr) !Broadcast phfrq and displ
1568      call xmpi_bcast (displ,master,mpi_enreg%comm_cell,ierr) !to all processes
1569 
1570      if(dtset%ieig2rf == 3 .or. dtset%ieig2rf == 4 .or. dtset%ieig2rf == 5 ) then
1571        bdeigrf = dtset%bdeigrf
1572        if(dtset%bdeigrf == -1) bdeigrf = dtset%mband
1573 !      if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1574 !      &         (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
1575 !      write(std_out,*)'Reading the dense grid WF file'
1576 !      call wfk_read_eigenvalues(dtfil%fnameabi_wfkfine,eigenq_fine,hdr_fine,mpi_enreg%comm_world)
1577 !      ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband")
1578 !      endif
1579        if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%nsym==1)then
1580          write(std_out,*) 'Entering: eig2tot'
1581          if(dtset%smdelta>0)then
1582            if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1583 &           (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
1584              call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1585 &             eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1586 &             dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1587 &             occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine)
1588            else
1589              call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1590 &             eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1591 &             dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1592 &             occ_rbz_pert,hdr0,eigbrd)
1593            end if
1594          else
1595            if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1596 &           (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
1597              call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1598 &             eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1599 &             dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1600 &             occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine)
1601            else
1602              call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,&
1603 &             eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,&
1604 &             dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,&
1605 &             occ_rbz_pert,hdr0)
1606            end if
1607          end if
1608          write(std_out,*) 'Leaving: eig2tot'
1609        end if
1610      end if
1611      if (dtset%ieig2rf > 0) then
1612        ABI_DEALLOCATE(eigen0_pert)
1613        ABI_DEALLOCATE(eigenq_pert)
1614        ABI_DEALLOCATE(occ_rbz_pert)
1615        ABI_DEALLOCATE(eigen1_pert)
1616        call hdr_free(hdr0)
1617        if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
1618 &       (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) )  then
1619 !         call hdr_free(hdr0)
1620          call hdr_free(hdr_fine)
1621          ABI_DEALLOCATE(eigenq_fine)
1622        end if
1623      end if ! ieig2rf == 3  or %ieig2rf == 4 or %ieig2rf == 5
1624    end if ! rfphon==1
1625  end if
1626  if(dtset%efmas>0) then
1627    ABI_DEALLOCATE(eigen0_pert)
1628    ABI_DEALLOCATE(eigen1_pert)
1629  end if
1630  ABI_DEALLOCATE(doccde)
1631 
1632 
1633  if(me==0)then
1634    if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and.rfuser==0 &
1635 &   .and. rfmagn==0) )then
1636      if(rfphon==1)then
1637 !      Compute and print the T=0 Fan, and possibly DDW contributions to the eigenenergies.
1638        if(dtset%ieig2rf > 0) then
1639          write(message, '(80a,9a)' ) ('=',mu=1,80),ch10,ch10,&
1640 &         ' ---- T=0 shift of eigenenergies due to electron-phonon interation at q ---- ',ch10,&
1641 &         ' Warning : the total shift must be computed through anaddb,                  ',ch10,&
1642 &         ' here, only the contribution of one q point is printed.                      ',ch10,&
1643 &         ' Print first the electronic eigenvalues, then the q-dependent Fan shift of eigenvalues.'
1644          call wrtout(ab_out,message,'COLL')
1645          call wrtout(std_out,  message,'COLL')
1646 
1647          if(qeq0)then
1648            write(message, '(a)' )&
1649 &           ' Phonons at gamma, also compute the Diagonal Debye-Waller shift of eigenvalues.'
1650            call wrtout(ab_out,message,'COLL')
1651            call wrtout(std_out,message,'COLL')
1652          end if
1653 
1654          write(message, '(a)' ) ' '
1655          call wrtout(ab_out,message,'COLL')
1656          call wrtout(std_out,message,'COLL')
1657 
1658          call prteigrs(eigen0,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1659 &         dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,3,0,dtset%prtvol,&
1660 &         eigen0,zero,zero,dtset%wtk)
1661 
1662          write(message, '(a)' ) ch10
1663          call wrtout(ab_out,message,'COLL')
1664          call wrtout(std_out,message,'COLL')
1665 
1666 !        Compute and print Fan contribution
1667          ABI_ALLOCATE(eigen_fan,(dtset%mband*dtset%nkpt*dtset%nsppol))
1668          ABI_ALLOCATE(eigen_fan_mean,(dtset%mband*dtset%nkpt*dtset%nsppol))
1669          call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_fan,gprimd,&
1670 &         dtset%mband,natom,dtset%nkpt,dtset%nsppol,1,phfrq,dtset%prtvol)
1671          call eigen_meandege(eigen0,eigen_fan,eigen_fan_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2)
1672          call prteigrs(eigen_fan_mean,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1673 &         dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,5,0,dtset%prtvol,&
1674 &         eigen0,zero,zero,dtset%wtk)
1675 
1676          if(qeq0 .or. dtset%getgam_eig2nkq>0)then
1677 
1678            write(message, '(a)' ) ch10
1679            call wrtout(ab_out,message,'COLL')
1680            call wrtout(std_out,message,'COLL')
1681 
1682 !          Compute and print Diagonal Debye-Waller contribution
1683            ABI_ALLOCATE(eigen_ddw,(dtset%mband*dtset%nkpt*dtset%nsppol))
1684            ABI_ALLOCATE(eigen_ddw_mean,(dtset%mband*dtset%nkpt*dtset%nsppol))
1685            if(qeq0)then
1686              call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_ddw,gprimd,&
1687 &             dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol)
1688              if(results_respfn%gam_jdtset == -dtset%jdtset)then
1689                sz1=dtset%mband*dtset%nsppol
1690                sz2=natom*dim_eig2nkq
1691                ABI_ALLOCATE(results_respfn%gam_eig2nkq,(2,sz1,dtset%nkpt,3,natom,3,sz2))
1692                results_respfn%gam_eig2nkq(:,:,:,:,:,:,:)=eig2nkq(:,:,:,:,:,:,:)
1693                results_respfn%gam_jdtset=dtset%jdtset
1694              end if
1695            else if(dtset%getgam_eig2nkq>0)then
1696              if(results_respfn%gam_jdtset==dtset%getgam_eig2nkq)then
1697                call elph2_fanddw(dim_eig2nkq,displ,results_respfn%gam_eig2nkq,eigen_ddw,&
1698 &               gprimd,dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol)
1699              else
1700                write(message,'(a,i0,2a,i0,2a)')&
1701 &               'results_respfn%gam_jdtset=',results_respfn%gam_jdtset,ch10,&
1702 &               'dtset%getgam_eig2nkq=',dtset%getgam_eig2nkq,ch10,&
1703 &               'So, it seems eig2nkq at gamma has not yet been computed, while it is needed now.'
1704                MSG_BUG(message)
1705              end if
1706            end if
1707            call eigen_meandege(eigen0,eigen_ddw,eigen_ddw_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2)
1708            call prteigrs(eigen_ddw_mean,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1709 &           dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,6,0,dtset%prtvol,&
1710 &           eigen0,zero,zero,dtset%wtk)
1711 
1712            write(message, '(a)' ) ch10
1713            call wrtout(ab_out,message,'COLL')
1714            call wrtout(std_out,message,'COLL')
1715 
1716 !          Print sum of mean Fan and DDW
1717            ABI_ALLOCATE(eigen_fanddw,(dtset%mband*dtset%nkpt*dtset%nsppol))
1718            eigen_fanddw=eigen_fan_mean+eigen_ddw_mean
1719            call prteigrs(eigen_fanddw,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,&
1720 &           dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,7,0,dtset%prtvol,&
1721 &           eigen0,zero,zero,dtset%wtk)
1722 
1723            ABI_DEALLOCATE(eigen_ddw)
1724            ABI_DEALLOCATE(eigen_ddw_mean)
1725            ABI_DEALLOCATE(eigen_fanddw)
1726 
1727          end if
1728 
1729          ABI_DEALLOCATE(eigen_fan)
1730          ABI_DEALLOCATE(eigen_fan_mean)
1731        end if
1732 
1733 !      In case of a non-analytical part,
1734 !      get the phonon frequencies for three different directions (in cartesian coordinates)
1735        if(analyt==0)then
1736          qphnrm=zero
1737          do idir=1,3
1738 !          Need to know the corresponding dielectric constant
1739            if(carflg(idir,natom+2,idir,natom+2)==1)then
1740              qphon(:)=zero ; qphon(idir)=one
1741 !            Get the phonon frequencies
1742              call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,&
1743 &             dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,&
1744 &             dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol)
1745 !            Print the phonon frequencies
1746              call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon)
1747 !            Check the completeness of the dynamical matrix
1748 !            and eventually send a warning
1749              call chkph3(carflg,idir,mpert,natom)
1750            end if
1751          end do
1752          if (idir < 4) then
1753            qphon(idir)=zero
1754          end if
1755        end if
1756 
1757        ABI_DEALLOCATE(displ)
1758        ABI_DEALLOCATE(eigval)
1759        ABI_DEALLOCATE(eigvec)
1760        ABI_DEALLOCATE(phfrq)
1761      end if ! rfphon == 1
1762      ABI_DEALLOCATE(carflg)
1763      ABI_DEALLOCATE(d2cart)
1764      ABI_DEALLOCATE(d2matr)
1765    end if ! End condition on if.not.
1766  end if ! master node
1767 
1768 !Deallocate arrays
1769  if (allocated(displ)) then
1770    ABI_DEALLOCATE(displ)
1771  end if
1772  if (allocated(eigval)) then
1773    ABI_DEALLOCATE(eigval)
1774  end if
1775  if (allocated(eigvec)) then
1776    ABI_DEALLOCATE(eigvec)
1777  end if
1778  if (allocated(phfrq)) then
1779    ABI_DEALLOCATE(phfrq)
1780  end if
1781 
1782  ABI_DEALLOCATE(clflg)
1783  ABI_DEALLOCATE(atindx)
1784  ABI_DEALLOCATE(atindx1)
1785  ABI_DEALLOCATE(blkflg)
1786  ABI_DEALLOCATE(dyew)
1787  ABI_DEALLOCATE(dyewq0)
1788  ABI_DEALLOCATE(dyfrlo)
1789  ABI_DEALLOCATE(dyfrnl)
1790  ABI_DEALLOCATE(dyfrwf)
1791  ABI_DEALLOCATE(dyfrx1)
1792  ABI_DEALLOCATE(dyfrx2)
1793  ABI_DEALLOCATE(dyvdw)
1794  ABI_DEALLOCATE(d2bbb)
1795  ABI_DEALLOCATE(d2cart_bbb)
1796  ABI_DEALLOCATE(d2eig0)
1797  ABI_DEALLOCATE(d2k0)
1798  ABI_DEALLOCATE(d2lo)
1799  ABI_DEALLOCATE(d2loc0)
1800  ABI_DEALLOCATE(d2nfr)
1801  ABI_DEALLOCATE(d2nl)
1802  ABI_DEALLOCATE(d2nl0)
1803  ABI_DEALLOCATE(d2nl1)
1804  ABI_DEALLOCATE(d2ovl)
1805  ABI_DEALLOCATE(d2vn)
1806  ABI_DEALLOCATE(eigen0)
1807  ABI_DEALLOCATE(eig2nkq)
1808  ABI_DEALLOCATE(eigbrd)
1809  ABI_DEALLOCATE(eltcore)
1810  ABI_DEALLOCATE(elteew)
1811  ABI_DEALLOCATE(eltfrhar)
1812  ABI_DEALLOCATE(eltfrnl)
1813  ABI_DEALLOCATE(eltfrloc)
1814  ABI_DEALLOCATE(eltfrkin)
1815  ABI_DEALLOCATE(eltfrxc)
1816  ABI_DEALLOCATE(eltvdw)
1817  ABI_DEALLOCATE(becfrnl)
1818  ABI_DEALLOCATE(piezofrnl)
1819  call efmasdeg_free_array(efmasdeg)
1820  call efmasval_free_array(efmasval)
1821  ABI_DEALLOCATE(grxc)
1822  ABI_DEALLOCATE(indsym)
1823  ABI_DEALLOCATE(kxc)
1824  ABI_DEALLOCATE(nattyp)
1825  ABI_DEALLOCATE(pertsy)
1826  ABI_DEALLOCATE(rfpert)
1827  ABI_DEALLOCATE(rhog)
1828  ABI_DEALLOCATE(rhor)
1829  ABI_DEALLOCATE(symq)
1830  ABI_DEALLOCATE(symrec)
1831  ABI_DEALLOCATE(vtrial)
1832  ABI_DEALLOCATE(ylm)
1833  ABI_DEALLOCATE(ylmgr)
1834  call pawfgr_destroy(pawfgr)
1835  if (psps%usepaw==1) then
1836    call pawrhoij_free(pawrhoij)
1837    call paw_an_free(paw_an)
1838    call paw_ij_free(paw_ij)
1839    call pawfgrtab_free(pawfgrtab)
1840  end if
1841  ABI_DEALLOCATE(nhat)
1842  ABI_DEALLOCATE(nhatgr)
1843  ABI_DATATYPE_DEALLOCATE(pawrhoij)
1844  ABI_DATATYPE_DEALLOCATE(paw_an)
1845  ABI_DATATYPE_DEALLOCATE(paw_ij)
1846  ABI_DATATYPE_DEALLOCATE(pawfgrtab)
1847  if(rfphon==1.and.psps%n1xccc/=0)then
1848    ABI_DEALLOCATE(blkflgfrx1)
1849  end if
1850 
1851 !Clean the header
1852  call hdr_free(hdr)
1853 
1854 !Clean GPU data
1855 #if defined HAVE_GPU_CUDA
1856  if (dtset%use_gpu_cuda==1) then
1857    call dealloc_hamilt_gpu(0,dtset%use_gpu_cuda)
1858  end if
1859 #endif
1860 
1861  call status(0,dtfil%filstat,iexit,level,'exit')
1862 
1863  call timab(138,2,tsec)
1864  call timab(132,2,tsec)
1865 
1866  DBG_EXIT("COLL")
1867 
1868 end subroutine respfn

m_respfn_driver/wrtloctens [ Functions ]

[ Top ] [ m_respfn_driver ] [ Functions ]

NAME

 wrtloctens

FUNCTION

 Output of the localisation tensor

INPUTS

 blkflg = flags for each element of the 2DTE (=1 if computed)
 d2bbb  = band by band decomposition of second order derivatives
 d2nl   = non-local contributions to the 2DTEs
 mband  = maximum number of bands
 mpert  = maximum number of ipert
 natom  = number of atoms in unit cell
 prtbbb = if = 1, write the band by band decomposition of the localization tensor
 rprimd = dimensional primitive translations for real space (bohr)
 usepaw = flag for PAW

OUTPUT

  (only writing)

TODO

  The localization tensor cannot be defined in the metallic case. It should not be computed.

PARENTS

      respfn

CHILDREN

      wrtout

SOURCE

1903 subroutine wrtloctens(blkflg,d2bbb,d2nl,mband,mpert,natom,prtbbb,rprimd,usepaw)
1904 
1905 
1906 !This section has been created automatically by the script Abilint (TD).
1907 !Do not modify the following lines by hand.
1908 #undef ABI_FUNC
1909 #define ABI_FUNC 'wrtloctens'
1910 !End of the abilint section
1911 
1912  implicit none
1913 
1914 !Arguments ------------------------------------
1915 !scalars
1916  integer,intent(in) :: mband,mpert,natom,prtbbb,usepaw
1917 !arrays
1918  integer,intent(in) :: blkflg(3,mpert,3,mpert)
1919  real(dp),intent(in) :: rprimd(3,3)
1920  real(dp),intent(inout) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb)
1921  real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert)
1922 
1923 !Local variables ------------------------------
1924 !scalars
1925  integer :: flag,iband,idir,idir2,jband
1926  character(len=500) :: message
1927 !arrays
1928  real(dp) :: loctenscart(2,3,3)
1929  real(dp),allocatable :: loctenscart_bbb(:,:,:,:,:)
1930 
1931 ! *********************************************************************
1932 
1933 !This feature is disabled in the PAW case
1934  if (usepaw==1) return
1935 
1936  if(prtbbb==1)then
1937    ABI_ALLOCATE(loctenscart_bbb,(2,3,3,mband,mband*prtbbb))
1938    loctenscart_bbb(:,:,:,:,:)=zero
1939  end if
1940 
1941 !complete missing elements
1942  flag = 0
1943  do idir2 = 1,3
1944    do idir = 1,3
1945 
1946      if (blkflg(idir2,natom+1,idir,natom+1)==0) then
1947        if (blkflg(idir,natom+1,idir2,natom+1)==0) then
1948          flag = 1
1949        else
1950          d2nl(1,idir2,natom+1,idir,natom+1) = d2nl(1,idir,natom+1,idir2,natom+1)
1951          d2nl(2,idir2,natom+1,idir,natom+1) =-d2nl(2,idir,natom+1,idir2,natom+1)
1952          if(prtbbb==1)then
1953            d2bbb(1,idir2,idir,natom+1,:,:) = d2bbb(1,idir,idir2,natom+1,:,:)
1954            d2bbb(2,idir2,idir,natom+1,:,:) =-d2bbb(2,idir,idir2,natom+1,:,:)
1955          end if
1956 
1957        end if
1958      end if
1959 
1960    end do ! idir=1,3
1961  end do ! idir2=1,3
1962 
1963 !Transform the tensor to cartesian coordinates
1964 
1965  loctenscart(1,:,:) = matmul(rprimd,d2nl(1,:,natom+1,:,natom+1))
1966  loctenscart(2,:,:) = matmul(rprimd,d2nl(2,:,natom+1,:,natom+1))
1967 
1968  loctenscart(1,:,:) = matmul(loctenscart(1,:,:),transpose(rprimd))
1969  loctenscart(2,:,:) = matmul(loctenscart(2,:,:),transpose(rprimd))
1970 
1971  loctenscart(:,:,:) = loctenscart(:,:,:)/(two_pi**2)
1972 
1973  if (prtbbb == 1) then
1974 
1975    write(message,'(a,a)')ch10, &
1976 &   ' Band by band decomposition of the localisation tensor (bohr^2)'
1977    call wrtout(std_out,message,'COLL')
1978    call wrtout(ab_out,message,'COLL')
1979 
1980    do iband = 1,mband
1981      do jband = 1,mband
1982 
1983        loctenscart_bbb(1,:,:,iband,jband) = matmul(rprimd,d2bbb(1,:,:,natom+1,iband,jband))
1984        loctenscart_bbb(2,:,:,iband,jband) = matmul(rprimd,d2bbb(2,:,:,natom+1,iband,jband))
1985 
1986        loctenscart_bbb(1,:,:,iband,jband) = matmul(loctenscart_bbb(1,:,:,iband,jband),transpose(rprimd))
1987        loctenscart_bbb(2,:,:,iband,jband) = matmul(loctenscart_bbb(2,:,:,iband,jband),transpose(rprimd))
1988 
1989        loctenscart_bbb(:,:,:,iband,jband) = loctenscart_bbb(:,:,:,iband,jband)/(two_pi**2)
1990 
1991 
1992        write(message,'(a,a,i5,a,i5,a)')ch10, &
1993 &       ' Localisation tensor (bohr^2) for band ',iband,',',jband, &
1994 &       ' in cartesian coordinates'
1995        call wrtout(std_out,message,'COLL')
1996        call wrtout(ab_out,message,'COLL')
1997 
1998        write(ab_out,*)'     direction              matrix element'
1999        write(ab_out,*)'  alpha     beta       real part   imaginary part'
2000        do idir2 = 1,3
2001          do idir = 1,3
2002            write(ab_out,'(5x,i1,8x,i1,3x,2f16.10)')idir2,idir,&
2003 &           loctenscart_bbb(1,idir2,idir,iband,jband),&
2004 &           loctenscart_bbb(2,idir2,idir,iband,jband)
2005          end do
2006        end do
2007 
2008      end do !jband
2009    end do !iband
2010 
2011  end if  !prtbbb
2012 
2013  if (usepaw==0) then
2014    write(message,'(a,a,a,a)')ch10, &
2015 &   ' Total localisation tensor (bohr^2) in cartesian coordinates',ch10,&
2016 &   '  WARNING : still subject to testing - especially symmetries.'
2017  else
2018    write(message,'(a,a,a,a,a,a)')ch10, &
2019 &   ' Total localisation tensor (bohr^2) in cartesian coordinates',ch10,&
2020 &   '  WARNING : probably wrong for PAW (printing for testing purpose)',ch10,&
2021 &   '  WARNING : still subject to testing - especially symmetries.'
2022  end if
2023  call wrtout(std_out,message,'COLL')
2024  call wrtout(ab_out,message,'COLL')
2025 
2026  write(ab_out,*)'     direction              matrix element'
2027  write(ab_out,*)'  alpha     beta       real part   imaginary part'
2028  do idir2 = 1,3
2029    do idir = 1,3
2030      write(ab_out,'(5x,i1,8x,i1,3x,2f16.10)')idir2,idir,&
2031 &     loctenscart(1,idir2,idir),&
2032 &     loctenscart(2,idir2,idir)
2033    end do
2034  end do
2035 
2036  if (flag == 1) then
2037    write(message,'(6a)')ch10,&
2038 &   ' WARNING : Localization tensor calculation (this does not apply to other properties).',ch10,&
2039 &   '  Not all d/dk perturbations were computed. So the localization tensor in reciprocal space is incomplete,',ch10,&
2040 &   '  and transformation to cartesian coordinates may be wrong.'
2041    call wrtout(std_out,message,'COLL')
2042    call wrtout(ab_out,message,'COLL')
2043  end if
2044 
2045  if (prtbbb == 1) then
2046    ABI_DEALLOCATE(loctenscart_bbb)
2047  end if
2048 
2049 end subroutine wrtloctens