TABLE OF CONTENTS


ABINIT/dfpt_etot [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_etot

FUNCTION

 Assemble different contributions to the variational part of the
 2nd derivative of total energy

INPUTS

  berryopt= 4/14: electric field is on; berryopt = 6/7/16/17: electric displacement field is on;
  eberry=energy associated with Berry phase
  edocc=correction to 2nd-order total energy coming from changes of occupation
  ehart1=1st-order Hartree part of 2nd-order total energy
  eeig0=0th-order eigenenergies part of 2nd-order total energy
  eew=2nd derivative of Ewald energy (hartree)
  efrhar=contrib. from frozen-wavefunction, hartree energy, to the 2nd-derivative of total energy
  efrkin=contrib. from frozen-wavefunction, kinetic energy, to the 2nd-derivative of total energy
  efrloc=contrib. from frozen-wavefunction, local potential, to the 2nd-derivative of total energy
  efrnl=contribution from frozen-wavefunction, non-local potential, to the 2nd-derivative of total energy
  efrx1=contrib. from frozen-wavefunction, xc core correction(1), to the 2nd-derivative of total energy
  efrx2=contribution from frozen-wavefunction, xc core correction(2),
           to the second-derivative of total energy.
  ek0=0th-order kinetic energy part of 2nd-order total energy.
  ek1=1st-order kinetic energy part of 2nd-order total energy.
  eii=2nd derivative of pseudopotential core energy (hartree)
  eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
  elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
  enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
  enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
  epaw1=1st-order PAW on-sitew part of 2nd-order total energy.
  evdw=DFT-D semi-empirical part of 2nd-order total energy
  exc1=1st-order exchange-correlation part of 2nd-order total energy
  ipert=type of the perturbation
  natom=number of atoms
  optene=option for the computation of 2nd-order total energy
         (-1=no computation; 0=direct scheme; 1=double-counting scheme)

OUTPUT

  deltae=change in energy between the previous and present SCF cycle
         and previous SCF cycle.
  etotal=2nd-order total energy
  evar=variational part of the 2nd-order total energy

SIDE EFFECTS

 input/output
 elast=previous value of the 2nd-order total energy, needed to compute deltae,
      then updated (cannot simply be saved, because set to zero
      at each new call of dfpt_scfcv).

PARENTS

      dfpt_scfcv

CHILDREN

SOURCE

1595 subroutine dfpt_etot(berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,&
1596 &                efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,&
1597 &                enl0,enl1,epaw1,etotal,evar,evdw,exc1,ipert,natom,optene)
1598 
1599 
1600 !This section has been created automatically by the script Abilint (TD).
1601 !Do not modify the following lines by hand.
1602 #undef ABI_FUNC
1603 #define ABI_FUNC 'dfpt_etot'
1604 !End of the abilint section
1605 
1606  implicit none
1607 
1608 !Arguments ------------------------------------
1609 !scalars
1610  integer,intent(in) :: berryopt,ipert,natom,optene
1611  real(dp),intent(in) :: eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1
1612  real(dp),intent(in) :: efrx2,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,epaw1
1613  real(dp),intent(in) :: evdw,exc1
1614  real(dp),intent(inout) :: elast
1615  real(dp),intent(out) :: deltae,etotal,evar
1616 
1617 !Local variables-------------------------------
1618 !scalars
1619 ! character(len=500) :: message
1620 
1621 ! *********************************************************************
1622 
1623  if (optene==1) then
1624    MSG_BUG('Double-counting scheme not yet allowed!')
1625  end if
1626 
1627  if (optene>-1) then
1628 
1629 !  Compute 2nd-order variational energy by direct scheme
1630    if (optene==0) then
1631 
1632 !    Atomic displ. perturbation
1633      if ( ipert>=1 .and. ipert<=natom  ) then
1634        evar=ek0+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1
1635 
1636      else if (ipert==natom+1) then
1637        evar=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1
1638 
1639      else if (ipert==natom+10 .or. ipert==natom+11) then
1640        evar=ek0+edocc+eeig0+eloc0+enl0+ek1 ! here ek1 contains a lot of contributions
1641 
1642 !      For ipert==natom+2, some contributions vanishes, noticeably ek1
1643      else if (ipert==natom+2) then
1644        evar=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+epaw1
1645 
1646 !      All terms enter for strain perturbation
1647      else if ( ipert==natom+3 .or. ipert==natom+4 ) then
1648        evar=ek0+ek1+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1
1649 
1650 !    terms for Zeeman perturbation, SPr 2deb
1651      else if ( ipert==natom+5 ) then
1652        evar=ek0+edocc+eeig0+eloc0+ehart1+exc1+enl0+epaw1
1653      end if
1654    end if
1655 
1656 !  Compute energy residual
1657    deltae=evar-elast
1658    elast=evar
1659 
1660 !  Compute 2nd-order total energy by direct scheme
1661    if (optene==0) then
1662      if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17) then
1663        if (ipert<=natom) then
1664          etotal=evar+eew+evdw+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2+two*eberry
1665        else if (ipert==natom+2) then
1666          etotal=half*evar+eew+evdw+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2+two*eberry
1667        end if
1668      else
1669        if (ipert/=natom+10 .and. ipert/=natom+11) then
1670          etotal=evar+eew+evdw+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2
1671        else
1672          etotal=evar ! For 2nd order sternheimer equations, the total (4th order) energy is not used (yet)
1673        end if
1674      end if
1675    end if
1676 
1677  end if
1678 
1679 end subroutine dfpt_etot

ABINIT/dfpt_newvtr [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_newvtr

FUNCTION

 Compute new first-order trial potential by mixing new and old values.
 First, compute preconditioned residual first-order potential.
 Then, call one of the self-consistency drivers, and  update vtrial.

INPUTS

  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  dielar(7)=input parameters for dielectric matrix:
                diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
  dtset <type(dataset_type)>=all input variables in this dataset
   | isecur=level of security of the computation
   | mffmem=governs the number of FFT arrays which are fit in core memory
   |          it is either 1, in which case the array f_fftgr is used,
   |          or 0, in which case the array f_fftgr_disk is used
   | natom=number of atoms
   | nspden=number of spin-density components
   | paral_kgb=option for (kpt,g vectors,bands) parallelism
   | pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part
  etotal=the total energy obtained from the input vtrial
  ffttomix(nfft*(1-nfftmix/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
  initialized= if 0 the initialization of the RF run is not yet finished
   iscf=( <= 0 =>non-SCF), >0 => SCF)
    iscf =1 => determination of the largest eigenvalue of the SCF cycle
    iscf =2 => SCF cycle, simple mixing
    iscf =3 => SCF cycle, Anderson mixing
    iscf =4 => SCF cycle, Anderson mixing (order 2)
    iscf =5 => SCF cycle, CG based on the minimization of the energy
    iscf =7 => SCF cycle, Pulay mixing
  ispmix=1 if mixing is done in real space, 2 if mixing is done in reciprocal space
  istep= number of the step in the SCF cycle
  mixtofft(nfftmix*(1-nfftmix/nfft))=Index of the points of the FFT grid used for mixing (coarse) on the FFT (fine) grid
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  nfft=(effective) number of FFT grid points (for this processor)
  nfftmix=dimension of FFT grid used to mix the densities (used in PAW only)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ngfftmix(18)=contain all needed information about 3D FFT, for the grid corresponding to nfftmix
  npawmix=-PAW only- number of spherical part elements to be mixed
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
                                         Use here rhoij residuals (and gradients)
  rhor(cplex*nfft,nspden)=array for 1st-order electron density
    in electrons/bohr**3.
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vresid(cplex*nfft,nspden)=array for the residual of the potential
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.

SIDE EFFECTS

  vtrial(cplex*nfft,nspden)= at input, it is the "in" trial potential that gave vresid=(v_out-v_in)
       at output, it is an updated "mixed" trial potential
  ==== if usepaw==1
    pawrhoij(natom)%nrhoijsel,rhoijselect,rhoijp= several arrays
                containing new values of rhoij (augmentation occupancies)

NOTES

  In case of PAW calculations:
    Computations are done either on the fine FFT grid or the coarse grid (depending on dtset%pawmixdg)
    All variables (nfft,ngfft,mgfft) refer to the fine FFT grid.
    All arrays (densities/potentials...) are computed on this fine FFT grid.
    Developpers have to be careful when introducing others arrays:
      they have to be stored on the fine FFT grid.
  In case of norm-conserving calculations the FFT grid is the usual FFT grid.

  Subtility in PAW and non-collinear magnetism:
    Potentials are stored in (up-up,dn-dn,Re[up-dn],Im[up-dn]) format
    On-site occupancies (rhoij) are stored in (n,mx,my,mz)
    This is compatible provided that the mixing factors for n and m are identical
    and that the residual is not a combination of V_res and rhoij_res (pawoptmix=0).

PARENTS

      dfpt_scfcv

CHILDREN

      ab7_mixing_copy_current_step,ab7_mixing_eval,ab7_mixing_eval_allocate
      ab7_mixing_eval_deallocate,fourdp,metric,moddiel,timab

SOURCE

1952 subroutine dfpt_newvtr(cplex,dbl_nnsclo,dielar,dtset,etotal,ffttomix,&
1953 &          initialized,iscf,ispmix,istep,mix,mixtofft,&
1954 &          mpi_enreg,my_natom,nfft,nfftmix,ngfft,ngfftmix,npawmix,pawrhoij,&
1955 &          qphon,rhor,rprimd,usepaw,vresid,vtrial)
1956 
1957 
1958 !This section has been created automatically by the script Abilint (TD).
1959 !Do not modify the following lines by hand.
1960 #undef ABI_FUNC
1961 #define ABI_FUNC 'dfpt_newvtr'
1962 !End of the abilint section
1963 
1964  implicit none
1965 
1966 !Arguments-------------------------------
1967 !scalars
1968  integer,intent(in) :: cplex,initialized,iscf,ispmix,istep,my_natom,nfft
1969  integer,intent(in) :: nfftmix,npawmix,usepaw
1970  integer,intent(inout) :: dbl_nnsclo !vz_i
1971  real(dp),intent(in) :: etotal
1972  type(MPI_type),intent(in) :: mpi_enreg
1973  type(ab7_mixing_object), intent(inout) :: mix
1974  type(dataset_type),intent(in) :: dtset
1975 !arrays
1976  integer,intent(in) :: ffttomix(nfft*(1-nfftmix/nfft))
1977  integer,intent(in) :: mixtofft(nfftmix*(1-nfftmix/nfft)),ngfft(18)
1978  integer,intent(in) :: ngfftmix(18)
1979  real(dp),intent(in) :: dielar(7),qphon(3)
1980  real(dp), intent(in), target :: rhor(cplex*nfft,dtset%nspden)
1981  real(dp),intent(in) :: rprimd(3,3)
1982  real(dp),intent(inout) :: vresid(cplex*nfft,dtset%nspden)
1983  real(dp),intent(inout) :: vtrial(cplex*nfft,dtset%nspden)
1984  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
1985 
1986 !Local variables-------------------------------
1987 !scalars
1988  integer :: cplex_mix,cplex_rhoij,dplex,i_vresid1,i_vrespc1,iatom,ifft,indx
1989  integer :: irhoij,ispden,jfft,jrhoij,klmn,kmix,moved_atm_inside,nfftot,nselect
1990  integer :: mpicomm,errid
1991  logical :: mpi_summarize,reset
1992  real(dp) :: fact,mixfac,mixfac_eff,mixfacmag,ucvol
1993  character(len=500) :: msg
1994 !arrays
1995  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2)
1996  real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:)
1997  real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:)
1998  real(dp), pointer :: vtrial0(:,:),vpaw(:)
1999  real(dp),allocatable :: vtrialg(:,:,:)
2000 
2001 ! *************************************************************************
2002 
2003  DBG_ENTER("COLL")
2004 
2005  call timab(158,1,tsec)
2006 
2007 !Compatibility tests
2008  if(usepaw==1) then
2009    if(dtset%nspden==4.and.dtset%pawoptmix==1) then
2010      MSG_ERROR('pawoptmix=1 is not compatible with nspden=4 !')
2011    end if
2012    if (my_natom>0) then
2013      if (pawrhoij(1)%cplex<cplex) then
2014        MSG_ERROR('pawrhoij()%cplex must be >=cplex !')
2015      end if
2016    end if
2017  end if
2018 
2019  nfftot=ngfft(1)*ngfft(2)*ngfft(3)
2020  cplex_mix=max(cplex,ispmix)
2021  if (usepaw==1.and.my_natom>0) cplex_rhoij=pawrhoij(1)%cplex
2022 
2023 !Compute different geometric tensor, as well as ucvol, from rprimd
2024  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2025  moved_atm_inside=0
2026 
2027 !Select components of potential to be mixed
2028  ABI_ALLOCATE(vtrial0,(cplex_mix*nfftmix,dtset%nspden))
2029  ABI_ALLOCATE(vresid0,(cplex_mix*nfftmix,dtset%nspden))
2030  if (ispmix==1.and.nfft==nfftmix) then
2031    vtrial0=vtrial;vresid0=vresid
2032  else if (nfft==nfftmix) then
2033    do ispden=1,dtset%nspden
2034      call fourdp(cplex,vtrial0(:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
2035      call fourdp(cplex,vresid0(:,ispden),vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
2036    end do
2037  else
2038    ABI_ALLOCATE(vtrialg,(2,nfft,dtset%nspden))
2039    ABI_ALLOCATE(vreswk,(2,nfft))
2040    do ispden=1,dtset%nspden
2041      fact=dielar(4);if (ispden>1) fact=dielar(7)
2042      call fourdp(cplex,vtrialg(:,:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
2043      call fourdp(cplex,vreswk,vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
2044      do ifft=1,nfft
2045        if (ffttomix(ifft)>0) then
2046          jfft=2*ffttomix(ifft)
2047          vtrial0(jfft-1,ispden)=vtrialg(1,ifft,ispden)
2048          vtrial0(jfft  ,ispden)=vtrialg(2,ifft,ispden)
2049          vresid0(jfft-1,ispden)=vreswk(1,ifft)
2050          vresid0(jfft  ,ispden)=vreswk(2,ifft)
2051        else
2052          vtrialg(:,ifft,ispden)=vtrialg(:,ifft,ispden)+fact*vreswk(:,ifft)
2053        end if
2054      end do
2055    end do
2056    ABI_DEALLOCATE(vreswk)
2057  end if
2058 
2059 !Precondition the potential residual:
2060 !Use a model dielectric function preconditioning, or simple mixing
2061  ABI_ALLOCATE(vrespc,(cplex_mix*nfftmix,dtset%nspden))
2062  call moddiel(cplex_mix,dielar,mpi_enreg,nfftmix,ngfftmix,dtset%nspden,ispmix,0,&
2063 & dtset%paral_kgb,qphon,rprimd,vresid0,vrespc)
2064 
2065 !PAW only : precondition the rhoij quantities (augmentation occupancies) residuals.
2066 !Use a simple preconditionning with the same mixing factor
2067 !as the model dielectric function.
2068  if (usepaw==1.and.my_natom>0) then
2069    ABI_ALLOCATE(rhoijrespc,(npawmix))
2070    mixfac=dielar(4);mixfacmag=abs(dielar(7))
2071    if (cplex_rhoij==1) then
2072      indx=0
2073      do iatom=1,my_natom
2074        do ispden=1,pawrhoij(iatom)%nspden
2075          mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
2076          do kmix=1,pawrhoij(iatom)%lmnmix_sz
2077            indx=indx+1;klmn=pawrhoij(iatom)%kpawmix(kmix)
2078            rhoijrespc(indx)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden)
2079          end do
2080        end do
2081      end do
2082    else
2083      indx=-1
2084      do iatom=1,my_natom
2085        do ispden=1,pawrhoij(iatom)%nspden
2086          mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag
2087          do kmix=1,pawrhoij(iatom)%lmnmix_sz
2088            indx=indx+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1
2089            rhoijrespc(indx:indx+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden)
2090          end do
2091        end do
2092      end do
2093    end if
2094  end if
2095 
2096 !------Compute new vtrial
2097 
2098  i_vresid1=mix%i_vresid(1)
2099  i_vrespc1=mix%i_vrespc(1)
2100 
2101 !Initialise working arrays for the mixing object.
2102  call ab7_mixing_eval_allocate(mix, istep)
2103 
2104 !Copy current step arrays.
2105  call ab7_mixing_copy_current_step(mix, vresid0, errid, msg, arr_respc = vrespc)
2106 
2107  if (errid /= AB7_NO_ERROR) then
2108    MSG_ERROR(msg)
2109  end if
2110 
2111  ABI_DEALLOCATE(vrespc)
2112  ABI_DEALLOCATE(vresid0)
2113 
2114 !PAW: either use the array f_paw or the array f_paw_disk
2115  ABI_ALLOCATE(vpaw,(npawmix*usepaw))
2116  if (usepaw==1.and.my_natom>0) then
2117    dplex=cplex_rhoij-1
2118    indx=-dplex
2119    do iatom=1,my_natom
2120      do ispden=1,pawrhoij(iatom)%nspden
2121        ABI_ALLOCATE(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,1))
2122        rhoijtmp=zero
2123        jrhoij=1
2124        do irhoij=1,pawrhoij(iatom)%nrhoijsel
2125          klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
2126          rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
2127          jrhoij=jrhoij+cplex_rhoij
2128        end do
2129        do kmix=1,pawrhoij(iatom)%lmnmix_sz
2130          indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex
2131          vpaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
2132          mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
2133          mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex)
2134        end do
2135        ABI_DEALLOCATE(rhoijtmp)
2136      end do
2137    end do
2138  end if
2139 
2140 !Unlike for GS, no need to modify the mean of vtrial
2141 
2142  mpicomm=0;mpi_summarize=.false.
2143  reset=.false.;if (initialized==0) reset=.true.
2144  call ab7_mixing_eval(mix, vtrial0, istep, nfftot, ucvol, &
2145 & mpicomm, mpi_summarize, errid, msg, &
2146 & reset = reset, isecur = dtset%isecur, &
2147 & pawopt = dtset%pawoptmix, response = 1, pawarr = vpaw, &
2148 & etotal = etotal, potden = rhor, comm_atom=mpi_enreg%comm_atom)
2149 
2150  if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
2151    dbl_nnsclo = 1
2152  else if (errid /= AB7_NO_ERROR) then
2153    ! MG FIXME, Why this?
2154    ! One should propagate the error so that we can handle it
2155    ! in the caller!
2156    MSG_ERROR(msg)
2157  end if
2158 
2159 !Do here the mixing of the potential
2160  if(iscf==2 .or. iscf==3 .or. iscf==7)then
2161 !  PAW: restore rhoij from compact storage
2162    if (usepaw==1.and.my_natom>0) then
2163      dplex=cplex_rhoij-1
2164      indx=-dplex
2165      do iatom=1,my_natom
2166        ABI_ALLOCATE(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
2167        rhoijtmp=zero
2168        if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
2169          do ispden=1,pawrhoij(iatom)%nspden
2170            jrhoij=1
2171            do irhoij=1,pawrhoij(iatom)%nrhoijsel
2172              klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
2173              rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
2174              jrhoij=jrhoij+cplex_rhoij
2175            end do
2176          end do
2177        end if
2178        do ispden=1,pawrhoij(iatom)%nspden
2179          do kmix=1,pawrhoij(iatom)%lmnmix_sz
2180            indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex
2181            rhoijtmp(klmn:klmn+dplex,ispden)=vpaw(indx:indx+dplex)
2182          end do
2183        end do
2184        nselect=0
2185        do klmn=1,pawrhoij(iatom)%lmn2_size
2186          if (any(abs(rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,:))>tol10)) then
2187            nselect=nselect+1
2188            pawrhoij(iatom)%rhoijselect(nselect)=klmn
2189            do ispden=1,pawrhoij(iatom)%nspden
2190              pawrhoij(iatom)%rhoijp(cplex_rhoij*nselect-dplex:cplex_rhoij*nselect,ispden)=&
2191 &             rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,ispden)
2192            end do
2193          end if
2194        end do
2195        pawrhoij(iatom)%nrhoijsel=nselect
2196        ABI_DEALLOCATE(rhoijtmp)
2197      end do
2198    end if
2199 
2200  else if(iscf==5 .or. iscf==6)then
2201    if(ispmix/=1) then
2202      MSG_ERROR('Mixing on reciprocal space not allowed with iscf=5 or 6.')
2203    end if
2204 !  PAW: apply a simple mixing to rhoij (this is temporary)
2205    if (usepaw==1.and.my_natom>0) then
2206      indx=1-cplex_rhoij
2207      do iatom=1,my_natom
2208        ABI_ALLOCATE(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
2209        rhoijtmp=zero
2210        if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
2211          do ispden=1,pawrhoij(iatom)%nspden
2212            do kmix=1,pawrhoij(iatom)%lmnmix_sz
2213              indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex
2214              rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) &
2215 &             -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
2216            end do
2217          end do
2218        end if
2219        do ispden=1,pawrhoij(iatom)%nspden
2220          jrhoij=1
2221          do irhoij=1,pawrhoij(iatom)%nrhoijsel
2222            klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
2223            rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
2224 &           +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
2225            jrhoij=jrhoij+cplex_rhoij
2226          end do
2227        end do
2228        nselect=0
2229        do klmn=1,pawrhoij(iatom)%lmn2_size
2230          if (any(abs(rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,:))>tol10)) then
2231            nselect=nselect+1
2232            pawrhoij(iatom)%rhoijselect(nselect)=klmn
2233            do ispden=1,pawrhoij(iatom)%nspden
2234              pawrhoij(iatom)%rhoijp(cplex_rhoij*nselect-dplex:cplex_rhoij*nselect,ispden)=&
2235 &             rhoijtmp(cplex_rhoij*klmn-dplex:cplex_rhoij*klmn,ispden)
2236            end do
2237          end if
2238        end do
2239        pawrhoij(iatom)%nrhoijsel=nselect
2240        ABI_DEALLOCATE(rhoijtmp)
2241      end do
2242    end if
2243  end if
2244 
2245  ABI_DEALLOCATE(vpaw)
2246  if (usepaw==1.and.my_natom>0)  then
2247    ABI_DEALLOCATE(rhoijrespc)
2248  end if
2249 
2250 !Eventually write the data on disk and deallocate f_fftgr_disk
2251  call ab7_mixing_eval_deallocate(mix)
2252 
2253 !Restore potential
2254  if (ispmix==1.and.nfft==nfftmix) then
2255    vtrial=vtrial0
2256  else if (nfft==nfftmix) then
2257    do ispden=1,dtset%nspden
2258      call fourdp(cplex,vtrial0(:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
2259    end do
2260  else
2261    do ispden=1,dtset%nspden
2262      do ifft=1,nfftmix
2263        jfft=mixtofft(ifft)
2264        vtrialg(1,jfft,ispden)=vtrial0(2*ifft-1,ispden)
2265        vtrialg(2,jfft,ispden)=vtrial0(2*ifft  ,ispden)
2266      end do
2267      call fourdp(cplex,vtrialg(:,:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
2268    end do
2269    ABI_DEALLOCATE(vtrialg)
2270  end if
2271  ABI_DEALLOCATE(vtrial0)
2272 
2273  call timab(158,2,tsec)
2274 
2275  DBG_ENTER("COLL")
2276 
2277 end subroutine dfpt_newvtr

ABINIT/dfpt_nselt [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nselt

FUNCTION

 This routine compute the non-stationary expression for the
 second derivative of the total energy, wrt strain for a whole row of
 mixed strain derivatives.

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL,
    if 2, COMPLEX
  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)
  effmass_free=effective mass for electrons (1. in common case)
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  idir=direction of the perturbation
  ipert=type of the perturbation
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the
     storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ
  kxc(nfft,nkxc)=exchange and correlation kernel
  mband=maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem =number of k points treated by this node.
  mk1mem =number of k points treated by this node  (RF data).
  mpert =maximum number of ipert
  mpi_enreg=information about MPI parallelization
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  maximum dimension for q points in grids for nonlocal form factors
  natom=number of atoms in cell.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
  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
  nkpt_rbz=number of k points in the reduced BZ for this perturbation
  nkxc=second dimension of the kxc array. If /=0,
   the exchange-correlation kernel must be computed.
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym1=number of symmetry elements in space group consistent with
    perturbation
  ntypat=number of types of atoms in unit cell.
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band
   and k in the reduced Brillouin zone (usually =2)
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  prtbbb=if 1, band-by-band decomposition (also dim of d2bbb)
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  qphon(3)=reduced coordinates for the phonon wavelength
  rhog(2,nfft)=array for Fourier transform of GS electron density
  rhor(nfft,nspden)=GS electron density in electrons/bohr**3.
  rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3.
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  typat(natom)=type integer for each atom in cell
  ucvol=unit cell volume in bohr**3.
  wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang)= real spherical harmonics for each G and k+q point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point
  ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k+g point

OUTPUT

  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
       second order derivatives
  d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
  d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs

PARENTS

      dfpt_scfcv

CHILDREN

      destroy_hamiltonian,dfpt_mkcore,dfpt_mkvxcstr,dfpt_nsteltwf,dotprod_vn
      hartrestr,init_hamiltonian,stresssym,vlocalstr,wrtout,xmpi_barrier

SOURCE

2372 subroutine dfpt_nselt(blkflg,cg,cg1,cplex,&
2373 & d2bbb,d2lo,d2nl,ecut,ecutsm,effmass_free,&
2374 & gmet,gprimd,gsqcut,idir,&
2375 & ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband,mgfft,&
2376 & mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,mpw1,&
2377 & natom,nband_rbz,nfft,ngfft,&
2378 & nkpt_rbz,nkxc,nloalg,npwarr,npwar1,nspden,nspinor,nsppol,&
2379 & nsym1,ntypat,occ_rbz,&
2380 & paral_kgb,ph1d,prtbbb,psps,qphon,rhog,&
2381 & rhor,rhor1,rmet,rprimd,symrc1,typat,ucvol,&
2382 & wtk_rbz,&
2383 & xred,ylm,ylm1,ylmgr,ylmgr1)
2384 
2385  use m_hamiltonian,only : init_hamiltonian,destroy_hamiltonian,gs_hamiltonian_type
2386  use m_spacepar,   only : hartrestr
2387  use m_mkcore,      only : dfpt_mkcore
2388 
2389 !This section has been created automatically by the script Abilint (TD).
2390 !Do not modify the following lines by hand.
2391 #undef ABI_FUNC
2392 #define ABI_FUNC 'dfpt_nselt'
2393 !End of the abilint section
2394 
2395  implicit none
2396 
2397 !Arguments -------------------------------
2398 !scalars
2399  integer,intent(in) :: cplex,idir,ipert,mband,mgfft,mk1mem
2400  integer,intent(in) :: mkmem,mpert,mpsang,mpw,mpw1,natom,nfft,nkpt_rbz
2401  integer,intent(in) :: nkxc,nspden,nspinor,nsppol,nsym1,ntypat
2402  integer,intent(in) :: paral_kgb,prtbbb
2403  real(dp),intent(in) :: ecut,ecutsm,effmass_free,gsqcut,ucvol
2404  type(MPI_type),intent(in) :: mpi_enreg
2405  type(pseudopotential_type),intent(in) :: psps
2406 !arrays
2407  integer,intent(in) :: istwfk_rbz(nkpt_rbz)
2408  integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
2409  integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18)
2410  integer,intent(in) :: nloalg(3),npwar1(nkpt_rbz),npwarr(nkpt_rbz)
2411  integer,intent(in) :: symrc1(3,3,nsym1),typat(natom)
2412  integer,intent(inout) :: blkflg(3,mpert,3,mpert)
2413  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
2414  real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)
2415  real(dp),intent(in) :: gmet(3,3)
2416  real(dp),intent(in) :: gprimd(3,3),kpt_rbz(3,nkpt_rbz),kxc(nfft,nkxc)
2417  real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
2418  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qphon(3),rhog(2,nfft)
2419  real(dp),intent(in) :: rhor(nfft,nspden)
2420  real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3)
2421  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom)
2422  real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
2423  real(dp),intent(in) :: ylm1(mpw1*mk1mem,mpsang*mpsang*psps%useylm)
2424  real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm)
2425  real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*psps%useylm)
2426  real(dp),intent(out) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb)
2427  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert)
2428  real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert)
2429 
2430 !Local variables-------------------------------
2431 !scalars
2432  integer :: ban2tot,bantot,bd2tot_index,bdtot_index
2433  integer :: icg,icg1,idir1,ifft,ii,ikg,ikg1,ikpt,comm
2434  integer :: ilm,ipert1,ispden,isppol,istr1,istwf_k
2435  integer :: mbd2kpsp,mbdkpsp,me,n1,n2,n3,n3xccc,n4,n5,n6
2436  integer :: nband_k,nfftot,npw1_k,npw_k,option
2437  real(dp) :: doti,dotr
2438  real(dp) :: wtk_k
2439  character(len=500) :: message
2440  type(gs_hamiltonian_type) :: gs_hamk
2441 !arrays
2442  integer :: ikpt_fbz(3)
2443  integer,allocatable :: kg1_k(:,:),kg_k(:,:)
2444  real(dp) :: kpoint(3),restr(6),dummy(0,0)
2445  real(dp),allocatable :: d2bbb_k(:,:,:,:),d2nl_k(:,:,:)
2446  real(dp),allocatable :: occ_k(:)
2447  real(dp),allocatable :: vhartr01(:),vpsp1(:),vxc1(:,:),xccc3d1(:),ylm1_k(:,:)
2448  real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:),ylmgr_k(:,:,:)
2449  type(pawtab_type) :: pawtab_dum(0)
2450 
2451 ! *********************************************************************
2452 
2453 !Init me
2454  comm = mpi_enreg%comm_cell
2455  me   = mpi_enreg%me_kpt
2456 
2457 !Unit numbers
2458 
2459 !Zero only portion of nonlocal matrix to be computed here
2460  d2nl(:,:,natom+3:natom+4,idir,ipert)=zero
2461  bdtot_index=0
2462  bd2tot_index=0
2463  icg=0
2464  icg1=0
2465  mbdkpsp=mband*nkpt_rbz*nsppol
2466  mbd2kpsp=2*mband**2*nkpt_rbz*nsppol
2467 
2468 !Update list of computed matrix elements
2469  if((ipert==natom+3) .or. (ipert==natom+4)) then
2470 !  Eventually expand when strain coupling to other perturbations is
2471 !  implemented
2472    do ipert1=natom+3,natom+4
2473      do idir1=1,3
2474        blkflg(idir1,ipert1,idir,ipert)=1
2475      end do
2476    end do
2477  end if
2478 
2479 !allocate(enl1nk(mbdkpsp))
2480  ABI_ALLOCATE(d2bbb_k,(2,3,mband,mband*prtbbb))
2481  ABI_ALLOCATE(d2nl_k,(2,3,mpert))
2482 
2483 
2484  ABI_ALLOCATE(kg_k,(3,mpw))
2485  ABI_ALLOCATE(kg1_k,(3,mpw1))
2486 
2487 
2488  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2489  n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
2490  nfftot=n1*n2*n3
2491 
2492 !Initialize Hamiltonian (k-independent terms) - NCPP only
2493  call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,&
2494 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,ph1d=ph1d)
2495 
2496  bantot = 0
2497  ban2tot = 0
2498 
2499 !LOOP OVER SPINS
2500  do isppol=1,nsppol
2501 
2502    if (nsppol/=1) then
2503      write(message,*)' ****  In dfpt_nselt for isppol=',isppol
2504      call wrtout(std_out,message,'COLL')
2505    end if
2506 
2507    ikg=0
2508    ikg1=0
2509 
2510    ikpt_fbz(1:3)=0
2511 
2512 !  BIG FAT k POINT LOOP
2513    do ikpt=1,nkpt_rbz
2514 
2515      nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
2516      istwf_k=istwfk_rbz(ikpt)
2517      npw_k=npwarr(ikpt)
2518      npw1_k=npwar1(ikpt)
2519      kpoint(:)=kpt_rbz(:,ikpt)
2520 
2521      bantot = bantot + nband_k
2522      ban2tot = ban2tot + 2*nband_k**2
2523 
2524      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
2525        bdtot_index=bdtot_index+nband_k
2526        bd2tot_index=bd2tot_index+2*nband_k**2
2527 !      Skip the rest of the k-point loop
2528        cycle
2529      end if
2530 
2531      ABI_ALLOCATE(occ_k,(nband_k))
2532 
2533      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang))
2534      ABI_ALLOCATE(ylm1_k,(npw1_k,mpsang*mpsang))
2535      if (ipert==natom+3.or.ipert==natom+4) then
2536        ABI_ALLOCATE(ylmgr_k,(npw_k,3,mpsang*mpsang))
2537        ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,mpsang*mpsang))
2538      end if
2539 
2540 !    enl1_k(:)=zero
2541      d2nl_k(:,:,:)=zero
2542      if(prtbbb==1)d2bbb_k(:,:,:,:)=zero
2543      occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
2544 
2545      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2546      if (psps%useylm==1) then
2547        do ilm=1,mpsang*mpsang
2548          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
2549        end do
2550        if (ipert==natom+3.or.ipert==natom+4) then
2551          do ilm=1,mpsang*mpsang
2552            do ii=1,3
2553              ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
2554            end do
2555          end do
2556        end if
2557      end if
2558 
2559      wtk_k=wtk_rbz(ikpt)
2560 
2561      kg1_k(:,:) = 0
2562 
2563      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
2564      if (psps%useylm==1) then
2565        do ilm=1,mpsang*mpsang
2566          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
2567        end do
2568        if (ipert==natom+3.or.ipert==natom+4) then
2569          do ilm=1,mpsang*mpsang
2570            do ii=1,3
2571              ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
2572            end do
2573          end do
2574        end if
2575      end if
2576 
2577 !    Compute the eigenvalues, wavefunction,
2578 !    contributions to kinetic energy, nonlocal energy, forces,
2579 !    and update of rhor1 to this k-point and this spin polarization.
2580 
2581 !    Note that dfpt_nsteltwf is called with kpoint, while kpt is used inside dfpt_vtowfk
2582      call dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,&
2583 &     istwf_k,kg_k,kg1_k,kpoint,mband,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,&
2584 &     npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm_k,ylmgr_k)
2585      d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:)
2586      if(prtbbb==1)then
2587        d2bbb(:,:,idir,ipert,:,:) = d2bbb(:,:,idir,ipert,:,:) + &
2588 &       d2bbb_k(:,:,:,:)
2589      end if
2590 
2591      ABI_DEALLOCATE(occ_k)
2592 
2593 !    Keep track of total number of bands (all k points so far, even for
2594 !    k points not treated by me)
2595      bdtot_index=bdtot_index+nband_k
2596      bd2tot_index=bd2tot_index+2*nband_k**2
2597 
2598 !    Shift array memory
2599      if (mkmem/=0) then
2600        icg=icg+npw_k*nspinor*nband_k
2601        ikg=ikg+npw_k
2602      end if
2603      if (mk1mem/=0) then
2604        icg1=icg1+npw1_k*nspinor*nband_k
2605        ikg1=ikg1+npw1_k
2606      end if
2607      ABI_DEALLOCATE(ylm_k)
2608      ABI_DEALLOCATE(ylm1_k)
2609      if (ipert==natom+3.or.ipert==natom+4)  then
2610        ABI_DEALLOCATE(ylmgr_k)
2611        ABI_DEALLOCATE(ylmgr1_k)
2612      end if
2613 
2614    end do ! End big k point loop
2615  end do ! End loop over spins
2616 
2617  if(xmpi_paral==1)then
2618    call xmpi_barrier(comm)
2619    call wrtout(std_out,' dfpt_nselt: loop on k-points and spins done in parallel','COLL')
2620  end if
2621 
2622 !Treat now varying occupation numbers
2623 !if(occopt>=3 .and. occopt <=8) then
2624 !SUPPRESSED metallic coding of vtorho
2625 
2626 !Treat fixed occupation numbers
2627 !else
2628 
2629 !Accumulation over parallel processed now carried out for all terms
2630 !in dfpt_nstdy.f
2631 
2632 !End of test on varying or fixed occupation numbers
2633 !end if
2634 
2635 !The imaginary part of d2nl will be must be set to zero here since
2636 !time-reversal symmetry will always be true for the strain peturbation.
2637 !The symmetry-reduced kpt set will leave a non-zero imaginary part.
2638 
2639  d2nl(2,:,natom+3:natom+4,idir,ipert)=zero
2640 
2641 !Symmetrize the non-local contributions,
2642 !as was needed for the stresses in a ground-state calculation
2643 
2644  if (nsym1>1) then
2645 !  Pack like symmetric-storage cartesian stress tensor
2646    ii=0
2647    do ipert1=natom+3,natom+4
2648      do idir1=1,3
2649        ii=ii+1
2650        restr(ii)=d2nl(1,idir1,ipert1,idir,ipert)
2651      end do
2652    end do
2653 !  Do the symmetrization using the ground state routine
2654    call stresssym(gprimd,nsym1,restr,symrc1)
2655 !  Unpack symmetrized stress tensor
2656    ii=0
2657    do ipert1=natom+3,natom+4
2658      do idir1=1,3
2659        ii=ii+1
2660        d2nl(1,idir1,ipert1,idir,ipert)=restr(ii)
2661      end do
2662    end do
2663  end if !nsym>1
2664 
2665 !----------------------------------------------------------------------------
2666 !Now, treat the local contribution
2667 
2668  ABI_ALLOCATE(vpsp1,(cplex*nfft))
2669  n3xccc=0
2670  if(psps%n1xccc/=0)n3xccc=nfft
2671  ABI_ALLOCATE(xccc3d1,(cplex*n3xccc))
2672  ABI_ALLOCATE(vxc1,(cplex*nfft,nspden))
2673  ABI_ALLOCATE(vhartr01,(nfft))
2674  xccc3d1(:)=zero
2675 
2676 !Double loop over strain perturbations
2677  do ipert1=natom+3,natom+4
2678    do idir1=1,3
2679      if(ipert1==natom+3) then
2680        istr1=idir1
2681      else
2682        istr1=idir1+3
2683      end if
2684 
2685 !    Get first-order local potential.
2686      call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfft,mpi_enreg,&
2687 &     psps%mqgrid_vl,natom,gs_hamk%nattyp,nfft,ngfft,ntypat,paral_kgb,ph1d,psps%qgrid_vl,&
2688 &     ucvol,psps%vlspl,vpsp1)
2689 
2690 !    Get first-order hartree potential.
2691      call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,natom,nfft,ngfft,&
2692 &     paral_kgb,rhog,rprimd,vhartr01)
2693 
2694 !    Get first-order exchange-correlation potential
2695      if(psps%n1xccc/=0)then
2696        call dfpt_mkcore(cplex,idir1,ipert1,natom,ntypat,n1,psps%n1xccc,&
2697 &       n2,n3,qphon,rprimd,typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
2698      end if ! psps%n1xccc/=0
2699 
2700      option=0
2701      call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,natom,nfft,ngfft,&
2702 &     dummy,dummy,nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor,rhor1,rprimd,&
2703 &     0,0,vxc1,xccc3d1)
2704 
2705 !    Combines density j2 with local potential j1
2706      do ispden=1,min(nspden,2)
2707        do ifft=1,cplex*nfft
2708          vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft)+vhartr01(ifft)
2709        end do
2710      end do
2711      call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol)
2712      write(std_out,*)
2713      d2lo(1,idir1,ipert1,idir,ipert)=dotr
2714      d2lo(2,idir1,ipert1,idir,ipert)=doti
2715    end do ! istr1
2716  end do ! ipert1
2717 
2718  call destroy_hamiltonian(gs_hamk)
2719 
2720  ABI_DEALLOCATE(vxc1)
2721  ABI_DEALLOCATE(xccc3d1)
2722  ABI_DEALLOCATE(vhartr01)
2723 
2724  ABI_DEALLOCATE(d2bbb_k)
2725  ABI_DEALLOCATE(d2nl_k)
2726  ABI_DEALLOCATE(kg_k)
2727  ABI_DEALLOCATE(kg1_k)
2728  ABI_DEALLOCATE(vpsp1)
2729 
2730 end subroutine dfpt_nselt

ABINIT/dfpt_nstdy [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nstdy

FUNCTION

 This routine compute the non-stationary expression for the
 second derivative of the total energy, for a whole row of
 mixed derivatives.
 Only for norm-conserving pseudopotentials (no PAW)

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
  eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  idir=direction of the perturbation
  indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  ipert=type of the perturbation
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ
  kxc(nfft,nkxc)=exchange and correlation kernel
  mkmem =number of k points treated by this node (GS data)
  mk1mem =number of k points treated by this node (RF data)
  mpert =maximum number of ipert
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  nattyp(ntypat)= # atoms of each type.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
  nfft=(effective) number of FFT grid points (for this proc)
  ngfft(1:18)=integer array with FFT box dimensions and other
  nkpt=number of k points in the full BZ
  nkpt_rbz=number of k points in the reduced BZ for this perturbation
  nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym1=number of symmetry elements in space group consistent with i perturbation
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band
   and k in the reduced Brillouin zone (usually =2)
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3.
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  ucvol=unit cell volume in bohr**3.
  wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point

OUTPUT

  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
                                        second order derivatives
  d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
  d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs

NOTES

 Note that the ddk perturbation should not be treated here.

PARENTS

      dfpt_scfcv

CHILDREN

      appdig,destroy_hamiltonian,dfpt_mkcore,dfpt_mkvxc,dfpt_mkvxc_noncoll
      dfpt_nstwf,dfpt_sygra,dfpt_vlocal,dotprod_vn,init_hamiltonian
      load_spin_hamiltonian,mati3inv,timab,wfk_close,wfk_open_read,wrtout
      xmpi_sum

SOURCE

3065 subroutine dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,&
3066 &          gmet,gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,&
3067 &          mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,nfft,ngfft,nkpt,nkpt_rbz,nkxc,&
3068 &          npwarr,npwar1,nspden,nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,&
3069 &          symrc1,ucvol,wtk_rbz,xred,ylm,ylm1,rhor,vxc)
3070 
3071  use m_hamiltonian
3072 
3073  use m_symtk,     only : mati3inv
3074  use m_dynmat,    only : dfpt_sygra
3075  use m_dfpt_mkvxc,     only : dfpt_mkvxc
3076  use m_mkcore,         only : dfpt_mkcore
3077 
3078 !This section has been created automatically by the script Abilint (TD).
3079 !Do not modify the following lines by hand.
3080 #undef ABI_FUNC
3081 #define ABI_FUNC 'dfpt_nstdy'
3082 !End of the abilint section
3083 
3084  implicit none
3085 
3086 !Arguments -------------------------------
3087 !scalars
3088  integer,intent(in) :: cplex,idir,ipert,mk1mem,mkmem,mpert,mpw,mpw1,nfft,nkpt,nkpt_rbz,nkxc,nspden,nsppol,nsym1
3089  real(dp),intent(in) :: gsqcut,ucvol
3090  type(MPI_type),intent(in) :: mpi_enreg
3091  type(datafiles_type),intent(in) :: dtfil
3092  type(dataset_type),intent(in) :: dtset
3093  type(pseudopotential_type),intent(in) :: psps
3094 !arrays
3095  integer,intent(in) :: atindx(dtset%natom),indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom)
3096  integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
3097  integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol),ngfft(18)
3098  integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz),symrc1(3,3,nsym1)
3099  integer,intent(inout) :: blkflg(3,mpert,3,mpert) !vz_i
3100  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
3101  real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*nsppol)
3102  real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*nsppol)
3103  real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol)
3104  real(dp),intent(in) :: gmet(3,3),kpt_rbz(3,nkpt_rbz)
3105  real(dp),intent(in) :: kxc(nfft,nkxc),occ_rbz(dtset%mband*nkpt_rbz*nsppol)
3106  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
3107  real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3)
3108  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom)
3109  real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
3110  real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
3111  real(dp),intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*dtset%prtbbb)!vz_i
3112  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i
3113 ! optional
3114  real(dp),optional,intent(in) :: rhor(nfft,nspden)
3115  real(dp),optional,intent(in) :: vxc(cplex*nfft,nspden)
3116 
3117 !Local variables-------------------------------
3118 !scalars
3119  integer,parameter :: formeig1=1
3120  integer :: ban2tot,bantot,bdtot_index,ddkcase,iband,icg,icg1,idir1
3121  integer :: ierr,ifft,ii,ikg,ikg1,ikpt,ilm,ipert1,ispden,isppol
3122  integer :: istwf_k,isym,jj,master,me,n1,n2,n3,n3xccc,n4,n5,n6
3123  integer :: nband_k,nfftot,npw1_k,npw_k,nspinor_,option,spaceworld,optnc
3124  real(dp) :: doti,dotr,wtk_k
3125  logical :: t_exist
3126  character(len=500) :: msg
3127  character(len=fnlen) :: fiwfddk
3128  type(gs_hamiltonian_type) :: gs_hamkq
3129 !arrays
3130  integer :: ddkfil(3)
3131  integer,allocatable :: kg1_k(:,:),kg_k(:,:),symrl1(:,:,:)
3132  real(dp) :: d2nl_elfd(2,3),d2nl_mgfd(2,3),kpoint(3),kpq(3),sumelfd(2),summgfd(2),tsec(2)
3133  real(dp),allocatable :: buffer1(:),buffer2(:),d2bbb_k(:,:,:,:),d2nl_k(:,:,:)
3134  real(dp),allocatable :: eig1_k(:),eig_k(:),occ_k(:)
3135  real(dp) :: rhodummy(0,0)
3136  real(dp),allocatable :: vpsp1(:),vxc1(:,:),work1(:,:,:),xccc3d1(:),ylm1_k(:,:),ylm_k(:,:)
3137  type(pawtab_type) :: pawtab(dtset%ntypat*psps%usepaw)
3138  type(wfk_t) :: ddks(3)
3139 
3140 ! *********************************************************************
3141 
3142  ABI_UNUSED(nkpt)
3143 
3144  DBG_ENTER("COLL")
3145 
3146 !Not valid for PAW
3147  if (psps%usepaw==1) then
3148    msg='This routine cannot be used for PAW (use pawnst3 instead) !'
3149    MSG_BUG(msg)
3150  end if
3151 
3152 !Keep track of total time spent in dfpt_nstdy
3153  call timab(101,1,tsec)
3154 
3155 !Init parallelism
3156  spaceworld=mpi_enreg%comm_cell
3157  me=mpi_enreg%me_kpt
3158 
3159  master =0
3160 
3161 !Zero only portion of nonlocal matrix to be computed here
3162  d2nl(:,:,1:dtset%natom+2,idir,ipert)=zero
3163 
3164  ABI_ALLOCATE(d2bbb_k,(2,3,dtset%mband,dtset%mband*dtset%prtbbb))
3165  ABI_ALLOCATE(d2nl_k,(2,3,mpert))
3166  ABI_ALLOCATE(eig_k,(nsppol*dtset%mband))
3167  ABI_ALLOCATE(eig1_k,(2*nsppol*dtset%mband**2))
3168  ABI_ALLOCATE(kg_k,(3,mpw))
3169  ABI_ALLOCATE(kg1_k,(3,mpw1))
3170 
3171 !Do not try to open electric field file
3172  ddkfil(:)=0
3173 !The treatment of homogeneous electric field potential need the existence of d/dk files.
3174  do idir1=1,3
3175    ddkcase=idir1+dtset%natom*3
3176    call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk)
3177 
3178 !  Check that ddk file exists
3179    t_exist = file_exists(fiwfddk)
3180    if (.not. t_exist) then
3181      ! Try netcdf file.
3182      t_exist = file_exists(nctk_ncify(fiwfddk))
3183      if (t_exist) then
3184        fiwfddk = nctk_ncify(fiwfddk)
3185        write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name."
3186        call wrtout(std_out,msg,'COLL')
3187      end if
3188    end if
3189 
3190    if (t_exist) then
3191 !    Note the use of unit numbers 21, 22 and 23
3192      ddkfil(idir1)=20+idir1
3193      write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk)
3194      call wrtout(std_out,msg,'COLL')
3195      call wrtout(ab_out,msg,'COLL')
3196      call wfk_open_read(ddks(idir1),fiwfddk,formeig1,dtset%iomode,ddkfil(idir1),spaceworld)
3197    end if
3198  end do
3199 
3200 !Update list of computed matrix elements
3201  if (ipert /= dtset%natom + 1) then
3202    do ipert1=1,mpert
3203      do idir1=1,3
3204        if(ipert1 <= dtset%natom .or. ipert1==dtset%natom+2 &
3205 &       .and. ddkfil(idir1)/=0) then
3206          blkflg(idir1,ipert1,idir,ipert)=1
3207        end if
3208      end do
3209    end do
3210  else
3211    ipert1 = dtset%natom + 1
3212    do idir1=1,3
3213 !    If was already computed in another run or dataset, or if is to be computed in the present one
3214      if ((ddkfil(idir1) /= 0).or. (dtset%rfdir(idir1)/=0.and. idir1<=idir) ) then
3215 !      if ((ddkfil(idir1) /= 0).or. (idir1==idir) ) then
3216        blkflg(idir1,ipert1,idir,ipert)=1
3217      end if
3218    end do
3219  end if
3220 
3221  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
3222  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
3223  nspinor_=dtset%nspinor
3224 
3225  bantot = 0
3226  ban2tot = 0
3227 
3228 !==== Initialize most of the Hamiltonian ====
3229 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
3230 !2) Perform the setup needed for the non-local factors:
3231 !3) Constant kleimann-Bylander energies are copied from psps to gs_hamk.
3232  call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,dtset%natom,&
3233 & dtset%typat,xred,nfft,dtset%mgfft,ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
3234 & use_gpu_cuda=dtset%use_gpu_cuda)
3235 
3236 !LOOP OVER SPINS
3237  bdtot_index=0
3238  icg=0;icg1=0
3239  do isppol=1,nsppol
3240 
3241    ikg=0;ikg1=0
3242 
3243 !  Continue to initialize the Hamiltonian
3244    call load_spin_hamiltonian(gs_hamkq,isppol,with_nonlocal=.true.)
3245 
3246 !  BIG FAT k POINT LOOP
3247    do ikpt=1,nkpt_rbz
3248 
3249      nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
3250      istwf_k=istwfk_rbz(ikpt)
3251      npw_k=npwarr(ikpt)
3252      npw1_k=npwar1(ikpt)
3253 
3254      eig_k(1:nband_k) = eigen0(1+bantot:nband_k+bantot)
3255      eig1_k(1:2*nband_k**2) = eigen1(1+ban2tot:2*nband_k**2+ban2tot)
3256      bantot = bantot + nband_k
3257      ban2tot = ban2tot + 2*nband_k**2
3258 
3259      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
3260        bdtot_index=bdtot_index+nband_k
3261 !      The wavefunction blocks for ddk file is skipped elsewhere in the loop
3262 !      Skip the rest of the k-point loop
3263        cycle
3264      end if
3265 
3266      ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
3267      ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
3268 
3269 !    In case of electric field pert1, read ddk wfs file
3270 !    Note that the symmetries are not used for ddk, so read each k point
3271 !    Also take into account implicitely the parallelism over k points
3272 
3273      do idir1=1,3
3274        if (ddkfil(idir1)/=0) then
3275          ii = wfk_findk(ddks(idir1), kpt_rbz(:, ikpt))
3276          ABI_CHECK(ii == indkpt1(ikpt),  "ii !=  indkpt1")
3277        end if
3278      end do
3279 
3280      ABI_ALLOCATE(occ_k,(nband_k))
3281      occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
3282      kpoint(:)=kpt_rbz(:,ikpt)
3283      kpq(:)=kpoint(:)+dtset%qptn(:)
3284      wtk_k=wtk_rbz(ikpt)
3285      d2nl_k(:,:,:)=zero
3286      if(dtset%prtbbb==1)d2bbb_k(:,:,:,:)=zero
3287 
3288 !    Get plane-wave vectors and related data at k
3289      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
3290      if (psps%useylm==1) then
3291        do ilm=1,psps%mpsang*psps%mpsang
3292          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
3293        end do
3294      end if
3295 
3296 !    Get plane-wave vectors and related data at k+q
3297      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
3298      if (psps%useylm==1) then
3299        do ilm=1,psps%mpsang*psps%mpsang
3300          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
3301        end do
3302      end if
3303 
3304 !    Compute the eigenvalues, wavefunction,
3305 !    contributions to kinetic energy, nonlocal energy, forces,
3306 !    and update of rhor1 to this k-point and this spin polarization.
3307 !    Note that dfpt_nstwf is called with kpoint, while kpt is used inside dfpt_vtowfk
3308      call dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,&
3309 &     icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpoint,kpq,mkmem,mk1mem,mpert,&
3310 &     mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,&
3311 &     occ_k,psps,rmet,ddks,wtk_k,ylm_k,ylm1_k)
3312 
3313      d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:)
3314      if(dtset%prtbbb==1)d2bbb(:,:,idir,ipert,:,:)=d2bbb(:,:,idir,ipert,:,:)+d2bbb_k(:,:,:,:)
3315 
3316 !    Keep track of total number of bands
3317      bdtot_index=bdtot_index+nband_k
3318 
3319 !    Shift arrays memory
3320      if (mkmem/=0) then
3321        icg=icg+npw_k*dtset%nspinor*nband_k
3322        ikg=ikg+npw_k
3323      end if
3324      if (mk1mem/=0) then
3325        icg1=icg1+npw1_k*dtset%nspinor*nband_k
3326        ikg1=ikg1+npw1_k
3327      end if
3328 
3329      ABI_DEALLOCATE(occ_k)
3330      ABI_DEALLOCATE(ylm_k)
3331      ABI_DEALLOCATE(ylm1_k)
3332 
3333 !    End big k point loop
3334    end do
3335 
3336 !  End loop over spins
3337  end do
3338 
3339  call destroy_hamiltonian(gs_hamkq)
3340 
3341 !Treat fixed occupation numbers (as in vtorho)
3342  if(xmpi_paral==1)then
3343    ABI_ALLOCATE(buffer1,(2*3*mpert))
3344    ABI_ALLOCATE(buffer2,(2*3*mpert))
3345 !  Pack d2nl
3346    buffer1(1:2*3*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/2*3*mpert/))
3347 !  Build sum of everything
3348    call timab(48,1,tsec)
3349    call xmpi_sum(buffer1,buffer2,2*3*mpert,spaceworld,ierr)
3350    call timab(48,2,tsec)
3351 !  Unpack the final result
3352    d2nl(:,:,:,idir,ipert)=reshape(buffer2(:),(/2,3,mpert/))
3353    ABI_DEALLOCATE(buffer1)
3354    ABI_DEALLOCATE(buffer2)
3355    if(dtset%prtbbb==1)then
3356      ABI_ALLOCATE(buffer1,(2*3*dtset%mband*dtset%mband))
3357      ABI_ALLOCATE(buffer2,(2*3*dtset%mband*dtset%mband))
3358 !    Pack d2bbb
3359      buffer1(1:2*3*dtset%mband*dtset%mband)=reshape(d2bbb(:,:,idir,ipert,:,:),(/2*3*dtset%mband*dtset%mband/))
3360 !    Build sum of everything
3361      call timab(48,1,tsec)
3362      call xmpi_sum(buffer1,buffer2,2*3*dtset%mband*dtset%mband,spaceworld,ierr)
3363      call timab(48,2,tsec)
3364 !    Unpack the final result
3365      d2bbb(:,:,idir,ipert,:,:)=reshape(buffer2(:),(/2,3,dtset%mband,dtset%mband/))
3366      ABI_DEALLOCATE(buffer1)
3367      ABI_DEALLOCATE(buffer2)
3368    end if
3369  end if ! xmpi_paral==1
3370 
3371 !In the case of the strain perturbation time-reversal symmetry will always
3372 !be true so imaginary part of d2nl will be must be set to zero here since
3373 !the symmetry-reduced kpt set will leave a non-zero imaginary part.
3374  if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) d2nl(2,:,:,idir,ipert)=zero
3375 
3376 !In case of electric field ipert1, close the ddk wf files
3377  do idir1=1,3
3378    if (ddkfil(idir1)/=0)then
3379      call wfk_close(ddks(idir1))
3380    end if
3381  end do
3382 
3383 !Symmetrize the non-local contributions,
3384 !as was needed for the forces in a ground-state calculation
3385 !However, here the quantity is complex, and there are phases !
3386 
3387 !Do the transform
3388  ABI_ALLOCATE(work1,(2,3,dtset%natom))
3389  do ipert1=1,dtset%natom
3390    do idir1=1,3
3391      work1(1,idir1,ipert1)=d2nl(1,idir1,ipert1,idir,ipert)
3392      work1(2,idir1,ipert1)=d2nl(2,idir1,ipert1,idir,ipert)
3393    end do
3394  end do
3395  call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work1,indsy1,ipert,nsym1,dtset%qptn,symrc1)
3396  ABI_DEALLOCATE(work1)
3397 
3398 !Must also symmetrize the electric/magnetic field perturbation response !
3399 !(XG 000803 This was not implemented until now)
3400  if(sum(ddkfil(:))/=0)then
3401 !  Get the symmetry matrices in terms of real space basis
3402    ABI_ALLOCATE(symrl1,(3,3,nsym1))
3403    do isym=1,nsym1
3404      call mati3inv(symrc1(:,:,isym),symrl1(:,:,isym))
3405    end do
3406 !  There should not be any imaginary part, but stay general (for debugging)
3407    d2nl_elfd(:,:)=d2nl(:,:,dtset%natom+2,idir,ipert)
3408    do ii=1,3
3409      sumelfd(:)=zero
3410      summgfd(:)=zero
3411      do isym=1,nsym1
3412        do jj=1,3
3413          if(symrl1(ii,jj,isym)/=0)then
3414            if(ddkfil(jj)==0)then
3415              blkflg(ii,dtset%natom+2,idir,ipert)=0
3416            end if
3417          end if
3418        end do
3419        sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+&
3420 &       dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+&
3421 &       dble(symrl1(ii,3,isym))*d2nl_elfd(:,3)
3422        summgfd(:)=summgfd(:)+dble(symrl1(ii,1,isym))*d2nl_mgfd(:,1)+&
3423 &       dble(symrl1(ii,2,isym))*d2nl_mgfd(:,2)+&
3424 &       dble(symrl1(ii,3,isym))*d2nl_mgfd(:,3)
3425      end do
3426      d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1)
3427    end do
3428 
3429    if ((dtset%prtbbb==1).and.(ipert<=dtset%natom)) then
3430      do iband = 1,dtset%mband
3431        d2nl_elfd(:,:)=d2bbb(:,:,idir,ipert,iband,iband)
3432        do ii=1,3
3433          sumelfd(:)=zero
3434          do isym=1,nsym1
3435            sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+&
3436 &           dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+&
3437 &           dble(symrl1(ii,3,isym))*d2nl_elfd(:,3)
3438          end do
3439          d2bbb(:,ii,idir,ipert,iband,iband)=sumelfd(:)/dble(nsym1)
3440        end do
3441      end do  !iband
3442    end if
3443 
3444    ABI_DEALLOCATE(symrl1)
3445  end if
3446 
3447 !----------------------------------------------------------------------------
3448 !Now, treat the local contribution
3449 
3450  nfftot=ngfft(1)*ngfft(2)*ngfft(3)
3451  ABI_ALLOCATE(vpsp1,(cplex*nfft))
3452  if (ipert /= dtset%natom + 1) then
3453    n3xccc=0;if(psps%n1xccc/=0) n3xccc=nfft
3454    ABI_ALLOCATE(xccc3d1,(cplex*n3xccc))
3455    ABI_ALLOCATE(vxc1,(cplex*nfft,nspden))
3456 
3457    do ipert1=1,mpert
3458      do idir1=1,3
3459        if(ipert1 <= dtset%natom)then
3460 
3461 !        Get first-order local potential and first-order pseudo core density
3462          call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_ff,dtset%natom,&
3463 &         nattyp,nfft,ngfft,dtset%ntypat,n1,n2,n3,dtset%paral_kgb,ph1d,psps%qgrid_ff,&
3464 &         dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
3465          if(psps%n1xccc/=0)then
3466            call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,n1,psps%n1xccc,&
3467 &           n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
3468          end if
3469 
3470 !        Get first-order exchange-correlation potential (core-correction contribution only !)
3471          if(psps%n1xccc/=0)then
3472            option=0
3473 !FR SPr EB non-collinear magnetism
3474            if (nspden==4.and.present(rhor).and.present(vxc)) then
3475              optnc=1
3476              call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,rhodummy,0,&
3477 &             nkxc,nspden,n3xccc,optnc,option,dtset%paral_kgb,dtset%qptn,rhor,rhor1,&
3478 &             rprimd,0,vxc,vxc1,xccc3d1)
3479            else
3480              call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,&
3481 &             nkxc,nspden,n3xccc,option,dtset%paral_kgb,dtset%qptn,rhodummy,&
3482 &             rprimd,0,vxc1,xccc3d1)
3483            end if
3484          else
3485            vxc1(:,:)=zero
3486          end if
3487 
3488 !        Norm-conserving pseudpopotential case:
3489 !        Combines density j2 with local potential j1 (vpsp1 and vxc1)
3490 !        XG030514 : this is a first possible coding, however, each dotprod contains
3491 !        a parallel section (reduction), so it is better to use only one dotprod ...
3492 !        call dotprod_vn(cplex,rhor1,dr_psp1,di_psp1,mpi_enreg,nfft,nfftot,1,2,vpsp1,ucvol)
3493 !        call dotprod_vn(cplex,rhor1,dr_xc1,di_xc1,mpi_enreg,nfft,nfftot,nspden,2,vxc1,ucvol)
3494 !        dotr=dr_psp1+dr_xc1;doti=di_psp1+di_xc1... but then, one needs to overload vxc1
3495          do ispden=1,min(nspden,2)
3496            do ifft=1,cplex*nfft
3497              vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft)
3498            end do
3499          end do
3500          call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol)
3501 
3502 !        MVeithen 021212 : in case ipert = 2, these lines compute the local part
3503 !        of the Born effective charges from phonon and electric
3504 !        field type perturbations, see eq. 43 of
3505 !        X. Gonze and C. Lee, PRB 55, 10355 (1997) [[cite:Gonze1997a]]
3506 !        The minus sign is due to the fact that the effective charges
3507 !        are minus the second derivatives of the energy
3508          if (ipert == dtset%natom+2) then
3509            d2lo(1,idir1,ipert1,idir,ipert)=-dotr
3510            d2lo(2,idir1,ipert1,idir,ipert)=-doti
3511          else
3512            d2lo(1,idir1,ipert1,idir,ipert)=dotr
3513            d2lo(2,idir1,ipert1,idir,ipert)=doti
3514          end if
3515 !        Endif ipert1<=natom
3516        end if
3517      end do
3518    end do
3519 
3520    ABI_DEALLOCATE(vxc1)
3521    ABI_DEALLOCATE(xccc3d1)
3522 
3523  end if ! ipert /= natom +1
3524 
3525  ABI_DEALLOCATE(d2bbb_k)
3526  ABI_DEALLOCATE(d2nl_k)
3527  ABI_DEALLOCATE(kg_k)
3528  ABI_DEALLOCATE(kg1_k)
3529  ABI_DEALLOCATE(vpsp1)
3530  ABI_DEALLOCATE(eig_k)
3531  ABI_DEALLOCATE(eig1_k)
3532 
3533  call timab(101,2,tsec)
3534 
3535  DBG_EXIT("COLL")
3536 
3537 end subroutine dfpt_nstdy

ABINIT/dfpt_nsteltwf [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nsteltwf

FUNCTION

 This routine computes the non-local and kinetic contribution to the
 2DTE matrix elements, in the non-stationary formulation

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
    wavefunctions at k,q.
  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)  (NOT NEEDED !)
  effmass_free=effective mass for electrons (1. in common case)
  gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
  icg=shift to be applied on the location of data in the array cg
  icg1=shift to be applied on the location of data in the array cg1
  ikpt=number of the k-point
  isppol=1 for unpolarized, 2 for spin-polarized
  istwf_k=flag controlling the storage of WFs
  kg_k(3,npw_k)=reduced planewave coordinates.
  kg1_k(3,npw1_k)=reduced planewave coordinates at k+q, with RF k points
  kpoint(3)=k-point in reduced coordinates
  mband=maximum number of bands
  mkmem =number of k points treated by this node.
  mk1mem =number of k points treated by this node (RF data).
  mpert =maximum number of ipert
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  natom=number of atoms in cell.
  nband_k=number of bands at this k point for that spin polarization
  npw_k=number of plane waves at this k point
  npw1_k=number of plane waves at this k+q point
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat=number of types of atoms in unit cell.
  occ_k(nband_k)=occupation number for each band (usually 2) for each k.
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rmet(3,3)=real space metric (bohr**2)
  wtk_k=weight assigned to the k point.
  ylm(npw_k,mpsang*mpsang)= real spherical harmonics for each G and k point
  ylmgr(npw_k,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point

OUTPUT

  d2nl_k(2,3,mpert)=non-local contributions to
   non-stationary 2DTE, for the present k point, and perturbation idir, ipert

PARENTS

      dfpt_nselt

CHILDREN

      dotprod_g,kpgstr,load_k_hamiltonian,mkffnl,mkkin,nonlop

SOURCE

2790 subroutine dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,&
2791 &  istwf_k,kg_k,kg1_k,kpoint,mband,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,&
2792 &  npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm,ylmgr)
2793 
2794  use m_wffile
2795 
2796  use m_hamiltonian, only : gs_hamiltonian_type,load_k_hamiltonian
2797  use m_mkffnl,      only : mkffnl
2798  use m_nonlop,      only : nonlop
2799 
2800 !This section has been created automatically by the script Abilint (TD).
2801 !Do not modify the following lines by hand.
2802 #undef ABI_FUNC
2803 #define ABI_FUNC 'dfpt_nsteltwf'
2804 !End of the abilint section
2805 
2806  implicit none
2807 
2808 !Arguments ------------------------------------
2809 !scalars
2810  integer,intent(in) :: icg,icg1,ikpt,isppol,istwf_k,mband,mk1mem,mkmem,mpert,mpw,mpw1,natom
2811  integer,intent(in) :: nspinor,nsppol,ntypat
2812  integer,intent(inout) :: nband_k,npw1_k,npw_k
2813  real(dp),intent(in) :: ecut,ecutsm,effmass_free,wtk_k
2814  type(MPI_type),intent(in) :: mpi_enreg
2815  type(pseudopotential_type),intent(in) :: psps
2816 !arrays
2817  integer,intent(in) :: kg1_k(3,npw1_k),kg_k(3,npw_k)
2818  real(dp),intent(in) :: kpoint(3)
2819  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
2820  real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)
2821  real(dp),intent(in) :: occ_k(nband_k),rmet(3,3)
2822  real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang)
2823  real(dp),intent(in) :: ylmgr(npw_k,3,psps%mpsang*psps%mpsang)
2824  real(dp),intent(inout) :: d2nl_k(2,3,mpert)
2825 
2826 !Local variables-------------------------------
2827 !scalars
2828  integer :: choice,cpopt,dimffnl,dimffnl2,iband
2829  integer :: ider,idir0,idir1,ilmn,ipert1,ipw,ipws,ispinor,istr1,itypat
2830  integer :: nkpg,nnlout,paw_opt,signs,tim_nonlop
2831  real(dp) :: doti,dotr
2832  type(gs_hamiltonian_type) :: gs_hamk
2833 !arrays
2834  real(dp) :: enlout(6),dum_svectout(1,1),dum(1),kpg_dum(0,0)
2835  real(dp),allocatable :: cwave0(:,:),cwavef(:,:),dkinpw(:),eig2_k(:)
2836  real(dp),allocatable :: ffnl(:,:,:,:),ffnl_ylm(:,:,:,:),ghc(:,:)
2837  real(dp),allocatable :: gvnl1(:,:),gvnlc(:,:),kinpw1(:),ph3d(:,:,:)
2838  type(pawcprj_type) :: cprj_dum(0,0)
2839 
2840 ! *********************************************************************
2841 
2842 
2843 !DEBUG
2844 !write(std_out,*)' dfpt_nstwf : enter '
2845 !ENDDEBUG
2846 
2847 !Init me
2848  ABI_ALLOCATE(ghc,(2,npw1_k*nspinor))
2849  ABI_ALLOCATE(gvnlc,(2,npw1_k*nspinor))
2850  ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor))
2851  ABI_ALLOCATE(eig2_k,(2*nsppol*mband**2))
2852  ABI_ALLOCATE(kinpw1,(npw1_k))
2853  ABI_ALLOCATE(dkinpw,(npw_k))
2854  nkpg=0
2855 
2856 !Compute nonlocal form factors ffnl at (k+G), for all atoms
2857  dimffnl=2
2858  ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
2859  if (psps%useylm==0) then
2860    ider=1;idir0=0
2861    call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,ider,idir0,&
2862 &   psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
2863 &   npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr)
2864  else
2865    ider=1;idir0=-7;dimffnl2=7
2866    ABI_ALLOCATE(ffnl_ylm,(npw_k,dimffnl2,psps%lmnmax,ntypat))
2867    call mkffnl(psps%dimekb,dimffnl2,psps%ekb,ffnl_ylm,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
2868 &   ider,idir0,psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
2869 &   nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr)
2870    do itypat=1,ntypat
2871      do ilmn=1,psps%lmnmax
2872        ffnl(:,1,ilmn,itypat)=ffnl_ylm(:,1,ilmn,itypat)
2873      end do
2874    end do
2875  end if
2876 
2877 !Compute kinetic contributions (1/2) (2 Pi)**2 (k+G)**2:
2878 ! call mkkin(ecut,ecutsm,effmass_free,gs_hamk%gmet,kg1_k,kinpw1,kpoint,npw1_k)
2879  call mkkin(ecut,ecutsm,effmass_free,gs_hamk%gmet,kg1_k,kinpw1,kpoint,npw1_k,0,0)
2880 
2881 !Load k/k+q-dependent part in the Hamiltonian datastructure
2882  ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
2883  call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,ffnl_k=ffnl,&
2884 & ph3d_k=ph3d,compute_ph3d=.true.)
2885 
2886  ABI_ALLOCATE(cwave0,(2,npw_k*nspinor))
2887  ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor))
2888 
2889 !Loop over bands
2890  do iband=1,nband_k
2891 
2892    if(mpi_enreg%proc_distrb(ikpt, iband, isppol) /= mpi_enreg%me_kpt) then
2893 !    Skip the eigenvalue and the gvnl records of this band
2894      cycle
2895    end if
2896 
2897 !  Get ground-state and first-order wavefunctions
2898    cwave0(:,:)=cg(:,1+(iband-1)*npw_k*nspinor+icg:iband*npw_k*nspinor+icg)
2899    cwavef(:,:)=cg1(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)
2900 
2901 !  Double loop over strain perturbations
2902    do ipert1=natom+3,natom+4
2903      do idir1=1,3
2904        if (ipert1==natom+3) istr1=idir1
2905        if (ipert1==natom+4) istr1=idir1+3
2906 
2907 !      Compute the derivative of the kinetic operator vs strain in dkinpw
2908        call kpgstr(dkinpw,ecut,ecutsm,effmass_free,gs_hamk%gmet,gs_hamk%gprimd,istr1,&
2909 &       kg1_k,kpoint,npw1_k)
2910 
2911 !      Get |vnon-locj1|u0> :
2912 !      first-order non-local, applied to zero-order wavefunction
2913 !      (??) this routine gives MINUS the non-local contribution
2914 
2915 !      When using Ylms, load the correct ffnl derivative
2916        if (psps%useylm==1) then
2917          do itypat=1,ntypat
2918            do ilmn=1,psps%lmnmax
2919              ffnl(:,2,ilmn,itypat)=ffnl_ylm(:,1+istr1,ilmn,itypat)
2920            end do
2921          end do
2922        end if
2923 
2924        signs=2 ; choice=3 ; nnlout=6 ; paw_opt=0 ; cpopt=-1 ; tim_nonlop=5
2925        call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,istr1,dum,mpi_enreg,1,nnlout,paw_opt,&
2926 &       signs,dum_svectout,tim_nonlop,cwave0,gvnl1)
2927 !      <G|Vnl1|Cnk> is contained in gvnl1
2928 
2929 !      Kinetic contribution
2930        do ispinor=1,nspinor
2931          do ipw=1,npw1_k
2932            ipws=ipw+npw1_k*(ispinor-1)
2933            if(kinpw1(ipw)<huge(0.0_dp)*1.d-11)then
2934              gvnl1(1,ipws)=gvnl1(1,ipws)+dkinpw(ipw)*cwave0(1,ipws)
2935              gvnl1(2,ipws)=gvnl1(2,ipws)+dkinpw(ipw)*cwave0(2,ipws)
2936            else
2937              gvnl1(1,ipws)=0.0_dp
2938              gvnl1(2,ipws)=0.0_dp
2939            end if
2940          end do
2941        end do
2942 
2943 !      construct the matrix element (<uj2|vj1|u0>)complex conjug.
2944 !      and add it to the 2nd-order matrix
2945 !      imaginary term should be zero for strain-strain 2nd derivatives,
2946 !      but keep it as a test for now
2947        call dotprod_g(dotr,doti,gs_hamk%istwf_k,npw1_k*nspinor,2,cwavef,gvnl1,&
2948 &       mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2949 
2950        d2nl_k(1,idir1,ipert1)= d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*2.0_dp*dotr
2951        d2nl_k(2,idir1,ipert1)= d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*2.0_dp*doti
2952 
2953      end do !idir1
2954    end do !ipert1
2955 
2956 !  UNTIL NOW, DO NOT TAKE INTO ACCOUNT istwf_k
2957 
2958 !  End loop over bands
2959  end do
2960 
2961  ABI_DEALLOCATE(cwave0)
2962  ABI_DEALLOCATE(cwavef)
2963 
2964 !###################################################################
2965 
2966  ABI_DEALLOCATE(eig2_k)
2967  ABI_DEALLOCATE(ghc)
2968  ABI_DEALLOCATE(gvnlc)
2969  ABI_DEALLOCATE(gvnl1)
2970  ABI_DEALLOCATE(kinpw1)
2971  ABI_DEALLOCATE(dkinpw)
2972  ABI_DEALLOCATE(ffnl)
2973  ABI_DEALLOCATE(ph3d)
2974  if (psps%useylm==1)  then
2975    ABI_DEALLOCATE(ffnl_ylm)
2976  end if
2977 
2978 end subroutine dfpt_nsteltwf

ABINIT/dfpt_rhofermi [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_rhofermi

FUNCTION

 This routine computes the fixed contribution to the first-order
 Fermi energy for metallic occupation and Q=0, as well as the
 Fermi level charge density needed to compute the remainder of the
 first-order Fermi energy from the self-consistent local potential
 at each step in the iteration process.

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions.
  cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL; if 2, COMPLEX
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
              projected with non-local projectors: cprj=<p_i|Cnk>
  cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q
              projected with non-local projectors: cprjq=<p_i|Cnk+q>
  doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
  docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
  idir=direction of the perturbation
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  ipert=type of the perturbation
  irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
  mband=maximum number of bands
  mkmem =number of k points treated by this node (GS data)
  mkqmem =number of k+q points treatede by this node (GS data)
  mk1mem =number of k points treated by this node.
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  maximum dimension for q points in grids for nonlocal form factors
  natom=number of atoms in cell.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  nfftf= -PAW ONLY- number of FFT grid points for the fine grid
         (nfftf=nfft for norm-conserving potential runs - see comment in respfn.F90)
  ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid
  nhatfermi(nfft,nspden)=array for fermi-level compensation charge density (PAW only)
  nkpt_rbz=number of k points in the IBZ for this perturbation
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym1=number of symmetry elements in space group consistent with
    perturbation
  occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
   at each k+q point of the reduced Brillouin zone.
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k (usually 2)
  paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawang1 <type(pawang_type)>=pawang datastr. containing only symmetries preserving the perturbation
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  prtvol=control print volume and debugging output
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional real space primitive translations
  symaf1(nsym1)=(anti)ferromagnetic part of symmetry operations
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  symrl1(3,3,nsym1)=3x3 matrices of the group symmetries
  ucvol=volume of the unit cell
  usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
  useylmgr1= 1 if ylmgr1 array is allocated
  vtrial(nfftf,nspden)=GS potential (Hartree).
  vxc(nfftf,nspden)=XC potential (Hartree).
  wtk_rbz(nkpt_rbz)=weight assigned to each k point.
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= spherical harmonics for each G and k+g point
  ylmgr1(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q

OUTPUT

  eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues
   (hartree) - only digonal elements computed here
  fe1fixed=fixed contribution to the first-order Fermi energy
   (nonlocal and kinetic in the case of strain)
  nhatfermi(cplex*nfftf,nspden)=fermi-level compensation charge density (PAW only)
  rhorfermi(cplex*nfftf,nspden)=fermi-level electronic density

NOTES

  This routine will NOT work with nspden==4:
    at least the use of fftpac should be modified.

PARENTS

      dfpt_scfcv

CHILDREN

      destroy_hamiltonian,destroy_rf_hamiltonian,dfpt_wfkfermi,fftpac
      init_hamiltonian,init_rf_hamiltonian,kpgstr,load_k_hamiltonian
      load_k_rf_hamiltonian,load_kprime_hamiltonian,load_spin_hamiltonian
      load_spin_rf_hamiltonian,mkffnl,mkkin,mkkpg,occeig,paw_ij_free
      paw_ij_init,paw_ij_nullify,pawdijfr,pawmkrho,pawrhoij_alloc
      pawrhoij_free,pawrhoij_free_unpacked,pawrhoij_init_unpacked
      pawrhoij_mpisum_unpacked,status,symrhg,timab,xmpi_sum

SOURCE

3656 subroutine dfpt_rhofermi(cg,cgq,cplex,cprj,cprjq,&
3657 & doccde_rbz,docckqde,dtfil,dtset,eigenq,eigen0,eigen1,fe1fixed,gmet,gprimd,idir,&
3658 & indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,mkmem,mkqmem,mk1mem,mpi_enreg,&
3659 & mpw,mpw1,my_natom,natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1,nspden,&
3660 & nsppol,nsym1,occkq,occ_rbz,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,&
3661 & phnons1,ph1d,prtvol,psps,rhorfermi,rmet,rprimd,symaf1,symrc1,symrl1,&
3662 & ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1)
3663 
3664 
3665  use m_hamiltonian
3666 
3667  use m_occ,         only : occeig
3668  use m_pawrhoij,    only : pawrhoij_type, pawrhoij_init_unpacked, pawrhoij_gather, &
3669 &                          pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, &
3670 &                          pawrhoij_free_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_get_nspden
3671  use m_paw_mkrho,   only : pawmkrho
3672 
3673  use m_spacepar,    only : symrhg
3674  use m_mkffnl,      only : mkffnl
3675 
3676 !This section has been created automatically by the script Abilint (TD).
3677 !Do not modify the following lines by hand.
3678 #undef ABI_FUNC
3679 #define ABI_FUNC 'dfpt_rhofermi'
3680 !End of the abilint section
3681 
3682  implicit none
3683 
3684 !Arguments -------------------------------
3685 !scalars
3686  integer,intent(in) :: cplex,idir,ipert,mband,mk1mem,mkmem,mkqmem
3687  integer,intent(in) :: mpw,mpw1,my_natom,natom,ncpgr,nfftf,nkpt_rbz,nspden,nsppol,nsym1
3688  integer,intent(in) :: prtvol,usecprj,useylmgr1
3689  real(dp),intent(in) :: ucvol
3690  real(dp),intent(out) :: fe1fixed
3691  type(MPI_type),intent(in) :: mpi_enreg
3692  type(datafiles_type),intent(in) :: dtfil
3693  type(dataset_type),intent(in) :: dtset
3694  type(pawang_type),intent(in) :: pawang,pawang1
3695  type(pawfgr_type),intent(in) :: pawfgr
3696  type(pseudopotential_type),intent(in) :: psps
3697 !arrays
3698  integer,intent(in) :: indsy1(4,nsym1,natom)
3699  integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))
3700  integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
3701  integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfftf(18)
3702  integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz),symaf1(nsym1)
3703  integer,intent(in) :: symrc1(3,3,nsym1),symrl1(3,3,nsym1)
3704  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol)
3705  real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband*mkqmem*nsppol)
3706  real(dp),intent(in) :: doccde_rbz(mband*nkpt_rbz*nsppol)
3707  real(dp),intent(in) :: docckqde(mband*nkpt_rbz*nsppol)
3708  real(dp),intent(in) :: eigen0(mband*nkpt_rbz*nsppol)
3709  real(dp),intent(in) :: eigenq(mband*nkpt_rbz*nsppol),gmet(3,3),gprimd(3,3)
3710  real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz)
3711  real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol),occkq(mband*nkpt_rbz*nsppol)
3712  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*natom)
3713  real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))
3714  real(dp),intent(in) :: rmet(3,3),rprimd(3,3),vtrial(nfftf,nspden),vxc(nfftf,nspden),wtk_rbz(nkpt_rbz)
3715  real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
3716  real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
3717  real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
3718  real(dp),intent(out) :: eigen1(2*mband*mband*nkpt_rbz*nsppol)
3719  real(dp),intent(out) :: nhatfermi(:,:)
3720  real(dp),intent(out) :: rhorfermi(cplex*nfftf,nspden)
3721  type(pawcprj_type),intent(in) :: cprj (natom,dtset%nspinor*mband*mkmem *nsppol*usecprj)
3722  type(pawcprj_type),intent(in) :: cprjq(natom,dtset%nspinor*mband*mkqmem*nsppol*usecprj)
3723  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
3724  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
3725  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
3726  type(pawrhoij_type),target,intent(inout)::pawrhoijfermi(my_natom*psps%usepaw)!vz_i
3727  type(pawtab_type), intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
3728 
3729 !Local variables-------------------------------
3730 !scalars
3731  integer,parameter :: level=17
3732  integer :: bd2tot_index,bdtot_index,buffer_size,counter,cplex_rhoij
3733  integer :: dimffnl1,dimffnlk,iatom,iband,ibg,ibgq
3734  integer :: icg,icgq,ider,idir0,ierr,iexit,ii,ikg,ikg1,ikpt,ilm,ilmn,indx
3735  integer :: ispden,isppol,istr,istwf_k
3736  integer :: mbd2kpsp,mcgq,mcgq_disk,mcprjq,mcprjq_disk
3737  integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nkpg1,npw1_k,npw_k,nspden_rhoij
3738  integer :: optfr,spaceworld
3739  logical :: paral_atom,qne0
3740  real(dp) :: arg,fe1norm,invfe1norm,wtk_k
3741  type(gs_hamiltonian_type) :: gs_hamkq
3742  type(rf_hamiltonian_type) :: rf_hamkq
3743 !arrays
3744  integer,allocatable :: kg1_k(:,:),kg_k(:,:)
3745  real(dp) :: kpoint(3),kpq(3),tsec(2)
3746  real(dp) :: ylmgr_dum(1,1,1)
3747  real(dp),allocatable :: buffer1(:),dkinpw(:),doccde_k(:)
3748  real(dp),allocatable :: doccde_kq(:),eig0_k(:),eig0_kq(:),eig1_k(:)
3749  real(dp),allocatable :: fe1fixed_k(:),fe1norm_k(:)
3750  real(dp),allocatable :: ffnl1(:,:,:,:),ffnlk(:,:,:,:)
3751  real(dp),allocatable :: kinpw1(:),kpg1_k(:,:),kpg_k(:,:)
3752  real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:)
3753  real(dp),allocatable :: rhoaug(:,:,:),rhogfermi(:,:),rhowfr(:,:)
3754  real(dp),allocatable :: rhoaug4(:,:,:,:)
3755  real(dp),allocatable :: rocceig(:,:),ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:)
3756  type(paw_ij_type),allocatable :: paw_ij1fr(:)
3757  type(pawrhoij_type),pointer :: pawrhoijfermi_unsym(:)
3758 ! real(dp),allocatable :: vlocal1(:,:,:,:),vlocal_tmp(:,:,:,:)
3759 ! real(dp),allocatable :: v1zeeman(:,:),vtrial_tmp(:,:)
3760 
3761 
3762 ! *********************************************************************
3763 
3764  DBG_ENTER('COLL')
3765 
3766 !Check arguments validity
3767  if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4.and.ipert/=natom+5) then
3768    MSG_BUG('wrong ipert argument!')
3769  end if
3770  if (cplex/=1) then
3771    MSG_BUG('wrong cplex/=1 argument !')
3772  end if
3773 
3774 
3775 !Keep track of total time spent in this routine
3776  call timab(121,1,tsec)
3777  call timab(124,1,tsec)
3778 
3779 !Retrieve parallelism data
3780  spaceworld=mpi_enreg%comm_cell
3781  me=mpi_enreg%me_kpt
3782  paral_atom=(my_natom/=dtset%natom)
3783 
3784 !Initialize output variables
3785  fe1fixed=zero
3786  if (psps%usepaw==0) rhorfermi(:,:)=zero
3787 
3788 !Initialisations/allocation of temporary variables
3789  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
3790  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
3791  bdtot_index=0 ; bd2tot_index=0 ; ibg=0 ; ibgq=0 ; icg=0 ; icgq=0
3792  qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14)
3793  mbd2kpsp=2*mband**2*nkpt_rbz*nsppol
3794  fe1norm=zero
3795  if (nspden/=4) then
3796    ABI_ALLOCATE(rhoaug,(cplex*n4,n5,n6))
3797  else
3798    ABI_ALLOCATE(rhoaug4,(cplex*n4,n5,n6,nspden))
3799  end if
3800  ABI_ALLOCATE(kg_k,(3,mpw))
3801  ABI_ALLOCATE(kg1_k,(3,mpw1))
3802  if (psps%usepaw==1) then
3803    ABI_ALLOCATE(rhowfr,(cplex*dtset%nfft,dtset%nspden))
3804    rhowfr(:,:)=zero
3805  end if
3806 
3807  mcgq=mpw1*dtset%nspinor*mband*mkqmem*nsppol;mcgq_disk=0
3808 
3809 !Prepare RF PAW files for reading and writing if mkmem, mkqmem or mk1mem==0
3810  if (psps%usepaw==1) then
3811    mcprjq=dtset%nspinor*mband*mkqmem*nsppol*usecprj;mcprjq_disk=0
3812  else
3813    mcprjq=0;mcprjq_disk=0
3814  end if
3815 
3816 !PAW:has to compute frozen part of Dij^(1) (without Vpsp(1) contribution)
3817  if (psps%usepaw==1) then
3818    ABI_DATATYPE_ALLOCATE(paw_ij1fr,(my_natom))
3819    call paw_ij_nullify(paw_ij1fr)
3820    call paw_ij_init(paw_ij1fr,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,&
3821 &   dtset%natom,dtset%ntypat,dtset%typat,pawtab,has_dijfr=1,&
3822 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom )
3823    optfr=1
3824    ABI_ALLOCATE(buffer1,(0))
3825    call pawdijfr(cplex,gprimd,idir,ipert,my_natom,natom,nfftf,ngfftf,dtset%nspden,dtset%nsppol,&
3826 &   dtset%ntypat,optfr,paw_ij1fr,pawang,pawfgrtab,pawrad,pawtab,&
3827 &   dtset%qptn,rprimd,ucvol,buffer1,vtrial,vxc,xred,&
3828 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
3829    ABI_DEALLOCATE(buffer1)
3830  end if
3831 
3832 !PAW:allocate memory for non-symetrized occupancies matrix at EFermi (pawrhoijfermi)
3833  pawrhoijfermi_unsym => pawrhoijfermi
3834  if (psps%usepaw==1) then
3835    if (paral_atom) then
3836      ABI_DATATYPE_ALLOCATE(pawrhoijfermi_unsym,(natom))
3837      cplex_rhoij=max(cplex,dtset%pawcpxocc)
3838      nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
3839      call pawrhoij_alloc(pawrhoijfermi_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
3840 &     dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
3841    else
3842      call pawrhoij_init_unpacked(pawrhoijfermi_unsym)
3843    end if
3844  end if
3845 
3846 !Initialize most of the Hamiltonian (arrays and quantities that do not depend on k + nl form factors)
3847  call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,natom,&
3848 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
3849 & paw_ij=paw_ij,usecprj=usecprj,ph1d=ph1d,use_gpu_cuda=dtset%use_gpu_cuda,&
3850 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_spintab=mpi_enreg%my_isppoltab)
3851  call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,paw_ij1=paw_ij1fr,&
3852 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_spintab=mpi_enreg%my_isppoltab)
3853 
3854 
3855 
3856 !LOOP OVER SPINS
3857  do isppol=1,nsppol
3858    ikg=0;ikg1=0
3859 !  Continue to initialize the Hamiltonian at k+q
3860    call load_spin_hamiltonian(gs_hamkq,isppol,with_nonlocal=.true.)
3861 !  call load_spin_rf_hamiltonian(rf_hamkq,isppol,with_nonlocal=.true.)
3862 
3863    call load_spin_rf_hamiltonian(rf_hamkq,isppol,with_nonlocal=.true.)
3864 
3865 
3866 !  Nullify contribution to density at EFermi from this k-point
3867    if (nspden/=4) then
3868      rhoaug(:,:,:)=zero
3869    else
3870      rhoaug4(:,:,:,:)=zero
3871    end if
3872    call timab(125,1,tsec)
3873 
3874 !  BIG FAT k POINT LOOP
3875    do ikpt=1,nkpt_rbz
3876 
3877      counter=100*ikpt+isppol
3878      call status(counter,dtfil%filstat,iexit,level,'loop ikpt     ')
3879      nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
3880      istwf_k=istwfk_rbz(ikpt)
3881      npw_k=npwarr(ikpt)
3882      npw1_k=npwar1(ikpt)
3883      wtk_k=wtk_rbz(ikpt)
3884 
3885      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
3886        eigen1(1+bd2tot_index : 2*nband_k**2+bd2tot_index) = zero
3887        bdtot_index=bdtot_index+nband_k
3888        bd2tot_index=bd2tot_index+2*nband_k**2
3889 !      Skip the rest of the k-point loop
3890        cycle
3891      end if
3892 
3893      ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
3894      ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
3895      ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
3896 
3897 !    Continue to initialize the Hamiltonian at k+q
3898      kpoint(:)=kpt_rbz(:,ikpt)
3899      kpq(:)=kpoint(:)+dtset%qptn(1:3)
3900 
3901      ABI_ALLOCATE(doccde_k,(nband_k))
3902      ABI_ALLOCATE(doccde_kq,(nband_k))
3903      ABI_ALLOCATE(eig0_k,(nband_k))
3904      ABI_ALLOCATE(eig0_kq,(nband_k))
3905      ABI_ALLOCATE(eig1_k,(2*nband_k**2))
3906      ABI_ALLOCATE(fe1fixed_k,(nband_k))
3907      ABI_ALLOCATE(fe1norm_k,(nband_k))
3908      ABI_ALLOCATE(occ_k,(nband_k))
3909      ABI_ALLOCATE(occ_kq,(nband_k))
3910      ABI_ALLOCATE(rocceig,(nband_k,nband_k))
3911 
3912      eig1_k(:)=zero
3913      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
3914      eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
3915      occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
3916      occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index)
3917      doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index)
3918      doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index)
3919 
3920 !    For each pair of active bands (m,n), generates the ratios
3921 !    rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
3922 !    and decide to which band to attribute it.
3923      call occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,&
3924 &     dtset%occopt,occ_k,occ_kq,rocceig)
3925 
3926 !    Get plane-wave coeffs and related data at k
3927      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
3928      if (psps%useylm==1) then
3929        do ilm=1,psps%mpsang*psps%mpsang
3930          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
3931        end do
3932      end if
3933 
3934 !    Get plane-wave coeffs and related data at k+q
3935      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
3936      if (psps%useylm==1) then
3937        do ilm=1,psps%mpsang*psps%mpsang
3938          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
3939        end do
3940        if (useylmgr1==1) then
3941          do ilm=1,psps%mpsang*psps%mpsang
3942            do ii=1,3
3943              ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
3944            end do
3945          end do
3946        end if
3947      end if
3948 
3949 !    Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
3950 
3951 !    Compute (k+G) vectors
3952      nkpg=0;if(ipert>=1.and.ipert<=natom) nkpg=3*dtset%nloalg(3)
3953      ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
3954      if (nkpg>0) then
3955        call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
3956      end if
3957 
3958 !    Compute (k+q+G) vectors
3959      nkpg1=0;if(ipert>=1.and.ipert<=natom) nkpg1=3*dtset%nloalg(3)
3960      ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1))
3961      if (nkpg1>0) then
3962        call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k)
3963      end if
3964 
3965 !    ===== Preparation of non-local contributions
3966 
3967      dimffnlk=0;if (ipert<=natom) dimffnlk=1
3968      ABI_ALLOCATE(ffnlk,(npw_k,dimffnlk,psps%lmnmax,dtset%ntypat))
3969 
3970 !    Compute nonlocal form factors ffnlk at (k+G)
3971      if (ipert<=natom ) then
3972        ider=0;idir0=0
3973        call status(counter,dtfil%filstat,iexit,level,'call mkffnl(0)')
3974        call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,&
3975 &       gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
3976 &       psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,dtset%ntypat,&
3977 &       psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
3978      end if
3979 
3980 !    Compute nonlocal form factors ffnl1 at (k+q+G)
3981      !-- Atomic displacement perturbation
3982      if (ipert<=natom) then
3983        ider=0;idir0=0
3984      !-- Strain perturbation
3985      else if (ipert==natom+3.or.ipert==natom+4) then
3986        if (ipert==natom+3) istr=idir
3987        if (ipert==natom+4) istr=idir+3
3988        ider=1;idir0=-istr
3989      else if (ipert==natom+5) then !SPr deb rfmagn
3990        ider=0;idir0=0
3991      end if
3992      dimffnl1=1+ider;if (ider==1.and.idir0==0) dimffnl1=dimffnl1+2*psps%useylm
3993      ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,dtset%ntypat))
3994      call status(counter,dtfil%filstat,iexit,level,'call mkffnl(1)')
3995      call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,&
3996 &     psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
3997 &     npw1_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
3998 
3999 !    ===== Preparation of kinetic contributions
4000 
4001      ABI_ALLOCATE(dkinpw,(npw_k))
4002      ABI_ALLOCATE(kinpw1,(npw1_k))
4003 
4004 !    Compute the derivative of the kinetic operator vs strain in dkinpw
4005      if (ipert==natom+3.or.ipert==natom+4) then
4006        if (ipert==natom+3) istr=idir
4007        if (ipert==natom+4) istr=idir+3
4008        call status(counter,dtfil%filstat,iexit,level,'call kpgstr   ')
4009        call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr,&
4010 &       kg_k,kpoint,npw_k)
4011      end if
4012 
4013 !    Compute (1/2) (2 Pi)**2 (k+q+G)**2:
4014      call status(counter,dtfil%filstat,iexit,level,'call mkkin(1) ')
4015 !     call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k)
4016      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
4017 
4018 !    ===== Load the k/k+q dependent parts of the Hamiltonian
4019 
4020 !    Load k-dependent part in the Hamiltonian datastructure
4021      ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamkq%matblk))
4022      call load_k_hamiltonian(gs_hamkq,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,&
4023 &     ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
4024      if (size(ffnlk)>0) then
4025        call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnlk)
4026      else
4027        call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1)
4028      end if
4029 
4030 !    Load k+q-dependent part in the Hamiltonian datastructure
4031 !        Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
4032      call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
4033 &     kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,&
4034 &     compute_gbound=.true.)
4035      if (qne0) then
4036        ABI_ALLOCATE(ph3d1,(2,npw1_k,gs_hamkq%matblk))
4037        call load_kprime_hamiltonian(gs_hamkq,ph3d_kp=ph3d1,compute_ph3d=.true.)
4038      end if
4039 
4040 !    Load k-dependent part in the 1st-order Hamiltonian datastructure
4041      call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw)
4042 
4043 !    Compute fixed contributions to 1st-order Fermi energy
4044 !    and Fermi level charge density
4045      call status(counter,dtfil%filstat,iexit,level,'call dfpt_wfkfermi  ')
4046      fe1fixed_k(:)=zero ; fe1norm_k(:)=zero
4047 
4048 !    Note that dfpt_wfkfermi is called with kpoint, while kpt is used inside dfpt_wfkfermi
4049      if (nspden/=4) then
4050        call dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,dtfil,eig0_k,eig1_k,fe1fixed_k,&
4051 &       fe1norm_k,gs_hamkq,ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,dtset%kptopt,mband,&
4052 &       mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k,&
4053 &       pawrhoijfermi_unsym,prtvol,rf_hamkq,rhoaug,rocceig,wtk_k)
4054      else
4055        call dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,dtfil,eig0_k,eig1_k,fe1fixed_k,&
4056 &       fe1norm_k,gs_hamkq,ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,dtset%kptopt,mband,&
4057 &       mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k,&
4058 &       pawrhoijfermi_unsym,prtvol,rf_hamkq,rhoaug4,rocceig,wtk_k)
4059      end if
4060 !    Free temporary storage
4061      ABI_DEALLOCATE(kpg_k)
4062      ABI_DEALLOCATE(kpg1_k)
4063      ABI_DEALLOCATE(dkinpw)
4064      ABI_DEALLOCATE(ffnlk)
4065      ABI_DEALLOCATE(ffnl1)
4066      ABI_DEALLOCATE(kinpw1)
4067      ABI_DEALLOCATE(doccde_k)
4068      ABI_DEALLOCATE(doccde_kq)
4069      ABI_DEALLOCATE(eig0_k)
4070      ABI_DEALLOCATE(eig0_kq)
4071      ABI_DEALLOCATE(occ_kq)
4072      ABI_DEALLOCATE(rocceig)
4073      ABI_DEALLOCATE(ylm_k)
4074      ABI_DEALLOCATE(ylm1_k)
4075      ABI_DEALLOCATE(ylmgr1_k)
4076      ABI_DEALLOCATE(ph3d)
4077      if (allocated(ph3d1)) then
4078        ABI_DEALLOCATE(ph3d1)
4079      end if
4080      call status(counter,dtfil%filstat,iexit,level,'after dfpt_wfkfermi ')
4081 
4082 !    Save eigenvalues (hartree)
4083      eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:)
4084 
4085 !    Accumulate sum over k points for 1st-order Fermi energy components
4086      do iband=1,nband_k
4087        fe1fixed=fe1fixed+wtk_k*occ_k(iband)*fe1fixed_k(iband)
4088        fe1norm=fe1norm+wtk_k*occ_k(iband)*fe1norm_k(iband)
4089      end do
4090 
4091      ABI_DEALLOCATE(eig1_k)
4092      ABI_DEALLOCATE(occ_k)
4093      ABI_DEALLOCATE(fe1fixed_k)
4094      ABI_DEALLOCATE(fe1norm_k)
4095 
4096 !    Keep track of total number of bands
4097 !    (all k points so far, even for k points not treated by me)
4098      bdtot_index=bdtot_index+nband_k
4099      bd2tot_index=bd2tot_index+2*nband_k**2
4100 
4101 !    Shift array memory
4102      if (mkmem/=0) then
4103        ibg=ibg+nband_k
4104        icg=icg+npw_k*dtset%nspinor*nband_k
4105        ikg=ikg+npw_k
4106      end if
4107      if (mkqmem/=0) then
4108        ibgq=ibgq+dtset%nspinor*nband_k
4109        icgq=icgq+npw1_k*dtset%nspinor*nband_k
4110      end if
4111      if (mk1mem/=0) then
4112        ikg1=ikg1+npw1_k
4113      end if
4114 
4115 !    End big k point loop
4116    end do
4117 
4118    call timab(125,2,tsec)
4119 
4120 !  Transfer density on augmented fft grid to normal fft grid in real space
4121 !  Also take into account the spin.
4122    if (nspden/=4) then
4123      if (psps%usepaw==0) then
4124        call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhorfermi,rhoaug,1)
4125      else
4126        call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhowfr   ,rhoaug,1)
4127      end if
4128    else
4129      if (psps%usepaw==0) then
4130        do ispden=1,4
4131          call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhorfermi,rhoaug4(:,:,:,ispden),1)
4132        end do
4133      end if
4134    end if
4135 
4136  end do ! End loop over spins
4137 
4138  !if(xmpi_paral==1)then
4139  !  call timab(166,1,tsec)
4140  !  call wrtout(std_out,'dfpt_rhofermi: loop on k-points and spins done in parallel','COLL')
4141  !  call xmpi_barrier(spaceworld)
4142  !  call timab(166,2,tsec)
4143  !end if
4144 
4145 !More memory cleaning
4146  call destroy_hamiltonian(gs_hamkq)
4147  call destroy_rf_hamiltonian(rf_hamkq)
4148  if(psps%usepaw==1) then
4149    call paw_ij_free(paw_ij1fr)
4150    ABI_DATATYPE_DEALLOCATE(paw_ij1fr)
4151  end if
4152  if (nspden/=4) then
4153    ABI_DEALLOCATE(rhoaug)
4154  else
4155    ABI_DEALLOCATE(rhoaug4)
4156  end if
4157  ABI_DEALLOCATE(kg_k)
4158  ABI_DEALLOCATE(kg1_k)
4159 
4160  call status(0,dtfil%filstat,iexit,level,'after loops   ')
4161  call timab(124,2,tsec)
4162 
4163 
4164 !=== MPI communications ==================
4165  if(xmpi_paral==1)then
4166    call timab(129,1,tsec)
4167 
4168 !  Identify MPI buffer size
4169    buffer_size=cplex*dtset%nfft*nspden+2+mbd2kpsp
4170    ABI_ALLOCATE(buffer1,(buffer_size))
4171 
4172 !  Pack rhorfermi, fe1fixed, fe1norm
4173    indx=cplex*dtset%nfft*nspden
4174    if (psps%usepaw==0) then
4175      buffer1(1:indx)=reshape(rhorfermi,(/indx/))
4176    else
4177      buffer1(1:indx)=reshape(rhowfr,(/indx/))
4178    end if
4179    buffer1(indx+1)=fe1fixed ; buffer1(indx+2)=fe1norm
4180    indx=indx+2 ; bd2tot_index=0
4181    do isppol=1,nsppol
4182      do ikpt=1,nkpt_rbz
4183        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
4184        buffer1(indx+1:indx+2*nband_k**2)=eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2)
4185        bd2tot_index=bd2tot_index+2*nband_k**2
4186        indx=indx+2*nband_k**2
4187      end do
4188    end do
4189    if(indx<buffer_size)buffer1(indx+1:buffer_size)=zero
4190 
4191 !  Build sum of everything
4192    call timab(48,1,tsec)
4193    call xmpi_sum(buffer1,buffer_size,spaceworld,ierr)
4194    call timab(48,2,tsec)
4195 
4196 !  Unpack the final result
4197    indx=cplex*dtset%nfft*nspden
4198    if (psps%usepaw==0) then
4199      rhorfermi(:,:)=reshape(buffer1(1:indx),(/cplex*dtset%nfft,nspden/))
4200    else
4201      rhowfr(:,:)=reshape(buffer1(1:indx),(/cplex*dtset%nfft,nspden/))
4202    end if
4203    fe1fixed=buffer1(indx+1) ; fe1norm =buffer1(indx+2)
4204    indx=indx+2 ; bd2tot_index=0
4205    do isppol=1,nsppol
4206      do ikpt=1,nkpt_rbz
4207        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
4208        eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2)=buffer1(indx+1:indx+2*nband_k**2)
4209        bd2tot_index=bd2tot_index+2*nband_k**2
4210        indx=indx+2*nband_k**2
4211      end do
4212    end do
4213    ABI_DEALLOCATE(buffer1)
4214 
4215 !  Accumulate PAW occupancies
4216    if (psps%usepaw==1) then
4217      call pawrhoij_mpisum_unpacked(pawrhoijfermi_unsym,spaceworld)
4218    end if
4219 
4220    call timab(129,2,tsec)
4221  end if ! if kpt parallel
4222 !=== End communications ==================
4223 
4224  call timab(127,1,tsec)
4225 
4226 !Normalize the fixed part of fermie1
4227  invfe1norm = zero ; if (abs(fe1norm) > tol10) invfe1norm=one/fe1norm
4228  fe1fixed=fe1fixed*invfe1norm
4229 
4230 
4231  if(nspden==4) then
4232 ! FR SPr symrhg will manage correctly this rearrangement
4233    rhorfermi(:,2)=rhorfermi(:,2)+(rhorfermi(:,1)+rhorfermi(:,4))    !(n+mx)
4234    rhorfermi(:,3)=rhorfermi(:,3)+(rhorfermi(:,1)+rhorfermi(:,4))    !(n+my)
4235    call timab(17,2,tsec)
4236  end if
4237 
4238 !Symmetrize the density
4239  call status(0,dtfil%filstat,iexit,level,'call symrhg   ')
4240 !In order to have the symrhg working in parallel on FFT coefficients, the size
4241 !of irzzon1 and phnons1 should be set to nfftot. Therefore, nsym\=1 does not work.
4242 !We also have the spin-up density, symmetrized, in rhorfermi(:,2).
4243  ABI_ALLOCATE(rhogfermi,(2,dtset%nfft))
4244  if (psps%usepaw==0) then
4245    call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,nspden,&
4246 &   nsppol,nsym1,dtset%paral_kgb,phnons1,rhogfermi,rhorfermi,rprimd,symaf1,symrl1)
4247  else
4248    call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,nspden,&
4249 &   nsppol,nsym1,dtset%paral_kgb,phnons1,rhogfermi,rhowfr   ,rprimd,symaf1,symrl1)
4250  end if
4251 
4252 !PAW: Build new rhoij quantities then symetrize them
4253 !Compute and add the compensation density to rhowfr to get the total density
4254  if (psps%usepaw == 1) then
4255    if (size(nhatfermi)>0) then
4256      call pawmkrho(1,arg,cplex,gprimd,0,indsy1,0,mpi_enreg,&
4257 &     my_natom,natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,pawfgr,&
4258 &     pawfgrtab,-10001,pawrhoijfermi,pawrhoijfermi_unsym,pawtab,dtset%qptn,&
4259 &     rhogfermi,rhowfr,rhorfermi,rprimd,symaf1,symrc1,dtset%typat,ucvol,&
4260 &     dtset%usewvl,xred,pawang_sym=pawang1,pawnhat=nhatfermi)
4261    else
4262      call pawmkrho(1,arg,cplex,gprimd,0,indsy1,0,mpi_enreg,&
4263 &     my_natom,natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,pawfgr,&
4264 &     pawfgrtab,-10001,pawrhoijfermi,pawrhoijfermi_unsym,pawtab,dtset%qptn,&
4265 &     rhogfermi,rhowfr,rhorfermi,rprimd,symaf1,symrc1,dtset%typat,ucvol,&
4266 &     dtset%usewvl,xred,pawang_sym=pawang1)
4267    end if
4268    ABI_DEALLOCATE(rhowfr)
4269    call pawrhoij_free_unpacked(pawrhoijfermi_unsym)
4270    if (paral_atom) then
4271      call pawrhoij_free(pawrhoijfermi_unsym)
4272      ABI_DATATYPE_DEALLOCATE(pawrhoijfermi_unsym)
4273    end if
4274  end if
4275  ABI_DEALLOCATE(rhogfermi)
4276 
4277 !Normalize the Fermi level charge density (and associated PAW occupancies)
4278  rhorfermi(:,:)=invfe1norm*rhorfermi(:,:)
4279  if (psps%usepaw==1) then
4280    if (size(nhatfermi)>0) nhatfermi(:,:)=invfe1norm*nhatfermi(:,:)
4281    do iatom=1,my_natom
4282      do ispden=1,nspden
4283        do ilmn=1,pawrhoijfermi(iatom)%nrhoijsel
4284          pawrhoijfermi(iatom)%rhoijp(ilmn,ispden)=&
4285 &         pawrhoijfermi(iatom)%rhoijp(ilmn,ispden)*invfe1norm
4286        end do
4287      end do
4288    end do
4289  end if
4290 
4291  call timab(127,2,tsec)
4292  call timab(121,2,tsec)
4293 
4294  DBG_EXIT('COLL')
4295 
4296 end subroutine dfpt_rhofermi

ABINIT/dfpt_scfcv [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_scfcv

FUNCTION

 Conducts set of passes or overall iterations of preconditioned
 conjugate gradient algorithm to converge wavefunctions to
 optimum and optionally to compute mixed derivatives of energy.

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  cg(2,mpw*nspinor*mband*mkmem*nsppol)=pw coefficients of GS wavefunctions at k.
  cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
              projected with non-local projectors: cprj=<p_i|Cnk>
  cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q
              projected with non-local projectors: cprjq=<p_i|Cnk+q>
  cpus= cpu time limit in seconds
  doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
  docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eew=2nd derivative of Ewald energy (hartree)
  efrhar=Contribution from frozen-wavefunction, hartree energy,
           to the second-derivative of total energy.
  efrkin=Contribution from frozen-wavefunction, kinetic energy,
           to the second-derivative of total energy.
  efrloc=Contribution from frozen-wavefunction, local potential,
           to the second-derivative of total energy.
  efrnl=Contribution from frozen-wavefunction, non-local potential,
           to the second-derivative of total energy.
  efrx1=Contribution from frozen-wavefunction, xc core correction(1),
           to the second-derivative of total energy.
  efrx2=Contribution from frozen-wavefunction, xc core correction(2),
           to the second-derivative of total energy.
  eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
  eii=2nd derivative of pseudopotential core energy (hartree)
  evdw=DFT-D semi-empirical part of 2nd-order total energy
  fermie=fermi energy (Hartree)
  hdr <type(hdr_type)>=the header of wf, den and pot files
  idir=direction of the current perturbation
  indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  ipert=type of the perturbation
  irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for RF symmetries
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates at k
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
  kxc(nfftf,nkxc)=exchange and correlation kernel (see rhotoxc.f)
  mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90)
  mkmem =number of k points treated by this node (GS data)
  mkqmem =number of k+q points which can fit in memory (GS data); 0 if use disk
  mk1mem =number of k points which can fit in memory (RF data); 0 if use disk
  mpert=maximum number of ipert
  mpw=maximum dimensioned size of npw for wfs at k.
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  nattyp(ntypat)= # atoms of each type.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point, for each polarization
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
  ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90)
  nkpt=number of k points in the full BZ
  nkpt_rbz=number of k points in the reduced BZ for this perturbation
  nkxc=second dimension of the kxc array.
  mpi_enreg=informations about MPI parallelization
  my_natom=number of atoms treated by current processor
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nspinor=number of spinorial components of the wavefunctions
  nsym1=number of symmetry elements in space group consistent with perturbation
  n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, nfftf
  occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
   at each k+q point of the reduced Brillouin zone.
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
   at each k point of the reduced Brillouin zone.
  paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS
  paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawang1 <type(pawang_type)>=pawang datastr. containing only symmetries preserving the perturbation
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  pertcase=fuill index of the perturbation
  phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic transl. phases, for RF symmetries
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information for the "fine" grid
  prtbbb=if 1, band-by-band decomposition (also dim of d2bbb)
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  qphon(3)=reduced coordinates for the phonon wavelength
  rhog(2,nfftf)=array for Fourier transform of GS electron density
  rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  symrl1(3,3,nsym1)=symmetry operations in real space in terms
   of primitive translations
  usecprj= 1 if cprj, cprjq arrays are stored in memory
  useylmgr = 1 if ylmgr  array is allocated
  useylmgr1= 1 if ylmgr1 array is allocated
  ddk<wfk_t>=ddk file
  vpsp1(cplex*nfftf)=first-order derivative of the ionic potential
  vtrial(nfftf,nspden)=GS potential (Hartree).
  vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree)
  wtk_rbz(nkpt_rbz)=weight for each k point in the reduced Brillouin zone
  xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm*useylmgr)= gradients of real spherical harmonics at k
  ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm*useylmgr1)= gradients of real spherical harmonics at k+q

OUTPUT

  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
  cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
    wavefunctions at k,q. They are orthogonalized to the active.
  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
       second order derivatives
  d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs
  d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs
  d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs
  eberry=energy associated with Berry phase
  edocc=correction to 2nd-order total energy coming from changes of occupation
  eeig0=0th-order eigenenergies part of 2nd-order total energy
  ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
    for strain perturbation only (zero otherwise, and not used)
  ehart1=1st-order Hartree part of 2nd-order total energy
  eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues (hartree)
  ek0=0th-order kinetic energy part of 2nd-order total energy.
  ek1=1st-order kinetic energy part of 2nd-order total energy.
  eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
  elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
  enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
  enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
  eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
        PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
  epaw1=1st-order PAW on-site part of 2nd-order total energy.
  etotal=total energy (sum of 7 contributions) (hartree)
  exc1=1st-order exchange-correlation part of 2nd-order total energy.
  gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
  gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}|\Psi^{(1)}>
      The wavefunction is orthogonal to the active space (for metals). It is not
      coherent with cg1.
  resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points
   of the reduced Brillouin zone, and spins
  residm=maximum value from resid array (except for nbdbuf highest bands)
  conv_retcode=return code, 0 if convergence was achieved.

SIDE EFFECTS

  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=updated wavefunctions (ortho. to occ. states);
  initialized= if 0 the initialization of the RF run is not yet finished
  mpi_enreg=informations about MPI parallelization
  rhog1(2,nfftf)=array for Fourier transform of RF electron density
  rhor1(cplex*nfftf,nspden)=array for RF electron density in electrons/bohr**3.
  === if psps%usepaw==1
    pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data

PARENTS

      dfpt_looppert

CHILDREN

      ab7_mixing_deallocate,ab7_mixing_new,ab7_mixing_use_disk_cache,appdig
      calcdensph,destroy_efield,dfpt_etot,dfpt_newvtr,dfpt_nselt,dfpt_nstdy
      dfpt_nstpaw,dfpt_rhofermi,dfpt_rhotov,dfpt_vtorho,dfptff_bec,dfptff_die
      dfptff_ebp,dfptff_edie,dfptff_initberry,fftdatar_write_from_hdr,fourdp
      getcut,hdr_update,metric,newfermie1,paw_an_free,paw_an_init
      paw_an_nullify,paw_an_reset_flags,paw_ij_free,paw_ij_init
      paw_ij_nullify,paw_ij_reset_flags,pawcprj_alloc,pawcprj_free
      pawcprj_getdim,pawdenpot,pawdij,pawdijfr,pawmknhat,pawnhatfr
      pawrhoij_alloc,pawrhoij_free,qmatrix,rf2_getidirs,scprqt,status,symdij
      timab,wfk_close,wrtout,xmpi_barrier,xmpi_isum,xmpi_sum,xmpi_wait

SOURCE

 271 subroutine dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,&
 272 &  dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset,&
 273 &  d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
 274 &  ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,&
 275 &  enl0,enl1,eovl1,epaw1,etotal,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,indkpt1,&
 276 &  indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
 277 &  kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,&
 278 &  mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,nattyp,nband_rbz,ncpgr,&
 279 &  nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,&
 280 &  nsym1,n3xccc,occkq,occ_rbz,&
 281 &  paw_an,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawrhoij1,pawtab,&
 282 &  pertcase,phnons1,ph1d,ph1df,&
 283 &  prtbbb,psps,qphon,resid,residm,rhog,rhog1,&
 284 &  rhor,rhor1,rprimd,symaf1,symrc1,symrl1,&
 285 &  usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,&
 286 &  wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,conv_retcode,&
 287 &  kramers_deg,&
 288 &  cg_mq,cg1_mq,cg1_active_mq,docckde_mq,eigen_mq,eigen1_mq,gh0c1_set_mq,gh1c_set_mq,&
 289 &  kg1_mq,npwar1_mq,occk_mq,resid_mq,residm_mq,rhog1_pq,rhog1_mq,rhor1_pq,rhor1_mq)
 290 
 291 
 292 !This section has been created automatically by the script Abilint (TD).
 293 !Do not modify the following lines by hand.
 294 #undef ABI_FUNC
 295 #define ABI_FUNC 'dfpt_scfcv'
 296 !End of the abilint section
 297 
 298  implicit none
 299 
 300 !Arguments ------------------------------------
 301  type(dataset_type),intent(in) :: dtset
 302  type(pseudopotential_type),intent(in) :: psps
 303  integer,intent(in) :: cplex,dim_eig2rf,idir,ipert,mgfftf,mk1mem,mkmem,mkqmem
 304  integer,intent(in) :: mpert,mpw,mpw1,my_natom,n3xccc,ncpgr,nfftf
 305  integer,intent(in) :: mpw1_mq !-q duplicate
 306  integer,intent(in) :: nkpt,nkpt_rbz,nkxc,nspden
 307  integer,intent(in) :: nsym1,pertcase,prtbbb,usecprj,useylmgr,useylmgr1
 308  logical,intent(in) :: kramers_deg
 309  integer,intent(inout) :: initialized
 310 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise
 311  integer,intent(in) :: atindx(dtset%natom)
 312  integer,intent(out) :: blkflg(3,mpert,3,mpert)
 313  integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom)
 314  integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 315  integer,intent(in) :: istwfk_rbz(nkpt_rbz)
 316  integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem),nattyp(psps%ntypat)
 317  integer,intent(in) :: nband_rbz(nkpt_rbz*dtset%nsppol)
 318  integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz)
 319  integer,optional,intent(in) :: npwar1_mq(nkpt_rbz)     !-q duplicate
 320  integer,optional,intent(in) :: kg1_mq(3,mpw1_mq*mk1mem)!
 321  integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1)
 322  integer,intent(out) :: conv_retcode
 323  real(dp),intent(in) :: cpus,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,eii
 324  real(dp),intent(out) :: eberry,edocc,eeig0,ehart01,ehart1,ek0,ek1,eloc0,elpsp1,enl0
 325  real(dp),intent(out) :: enl1,eovl1,epaw1,etotal,evdw,exc1,residm
 326  real(dp),optional,intent(out) :: residm_mq       !-q duplicate
 327  real(dp),intent(inout) :: fermie
 328  real(dp),intent(in) :: qphon(3)
 329 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise
 330  integer,intent(in) :: ngfftf(18)
 331  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*dtset%nsppol)
 332  real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol)
 333  real(dp),intent(out) :: cg1_active(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf)
 334  real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf)
 335  real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf)
 336  real(dp),intent(in)  :: cgq(2,mpw1*dtset%nspinor*dtset%mband*mkqmem*dtset%nsppol)
 337  real(dp),optional,intent(inout) :: cg1_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol)                  !start -q duplicates
 338  real(dp),optional,intent(out)   :: cg1_active_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf)!
 339  real(dp),optional,intent(out)   :: gh1c_set_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf)  !
 340  real(dp),optional,intent(out)   :: gh0c1_set_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*dim_eig2rf) !
 341  real(dp),optional,intent(in)    :: cg_mq(2,mpw1_mq*dtset%nspinor*dtset%mband*mkqmem*dtset%nsppol)                   !
 342  real(dp),optional,intent(in)    :: eigen_mq(dtset%mband*nkpt_rbz*dtset%nsppol)                                      !
 343  real(dp),optional,intent(in)    :: docckde_mq(dtset%mband*nkpt_rbz*dtset%nsppol)                                    !
 344  real(dp),optional,intent(out)   :: eigen1_mq(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)                       !
 345  real(dp),optional,intent(in)    :: occk_mq(dtset%mband*nkpt_rbz*dtset%nsppol)                                       !
 346  real(dp),optional,intent(out)   :: resid_mq(dtset%mband*nkpt_rbz*nspden)                                            !end
 347  real(dp),intent(out) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)
 348  real(dp),intent(out) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert)
 349  real(dp),intent(out) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw)
 350  real(dp),intent(in) :: dielt(3,3)
 351  real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*dtset%nsppol)
 352  real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*dtset%nsppol)
 353  real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*dtset%nsppol)
 354  real(dp),intent(out) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)
 355  real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*dtset%nsppol)
 356  real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),kxc(nfftf,nkxc)
 357  real(dp),intent(in) :: nhat(nfftf,dtset%nspden)
 358  real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*dtset%nsppol)
 359  real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol)
 360  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom),ph1df(2,3*(2*mgfftf+1)*dtset%natom)
 361  real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 362  real(dp),intent(out) :: resid(dtset%mband*nkpt_rbz*nspden)
 363  real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),rprimd(3,3)
 364  real(dp),intent(inout) :: rhog1(2,nfftf),rhor1(cplex*nfftf,nspden),xred(3,dtset%natom)
 365  real(dp),optional,intent(inout) :: rhog1_pq(2,nfftf),rhor1_pq(cplex*nfftf,nspden)                                 !+q/-q duplicates
 366  real(dp),optional,intent(inout) :: rhog1_mq(2,nfftf),rhor1_mq(cplex*nfftf,nspden)                                 !
 367  real(dp),target,intent(in) :: vtrial(nfftf,nspden)
 368  real(dp),intent(in) :: vpsp1(cplex*nfftf),vxc(nfftf,nspden)
 369  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xccc3d1(cplex*n3xccc)
 370  real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 371  real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
 372  real(dp),intent(in) :: ylmgr(mpw*mkmem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
 373  real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-dtset%natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
 374  real(dp),intent(in) :: zeff(3,3,dtset%natom)
 375  type(pawcprj_type),intent(in) :: cprj(dtset%natom,dtset%nspinor*dtset%mband*mkmem*dtset%nsppol*usecprj)
 376  type(pawcprj_type),intent(in) :: cprjq(dtset%natom,dtset%nspinor*dtset%mband*mkqmem*dtset%nsppol*usecprj)
 377  type(datafiles_type),intent(in) :: dtfil
 378  type(hdr_type),intent(inout) :: hdr
 379  type(pawang_type),intent(in) :: pawang,pawang1
 380  type(pawfgr_type),intent(in) :: pawfgr
 381  type(paw_an_type),intent(in) :: paw_an(my_natom*psps%usepaw)
 382  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
 383  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
 384  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 385  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw)
 386  type(pawrhoij_type),intent(inout) :: pawrhoij1(my_natom*psps%usepaw)
 387  type(pawtab_type), intent(in) :: pawtab(psps%ntypat*psps%usepaw)
 388  type(MPI_type),intent(inout) :: mpi_enreg
 389  type(wfk_t),intent(inout) :: ddk_f(4)
 390 
 391 !Local variables-------------------------------
 392 !scalars
 393  integer,parameter :: level=12,response=1
 394  integer :: afford,bantot_rbz,choice,cplex_rhoij,dbl_nnsclo
 395  integer :: has_dijfr,has_diju,iatom,ider,idir_dum,idir_paw1,ierr,iexit,errid,denpot
 396  integer :: iprcel,iscf10_mod,iscf_mod,ispden,ispmix
 397  integer :: istep,itypat,izero,lmn2_size,me,mgfftdiel,mvdum
 398  integer :: nfftdiel,nfftmix,nfftotf,nhat1grdim,npawmix,npwdiel,nspden_rhoij,nstep,nzlmopt
 399  integer :: optene,optfr,option,optres,prtfor,quit,quit_sum,qzero
 400  integer :: my_quit,quitsum_request,timelimit_exit,varid,ncerr,ncid
 401  integer ABI_ASYNC :: quitsum_async
 402  integer :: rdwrpaw,spaceComm,sz1,sz2,usexcnhat,Z_kappa
 403  integer :: dbl_nnsclo_mq,ifft !-q duplicate for dbl_nnsclo
 404 !integer :: pqmq ! pqmq = indicator for potential mixing
 405  logical :: need_fermie1,paral_atom,use_nhat_gga
 406  real(dp) :: wtime_step,now,prev
 407  real(dp) :: born,born_bar,boxcut,deltae,diffor,diel_q,dum,ecut,ecutf,elast
 408  real(dp) :: epawdc1_dum,evar,fe1fixed,fermie1,gsqcut,qphon_norm,maxfor,renorm,res2,res3,residm2
 409  real(dp) :: ucvol,vxcavg,elmag1
 410  real(dp) :: res2_mq,fe1fixed_mq,elast_mq
 411  real(dp) :: eberry_mq,edocc_mq,eeig0_mq,ehart01_mq,ehart1_mq,ek0_mq,ek1_mq,eloc0_mq,elpsp1_mq,enl0_mq
 412  real(dp) :: enl1_mq,eovl1_mq,epaw1_mq,exc1_mq,fermie1_mq,deltae_mq,elmag1_mq
 413  character(len=500) :: msg
 414  character(len=fnlen) :: fi1o
 415 !character(len=fnlen) :: fi1o_vtk
 416  integer  :: prtopt
 417  type(ab7_mixing_object) :: mix
 418  type(efield_type) :: dtefield
 419 !arrays
 420  integer :: ngfftmix(18)
 421  integer,allocatable :: dimcprj(:),pwindall(:,:,:)
 422  integer,pointer :: my_atmtab(:)
 423  real(dp) :: dielar(7)
 424  real(dp) :: favg(3),gmet(3,3),gprimd(3,3),q_cart(3),qphon2(3),qred2cart(3,3)
 425  real(dp) :: rmet(3,3),tollist(12),tsec(2)
 426  real(dp) :: zeff_red(3),zeff_bar(3,3)
 427  real(dp) :: intgden(dtset%nspden,dtset%natom),dentot(dtset%nspden)
 428 !real(dp) :: zdmc_red(3),zdmc_bar(3,3),mean_rhor1(1) !dynamic magnetic charges and mean density
 429  real(dp),allocatable :: dielinv(:,:,:,:,:)
 430  real(dp),allocatable :: fcart(:,:),nhat1(:,:),nhat1gr(:,:,:),nhatfermi(:,:),nvresid1(:,:),nvresid2(:,:)
 431  real(dp),allocatable :: qmat(:,:,:,:,:,:),resid2(:),rhog2(:,:),rhor2(:,:),rhorfermi(:,:)
 432  real(dp),allocatable :: susmat(:,:,:,:,:),vhartr1(:),vxc1(:,:)
 433  real(dp),allocatable :: vhartr1_tmp(:,:)
 434  real(dp),allocatable,target :: vtrial1(:,:),vtrial2(:,:)
 435  real(dp),allocatable :: vtrial1_pq(:,:),vtrial1_mq(:,:),rhorfermi_mq(:,:)
 436  real(dp),allocatable :: nvresid1_mq(:,:),vxc1_mq(:,:),vhartr1_mq(:)
 437  real(dp),pointer :: vtrial1_tmp(:,:)
 438  type(pawcprj_type),allocatable :: cprj1(:,:)
 439  type(paw_an_type),allocatable :: paw_an1(:)
 440  type(paw_ij_type),allocatable :: paw_ij1(:)
 441  type(pawrhoij_type),allocatable :: pawrhoijfermi(:)
 442 
 443 ! *********************************************************************
 444 
 445  DBG_ENTER("COLL")
 446 
 447  call timab(120,1,tsec)
 448  call timab(154,1,tsec)
 449 
 450  call status(0,dtfil%filstat,iexit,level,'init          ')
 451 
 452  ! intel 18 really needs this to be initialized
 453  maxfor = zero
 454 
 455  ! enable time limit handler if not done in callers.
 456  if (enable_timelimit_in(ABI_FUNC) == ABI_FUNC) then
 457    write(std_out,*)"Enabling timelimit check in function: ",trim(ABI_FUNC)," with timelimit: ",trim(sec2str(get_timelimit()))
 458  end if
 459 
 460 !Parallelism data
 461  spaceComm=mpi_enreg%comm_cell
 462  me=mpi_enreg%me_kpt
 463  paral_atom=(my_natom/=dtset%natom)
 464  my_atmtab=>mpi_enreg%my_atmtab
 465 
 466  _IBM6("XLF in dfpt_scfcv")
 467 
 468 !Save some variables from dataset definition
 469  ecut=dtset%ecut
 470  ecutf=ecut;if (psps%usepaw==1) ecutf=dtset%pawecutdg
 471  iprcel=dtset%iprcel
 472  tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr
 473  tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe
 474  tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff
 475  nfftotf=product(ngfftf(1:3))
 476  nstep=dtset%nstep
 477  iscf_mod=dtset%iscf
 478  iscf10_mod=mod(iscf_mod,10)
 479 
 480  qzero=0; if(qphon(1)**2+qphon(2)**2+qphon(3)**2 < tol14) qzero=1
 481 
 482  need_fermie1=((qzero==1.and.dtset%frzfermi==0.and.nstep>0).and.&
 483 & (dtset%occopt>=3.and.dtset%occopt<=8).and. &
 484 & (ipert<=dtset%natom.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or.ipert==dtset%natom+5))
 485 
 486 !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f
 487  if (ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11) iscf_mod=-3
 488 
 489 !Compute different geometric tensor, as well as ucvol, from rprimd
 490  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 491 
 492 !Compute large sphere cut-off gsqcut
 493  qphon2(:)=zero;if (psps%usepaw==1) qphon2(:)=qphon(:)
 494  call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,qphon2,ngfftf)
 495 
 496 !Some variables need to be initialized/nullify at start
 497  quit=0 ; dbl_nnsclo=0 ; elast=zero; conv_retcode = -1
 498  optres=merge(0,1,abs(iscf_mod)<10)
 499  usexcnhat=0
 500 !This might be taken away later
 501  edocc=zero ; eeig0=zero ; ehart01=zero ; ehart1=zero ; ek0=zero ; ek1=zero
 502  eloc0=zero ; elpsp1=zero ; enl0=zero ; enl1=zero ; eovl1=zero; exc1=zero
 503  deltae=zero ; fermie1=zero ; epaw1=zero ; eberry=zero ; elmag1=zero
 504  elast_mq=zero
 505  dbl_nnsclo_mq=0
 506 !This might be taken away later
 507  edocc_mq=zero ; eeig0_mq=zero ; ehart01_mq=zero ; ehart1_mq=zero ; ek0_mq=zero ; ek1_mq=zero
 508  eloc0_mq=zero ; elpsp1_mq=zero ; enl0_mq=zero ; enl1_mq=zero ; eovl1_mq=zero; exc1_mq=zero
 509  deltae_mq=zero ; fermie1_mq=zero ; epaw1_mq=zero ; eberry_mq=zero ; elmag1_mq=zero
 510  res2_mq=zero
 511 
 512 !Examine tolerance criteria, and eventually  print a line to the output
 513 !file (with choice=1, the only non-dummy arguments of scprqt are
 514 !nstep, tollist and iscf - still, diffor,res2,prtfor,fcart are here initialized to 0)
 515  choice=1 ; prtfor=0 ; diffor=zero ; res2=zero
 516  ABI_ALLOCATE(fcart,(3,dtset%natom))
 517 
 518  call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
 519 & etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),&
 520 & 1,iscf_mod,istep,kpt_rbz,maxfor,&
 521 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
 522 & nstep,occ_rbz,0,prtfor,0,&
 523 & quit,res2,resid,residm,response,&
 524 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode)
 525 
 526 !Allocations/initializations for PAW only
 527  if(psps%usepaw==1) then
 528    usexcnhat=maxval(pawtab(:)%usexcnhat)
 529    use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0)
 530 !  1st-order compensation density
 531    ABI_ALLOCATE(nhat1,(cplex*nfftf,dtset%nspden))
 532    nhat1=zero
 533 !  Projections of 1-st order WF on nl projectors
 534    ABI_DATATYPE_ALLOCATE(cprj1,(dtset%natom,dtset%nspinor*dtset%mband*mk1mem*dtset%nsppol*usecprj))
 535    if (usecprj==1.and.mk1mem/=0) then
 536      !cprj ordered by atom-type
 537      ABI_ALLOCATE(dimcprj,(dtset%natom))
 538      call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
 539      call pawcprj_alloc(cprj1,0,dimcprj)
 540      ABI_DEALLOCATE(dimcprj)
 541    end if
 542 !  1st-order arrays/variables related to the PAW spheres
 543    ABI_DATATYPE_ALLOCATE(paw_an1,(my_natom))
 544    ABI_DATATYPE_ALLOCATE(paw_ij1,(my_natom))
 545    call paw_an_nullify(paw_an1)
 546    call paw_ij_nullify(paw_ij1)
 547 
 548    has_dijfr=0;if (ipert/=dtset%natom+1.and.ipert/=dtset%natom+10) has_dijfr=1
 549    has_diju=0;if (dtset%usepawu==1) has_diju=1
 550    call paw_an_init(paw_an1,dtset%natom,dtset%ntypat,0,0,dtset%nspden,&
 551 &   cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,&
 552 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 553    call paw_ij_init(paw_ij1,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,&
 554 &   dtset%ntypat,dtset%typat,pawtab,&
 555 &   has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,has_dijU=has_diju,&
 556 &   mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom)
 557  else
 558    ABI_ALLOCATE(nhat1,(0,0))
 559    ABI_DATATYPE_ALLOCATE(cprj1,(0,0))
 560    ABI_DATATYPE_ALLOCATE(paw_an1,(0))
 561    ABI_DATATYPE_ALLOCATE(paw_ij1,(0))
 562  end if ! PAW
 563 
 564 !Various allocations (potentials)
 565  ABI_ALLOCATE(vhartr1,(cplex*nfftf))
 566  ABI_ALLOCATE(vtrial1,(cplex*nfftf,nspden))
 567  if(.not.kramers_deg) then
 568    ABI_ALLOCATE(vhartr1_mq,(cplex*nfftf))
 569    ABI_ALLOCATE(vtrial1_pq,(cplex*nfftf,nspden))
 570    ABI_ALLOCATE(vtrial1_mq,(cplex*nfftf,nspden))
 571  end if
 572 ! TODO: for non collinear case this should always be nspden, in NCPP case as well!!!
 573  ABI_ALLOCATE(vxc1,(cplex*nfftf,nspden*(1-usexcnhat))) ! Not always needed
 574  vtrial1_tmp => vtrial1   ! this is to avoid errors when vtrial1_tmp is unused
 575 
 576  if (.not.kramers_deg) then
 577    ABI_ALLOCATE(vxc1_mq,(cplex*nfftf,nspden*(1-usexcnhat)))
 578  end if
 579 
 580 !Several parameters and arrays for the SCF mixing:
 581 !These arrays are needed only in the self-consistent case
 582  if (iscf_mod>0.or.iscf_mod==-3) then
 583    ABI_ALLOCATE(nvresid1,(cplex*nfftf,dtset%nspden))
 584    if (nstep==0) nvresid1=zero
 585    if ((dtset%getddb .ne. 0 .or. dtset%irdddb .ne.0) .and. qzero .ne. 1) then
 586      ABI_ALLOCATE(nvresid2,(cplex*nfftf,dtset%nspden))
 587      if (nstep==0) nvresid2=zero
 588    end if
 589    if (.not.kramers_deg) then
 590      ABI_ALLOCATE(nvresid1_mq,(cplex*nfftf,dtset%nspden))
 591      if (nstep==0) nvresid1_mq=zero
 592    end if
 593  else
 594    ABI_ALLOCATE(nvresid1,(0,0))
 595    if(.not.kramers_deg) then
 596      ABI_ALLOCATE(nvresid1_mq,(0,0))
 597    end if
 598  end if
 599  if(nstep>0 .and. iscf_mod>0) then
 600    dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
 601    dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
 602    dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
 603    dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag
 604 !  Additional allocation for mixing within PAW
 605    npawmix=0
 606    if(psps%usepaw==1) then
 607      do iatom=1,my_natom
 608        itypat=pawrhoij1(iatom)%itypat
 609        lmn2_size=pawtab(itypat)%lmn2_size
 610        pawrhoij1(iatom)%use_rhoijres=1
 611        sz1=pawrhoij1(iatom)%cplex*lmn2_size;sz2=pawrhoij1(iatom)%nspden
 612        ABI_ALLOCATE(pawrhoij1(iatom)%rhoijres,(sz1,sz2))
 613        do ispden=1,pawrhoij1(iatom)%nspden
 614          pawrhoij1(iatom)%rhoijres(:,ispden)=zero
 615        end do
 616        ABI_ALLOCATE(pawrhoij1(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz))
 617        pawrhoij1(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz
 618        pawrhoij1(iatom)%kpawmix=pawtab(itypat)%kmix
 619        npawmix=npawmix+pawrhoij1(iatom)%nspden*pawtab(itypat)%lmnmix_sz*pawrhoij1(iatom)%cplex
 620      end do
 621    end if
 622    denpot = AB7_MIXING_POTENTIAL
 623    if (dtset%iscf > 10) denpot = AB7_MIXING_DENSITY
 624    if (psps%usepaw==1.and.dtset%pawmixdg==0) then
 625      ispmix=AB7_MIXING_FOURRIER_SPACE;nfftmix=dtset%nfft;ngfftmix(:)=dtset%ngfft(:)
 626    else
 627      ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:)
 628    end if
 629    if (iscf10_mod == 5 .or. iscf10_mod == 6) then
 630      call ab7_mixing_new(mix, iscf10_mod, denpot, cplex, &
 631 &     nfftf, dtset%nspden, npawmix, errid, msg, dtset%npulayit)
 632    else
 633      call ab7_mixing_new(mix, iscf10_mod, denpot, max(cplex, ispmix), &
 634 &     nfftmix, dtset%nspden, npawmix, errid, msg, dtset%npulayit)
 635    end if
 636    if (errid /= AB7_NO_ERROR) then
 637      MSG_ERROR(msg)
 638    end if
 639    if (dtset%mffmem == 0) then
 640      call ab7_mixing_use_disk_cache(mix, dtfil%fnametmp_fft)
 641    end if
 642  end if ! iscf, nstep
 643 
 644 !Here, allocate arrays for computation of susceptibility and dielectric matrix or for TDDFT
 645  if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then
 646 !  Here, for TDDFT, artificially set iprcel . Also set a variable to reduce the memory needs.
 647    afford=1
 648    if(iscf_mod==-1) then
 649      iprcel=21
 650      afford=0
 651    end if
 652    npwdiel=1
 653    mgfftdiel=1
 654    nfftdiel=1
 655 !  Now, performs allocation
 656 !  CAUTION : the dimensions are still those of GS, except for phnonsdiel
 657    ABI_ALLOCATE(dielinv,(2,npwdiel*afford,nspden,npwdiel,nspden))
 658    ABI_ALLOCATE(susmat,(2,npwdiel*afford,nspden,npwdiel,nspden))
 659  end if
 660 
 661 !Initialize Berry-phase related stuffs
 662  if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
 663 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
 664    ABI_ALLOCATE(pwindall,(max(mpw,mpw1)*mkmem,8,3))
 665    call dfptff_initberry(dtefield,dtset,gmet,kg,kg1,dtset%mband,mkmem,mpi_enreg,&
 666 &   mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,occ_rbz,pwindall,rprimd)
 667 !  calculate inverse of the overlap matrix
 668    ABI_ALLOCATE(qmat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3))
 669    call qmatrix(cg,dtefield,qmat,mpw,mpw1,mkmem,dtset%mband,npwarr,nkpt,dtset%nspinor,dtset%nsppol,pwindall)
 670  else
 671    ABI_ALLOCATE(pwindall,(0,0,0))
 672    ABI_ALLOCATE(qmat,(0,0,0,0,0,0))
 673  end if
 674 
 675  call timab(154,2,tsec)
 676 
 677 !######################################################################
 678 !PERFORM ELECTRONIC ITERATIONS
 679 !######################################################################
 680 
 681 !Offer option of computing 2nd-order total energy with existing
 682 !wavefunctions when nstep<=0, else do nstep iterations
 683 !Note that for non-self-consistent calculations, this loop will be exited
 684 !after the first call to dfpt_vtorho
 685 
 686 !Pass through the first routines even when nstep==0
 687 !write(std_out,*) 'dfpt_scfcv, nstep=', max(1,nstep)
 688 
 689  quitsum_request = xmpi_request_null; timelimit_exit = 0
 690 
 691  do istep=1,max(1,nstep)
 692 
 693    ! Handle time limit condition.
 694    if (istep == 1) prev = abi_wtime()
 695    if (istep  > 1) then
 696      now = abi_wtime()
 697      wtime_step = now - prev
 698      prev = now
 699      call wrtout(std_out,sjoin("dfpt_scfcv: previous iteration took ",sec2str(wtime_step)))
 700 
 701      if (have_timelimit_in(ABI_FUNC)) then
 702        if (istep > 2) then
 703          call xmpi_wait(quitsum_request,ierr)
 704          if (quitsum_async > 0) then
 705            write(msg,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),". Will exit istep loop in dfpt_scfcv."
 706            MSG_COMMENT(msg)
 707            call wrtout(ab_out, msg, "COLL")
 708            timelimit_exit = 1
 709            exit
 710          end if
 711        end if
 712 
 713        my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1
 714        call xmpi_isum(my_quit,quitsum_async,spacecomm,quitsum_request,ierr)
 715      end if
 716    end if
 717 
 718 !  ######################################################################
 719 !  The following steps are done once
 720 !  ----------------------------------------------------------------------
 721    if (istep==1)then
 722 
 723 !    PAW only: compute frozen part of 1st-order compensation density
 724 !    and frozen part of psp strengths Dij
 725 !    ----------------------------------------------------------------------
 726      if (psps%usepaw==1) then
 727        optfr=0
 728        idir_paw1 = idir
 729        if (ipert==dtset%natom+11) then
 730          call rf2_getidirs(idir,idir_dum,idir_paw1)
 731        end if
 732        call pawdijfr(cplex,gprimd,idir_paw1,ipert,my_natom,dtset%natom,nfftf,ngfftf,nspden,dtset%nsppol,&
 733 &       psps%ntypat,optfr,paw_ij1,pawang,pawfgrtab,pawrad,pawtab,qphon,&
 734 &       rprimd,ucvol,vpsp1,vtrial,vxc,xred,&
 735 &       mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 736 
 737        if ((iscf_mod>=0.or.usexcnhat==0).and.(dtset%pawstgylm/=0)) then
 738          ider=0;if ((ipert<=dtset%natom).and.(use_nhat_gga)) ider=1
 739          call pawnhatfr(ider,idir_paw1,ipert,my_natom,dtset%natom,nspden,psps%ntypat,&
 740 &         pawang,pawfgrtab,pawrhoij,pawtab,rprimd,&
 741 &         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 742        end if
 743      end if
 744 !    PAW only: we sometimes have to compute 1st-order compensation density
 745 !    and eventually add it to density from 1st-order WFs
 746 !    ----------------------------------------------------------------------
 747      nhat1grdim=0
 748      ABI_ALLOCATE(nhat1gr,(0,0,0))
 749      if (psps%usepaw==1.and.ipert/=dtset%natom+1.and.ipert/=dtset%natom+10) then
 750        call timab(564,1,tsec)
 751        nhat1grdim=0;if (dtset%xclevel==2) nhat1grdim=usexcnhat*dtset%pawnhatxc
 752        ider=2*nhat1grdim;izero=0
 753        if (nhat1grdim>0)   then
 754          ABI_DEALLOCATE(nhat1gr)
 755          ABI_ALLOCATE(nhat1gr,(cplex*nfftf,dtset%nspden,3*nhat1grdim))
 756        end if
 757        call pawmknhat(dum,cplex,ider,idir_paw1,ipert,izero,gprimd,my_natom,dtset%natom,&
 758 &       nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1,&
 759 &       pawrhoij1,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
 760 &       mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 761        if (dtfil%ireadwf/=0.and.dtset%get1den==0.and.dtset%ird1den==0.and.initialized==0) then
 762          rhor1(:,:)=rhor1(:,:)+nhat1(:,:)
 763          call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 764        end if
 765        call timab(564,2,tsec)
 766      end if
 767 !    Set initial guess for 1st-order potential
 768 !    ----------------------------------------------------------------------
 769      call status(istep,dtfil%filstat,iexit,level,'get vtrial1   ')
 770      option=1;optene=0;if (iscf_mod==-2) optene=1
 771      call dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,&
 772 &     dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,nhat1,nhat1gr,nhat1grdim,&
 773 &     nkxc,nspden,n3xccc,optene,option,dtset%paral_kgb,dtset%qptn,&
 774 &     rhog,rhog1,rhor,rhor1,rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,&
 775 &     nvresid1,res2,vtrial1,vxc,vxc1,xccc3d1,dtset%ixcrot)
 776 
 777      if(.not.kramers_deg) then
 778        vtrial1_pq=vtrial1 !save trial potential at +q
 779        !rhor1_mq=rhor1
 780        !rhog1_mq=rhog1
 781        !get initial guess for vtrial1 at -q
 782        do ifft=1,nfftf
 783          vtrial1_mq(2*ifft-1,1)=+vtrial1(2*ifft-1,1)
 784          vtrial1_mq(2*ifft-1,2)=+vtrial1(2*ifft-1,2)
 785          vtrial1_mq(2*ifft  ,1)=-vtrial1(2*ifft  ,1)
 786          vtrial1_mq(2*ifft  ,2)=-vtrial1(2*ifft  ,2)
 787          vtrial1_mq(2*ifft-1,3)= vtrial1(2*ifft  ,4) !Re[V^12]
 788          vtrial1_mq(2*ifft  ,3)= vtrial1(2*ifft-1,4) !Im[V^12],see definition of v(:,4) cplex=2 case
 789          vtrial1_mq(2*ifft  ,4)= vtrial1(2*ifft-1,3) !Re[V^21]=Re[V^12]
 790          vtrial1_mq(2*ifft-1,4)= vtrial1(2*ifft  ,3) !Re[V^21]=Re[V^12]
 791        end do
 792      end if
 793 
 794 !    For Q=0 and metallic occupation, initialize quantities needed to
 795 !    compute the first-order Fermi energy
 796 !    ----------------------------------------------------------------------
 797      if (need_fermie1) then
 798        ABI_ALLOCATE(rhorfermi,(cplex*nfftf,nspden))
 799        if(.not.kramers_deg) then
 800          ABI_ALLOCATE(rhorfermi_mq,(cplex*nfftf,nspden))
 801        end if
 802        if (psps%usepaw==1.and.usexcnhat==0) then
 803          ABI_ALLOCATE(nhatfermi,(cplex*nfftf,nspden))
 804        else
 805          ABI_ALLOCATE(nhatfermi,(0,0))
 806        end if
 807        ABI_DATATYPE_ALLOCATE(pawrhoijfermi,(my_natom*psps%usepaw))
 808        if (psps%usepaw==1) then
 809          cplex_rhoij=max(cplex,dtset%pawcpxocc)
 810          nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
 811          call pawrhoij_alloc(pawrhoijfermi,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
 812 &         dtset%nsppol,dtset%typat,pawtab=pawtab,mpi_atmtab=mpi_enreg%my_atmtab,&
 813 &         comm_atom=mpi_enreg%comm_atom)
 814        end if
 815 
 816        call dfpt_rhofermi(cg,cgq,cplex,cprj,cprjq,&
 817 &       doccde_rbz,docckqde,dtfil,dtset,eigenq,eigen0,eigen1,fe1fixed,gmet,gprimd,idir,&
 818 &       indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mkmem,mkqmem,mk1mem,mpi_enreg,&
 819 &       mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1,&
 820 &       nspden,dtset%nsppol,nsym1,occkq,occ_rbz,&
 821 &       paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,&
 822 &       phnons1,ph1d,dtset%prtvol,psps,rhorfermi,rmet,rprimd,symaf1,symrc1,symrl1,&
 823 &       ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1)
 824        if (.not.kramers_deg) then
 825          call dfpt_rhofermi(cg,cg_mq,cplex,cprj,cprjq,&
 826 &         doccde_rbz,docckde_mq,dtfil,dtset,eigen_mq,eigen0,eigen1_mq,fe1fixed_mq,gmet,gprimd,idir,&
 827 &         indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1_mq,kpt_rbz,dtset%mband,mkmem,mkqmem,mk1mem,mpi_enreg,&
 828 &         mpw,mpw1_mq,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1_mq,&
 829 &         nspden,dtset%nsppol,nsym1,occk_mq,occ_rbz,&
 830 &         paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,&
 831 &         phnons1,ph1d,dtset%prtvol,psps,rhorfermi_mq,rmet,rprimd,symaf1,symrc1,symrl1,&
 832 &         ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1)
 833        end if
 834 !        call calcdensph(gmet,mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,&
 835 !&        dtset%ntypat,ab_out,dtset%ratsph,rhorfermi,rprimd,dtset%typat,ucvol,xred,&
 836 !&        -1,cplex) !dbg
 837 
 838      end if
 839 
 840    end if ! End the condition of istep==1
 841 
 842 !  ######################################################################
 843 !  The following steps are done at every iteration
 844 !  ----------------------------------------------------------------------
 845 
 846    if (psps%usepaw==1)then
 847 !    Computation of "on-site" 2nd-order energy, first-order potentials, first-order densities
 848      nzlmopt=0;if (istep==2.and.dtset%pawnzlm>0) nzlmopt=-1
 849      if (istep>2) nzlmopt=dtset%pawnzlm
 850      call paw_an_reset_flags(paw_an1) ! Force the recomputation of on-site potentials
 851      call paw_ij_reset_flags(paw_ij1,self_consistent=.true.) ! Force the recomputation of Dij
 852      option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1
 853      call status(istep,dtfil%filstat,iexit,level,'call pawdenpot')
 854      call pawdenpot(dum,epaw1,epawdc1_dum,ipert,dtset%ixc,my_natom,dtset%natom,&
 855 &     dtset%nspden,psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an1,paw_an,paw_ij1,pawang,&
 856 &     dtset%pawprtvol,pawrad,pawrhoij1,dtset%pawspnorb,pawtab,dtset%pawxcdev,&
 857 &     dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, &
 858 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 859 
 860 !    First-order Dij computation
 861      call status(istep,dtfil%filstat,iexit,level,'call pawdij   ')
 862      call timab(561,1,tsec)
 863      if (has_dijfr>0) then
 864        !vpsp1 contribution to Dij already stored in frozen part of Dij
 865        ABI_ALLOCATE(vtrial1_tmp,(cplex*nfftf,nspden))
 866        vtrial1_tmp=vtrial1
 867        do ispden=1,min(dtset%nspden,2)
 868          vtrial1_tmp(:,ispden)=vtrial1_tmp(:,ispden)-vpsp1(:)
 869        end do
 870      else
 871        vtrial1_tmp => vtrial1
 872      end if
 873      call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,dtset%natom,&
 874 &     nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1,paw_ij1,pawang,&
 875 &     pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1,dtset%pawspnorb,pawtab,&
 876 &     dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%charge,vtrial1_tmp,vxc1,xred,&
 877 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 878      if (has_dijfr>0) then
 879        ABI_DEALLOCATE(vtrial1_tmp)
 880      end if
 881 
 882      call status(istep,dtfil%filstat,iexit,level,'call symdij   ')
 883      call symdij(gprimd,indsy1,ipert,my_natom,dtset%natom,nsym1,psps%ntypat,0,&
 884 &     paw_ij1,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, &
 885 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
 886 &     qphon=qphon)
 887      call timab(561,2,tsec)
 888    end if ! end usepaw section
 889 
 890 !  ######################################################################
 891 !  The following steps are done only when nstep>0
 892 !  ----------------------------------------------------------------------
 893    call status(istep,dtfil%filstat,iexit,level,'loop istep    ')
 894 
 895    if(iscf_mod>0.and.nstep>0)then
 896      write(msg, '(a,a,i4)' )ch10,' ITER STEP NUMBER  ',istep
 897      call wrtout(std_out,msg,'COLL')
 898    end if
 899 
 900 !  For Q=0 and metallic occupation, calculate the first-order Fermi energy
 901    if (need_fermie1) then
 902      call newfermie1(cplex,fermie1,fe1fixed,ipert,istep,dtset%ixc,my_natom,dtset%natom,&
 903 &     nfftf,nfftotf,nhatfermi,nspden,dtset%ntypat,dtset%occopt,paw_an,paw_an1,paw_ij1,pawang,&
 904 &     dtset%pawnzlm,pawrad,pawrhoij1,pawrhoijfermi,pawtab,dtset%pawxcdev,&
 905 &     dtset%prtvol,rhorfermi,ucvol,psps%usepaw,usexcnhat,vtrial1,vxc1,dtset%xclevel,&
 906 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 907      if (.not.kramers_deg) then
 908        !fermie1_mq is updated as well at "-q"
 909        call newfermie1(cplex,fermie1_mq,fe1fixed_mq,ipert,istep,dtset%ixc,my_natom,dtset%natom,&
 910 &       nfftf,nfftotf,nhatfermi,nspden,dtset%ntypat,dtset%occopt,paw_an,paw_an1,paw_ij1,pawang,&
 911 &       dtset%pawnzlm,pawrad,pawrhoij1,pawrhoijfermi,pawtab,dtset%pawxcdev,&
 912 &       dtset%prtvol,rhorfermi_mq,ucvol,psps%usepaw,usexcnhat,vtrial1_mq,vxc1_mq,dtset%xclevel,&
 913 &       mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 914      end if
 915    end if
 916 
 917 !  No need to continue and call dfpt_vtorho, when nstep==0
 918    if(nstep==0) exit
 919 
 920 !  #######################e1magh###############################################
 921 !  Compute the 1st-order density rho1 from the 1st-order trial potential
 922 !  ----------------------------------------------------------------------
 923 
 924    call dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,&
 925 &   dbl_nnsclo,dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,dtset%qptn,edocc,&
 926 &   eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,enl0,enl1,fermie1,gh0c1_set,gh1c_set,&
 927 &   gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,&
 928 &   mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,&
 929 &   nhat1,nkpt_rbz,npwarr,npwar1,res2,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid1,&
 930 &   occkq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,&
 931 &   pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid,residm,rhog1,&
 932 &   rhor1,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,usecprj,useylmgr1,ddk_f,&
 933 &   vtrial,vtrial1,wtk_rbz,xred,ylm,ylm1,ylmgr1)
 934    if (.not.kramers_deg) then
 935      rhor1_pq=rhor1 !at this stage rhor1_pq contains only one term of the 1st order density at +q
 936      rhog1_pq=rhog1 !same for rhog1_pq
 937      !get the second term related to 1st order wf at -q
 938      call dfpt_vtorho(cg,cg_mq,cg1_mq,cg1_active_mq,cplex,cprj,cprjq,cprj1,&
 939 &     dbl_nnsclo_mq,dim_eig2rf,doccde_rbz,docckde_mq,dtefield,dtfil,dtset,-dtset%qptn,edocc_mq,&
 940 &     eeig0_mq,eigen_mq,eigen0,eigen1_mq,ek0_mq,ek1_mq,eloc0_mq,enl0_mq,enl1_mq,fermie1_mq,gh0c1_set_mq,gh1c_set_mq,&
 941 &     gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1_mq,kpt_rbz,dtset%mband,&
 942 &     mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1_mq,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,&
 943 &     nhat1,nkpt_rbz,npwarr,npwar1_mq,res2_mq,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid1_mq,&
 944 &     occk_mq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,&
 945 &     pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid_mq,residm_mq,rhog1_mq,&
 946 &     rhor1_mq,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,usecprj,useylmgr1,ddk_f,&
 947 &     vtrial,vtrial1_mq,wtk_rbz,xred,ylm,ylm1,ylmgr1)
 948      !reconstruct the +q and -q densities, this might bug if fft parallelization is used, todo...
 949      do ifft=1,nfftf
 950        rhor1_pq(2*ifft-1,:) = half*(rhor1(2*ifft-1,:)+rhor1_mq(2*ifft-1,:))
 951        rhor1_pq(2*ifft  ,:) = half*(rhor1(2*ifft  ,:)-rhor1_mq(2*ifft  ,:))
 952        rhor1_mq(2*ifft-1,:) = rhor1_pq(2*ifft-1,:)
 953        rhor1_mq(2*ifft  ,:) =-rhor1_pq(2*ifft  ,:)
 954      end do
 955      rhor1=rhor1_pq
 956      call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 957      call fourdp(cplex,rhog1_mq,rhor1_mq(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 958    end if
 959 
 960    if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
 961 &   dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
 962 
 963 !    calculate \Omega E \cdot P term
 964      if (ipert<=dtset%natom) then
 965 !      phonon perturbation
 966        call  dfptff_ebp(cg,cg1,dtefield,eberry,dtset%mband,mkmem,&
 967 &       mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat)
 968      else if (ipert==dtset%natom+2) then
 969 !      electric field perturbation
 970        call  dfptff_edie(cg,cg1,dtefield,eberry,idir,dtset%mband,mkmem,&
 971 &       mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd)
 972      end if
 973    end if
 974 
 975 !  SPr: don't remove the following comments for debugging
 976 !  call calcdensph(gmet,mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,&
 977 !&   dtset%ntypat,ab_out,dtset%ratsph,rhor1,rprimd,dtset%typat,ucvol,xred,&
 978 !&   idir+1,cplex)
 979 !     write(*,*) ' n ( 1,2)',intgden(1,1),' ',intgden(1,2)
 980 !     write(*,*) ' mx( 1,2)',intgden(2,1),' ',intgden(2,2)
 981 !     write(*,*) ' my( 1,2)',intgden(3,1),' ',intgden(3,2)
 982 !     write(*,*) ' mz( 1,2)',intgden(4,1),' ',intgden(4,2)
 983 !  call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,&
 984 !&     efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,&
 985 !&     enl0,enl1,epaw1,etotal,evar,evdw,exc1,elmag1,ipert,dtset%natom,optene)
 986 !     write(*,*) 'SPr: ek1=',ek1,'  exc1=',exc1,' elmag1=',elmag
 987 !  if (ipert==dtset%natom+5) then
 988 !  !calculate 1st order magnetic potential contribution to the energy
 989 !    call dfpt_e1mag(e1mag,rhor1,rhog1);
 990 !  endif
 991 
 992 !  ######################################################################
 993 !  Skip out of step loop if non-SCF (completed)
 994 !  ----------------------------------------------------------------------
 995 
 996 !  Indeed, nstep loops have been done inside dfpt_vtorho
 997    if (iscf_mod<=0 .and. iscf_mod/=-3) exit
 998 
 999 !  ######################################################################
1000 !  In case of density mixing , compute the total 2nd-order energy,
1001 !  check the exit criterion, then mix the 1st-order density
1002 !  ----------------------------------------------------------------------
1003 
1004    if (iscf_mod>=10) then
1005      optene = 1 ! use double counting scheme
1006      call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,&
1007 &     efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,&
1008 &     enl0,enl1,epaw1,etotal,evar,evdw,exc1,ipert,dtset%natom,optene)
1009      call timab(152,1,tsec)
1010      choice=2
1011      call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
1012 &     etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),&
1013 &     1,iscf_mod,istep,kpt_rbz,maxfor,&
1014 &     mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
1015 &     nstep,occ_rbz,0,prtfor,0,&
1016 &     quit,res2,resid,residm,response,&
1017 &     tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode)
1018      call timab(152,2,tsec)
1019 
1020      if (istep==nstep) quit=1
1021 !    If criteria in scprqt say to quit, then exit the loop over istep
1022      quit_sum=quit
1023      call xmpi_sum(quit_sum,spaceComm,ierr)
1024 
1025      if (quit_sum>0) exit
1026      call status(istep,dtfil%filstat,iexit,level,'call newrho   ')
1027 !    INSERT HERE CALL TO NEWRHO3 : to be implemented
1028      if (psps%usepaw==1) then
1029        MSG_BUG("newrho3 not implemented: use potential mixing!")
1030      end if
1031      initialized=1
1032    end if
1033 
1034 !  ######################################################################
1035 !  Compute the new 1st-order potential from the 1st-order density
1036 !  ----------------------------------------------------------------------
1037 
1038    if (ipert<dtset%natom+10) then
1039      optene=1
1040      call status(istep,dtfil%filstat,iexit,level,'call dfpt_rhotov   ')
1041      call dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,&
1042 &     dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,&
1043 &     nspden,n3xccc,optene,optres,dtset%paral_kgb,dtset%qptn,rhog,rhog1,rhor,rhor1,&
1044 &     rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,nvresid1,res2,vtrial1,vxc,vxc1,xccc3d1,dtset%ixcrot)
1045    end if
1046 
1047 !  ######################################################################
1048 !  In case of potential mixing , compute the total 2nd-order energy,
1049 !  check the exit criterion, then mix the 1st-order potential
1050 !  ----------------------------------------------------------------------
1051 
1052    if (iscf_mod<10) then
1053 
1054 !    PAW: has to compute here the "on-site" 2nd-order energy
1055      if (psps%usepaw==1) then
1056        nzlmopt=0;if (istep==1.and.dtset%pawnzlm>0) nzlmopt=-1
1057        if (istep>1) nzlmopt=dtset%pawnzlm
1058        call paw_an_reset_flags(paw_an1) ! Force the recomputation of on-site potentials
1059        option=2
1060        call pawdenpot(dum,epaw1,epawdc1_dum,ipert,dtset%ixc,my_natom,dtset%natom,dtset%nspden,&
1061 &       psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an1,paw_an,paw_ij1,pawang,dtset%pawprtvol,&
1062 &       pawrad,pawrhoij1,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,&
1063 &       dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp,&
1064 &       mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1065      end if
1066 
1067      optene = 0 ! use direct scheme
1068      call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,&
1069 &     efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,&
1070 !&     enl0,enl1,epaw1,etotal,evar,evdw,exc1,elmag1,ipert,dtset%natom,optene)
1071 !    !debug: compute the d2E/d-qd+q energy, should be equal to the one from previous line
1072 !    if(.not.kramers_deg) then
1073 !      call dfpt_etot(dtset%berryopt,deltae_mq,eberry_mq,edocc_mq,eeig0_mq,eew,efrhar,efrkin,&
1074 !&       efrloc,efrnl,efrx1,efrx2,ehart1_mq,ek0_mq,ek1_mq,eii,elast_mq,eloc0_mq,elpsp1_mq,&
1075 !&       enl0_mq,enl1_mq,epaw1_mq,etotal_mq,evar_mq,evdw,exc1_mq,elmag1_mq,ipert,dtset%natom,optene)
1076 !     end if
1077 &     enl0,enl1,epaw1,etotal,evar,evdw,exc1,ipert,dtset%natom,optene)
1078 
1079      call timab(152,1,tsec)
1080      choice=2
1081      call status(istep,dtfil%filstat,iexit,level,'print info    ')
1082      call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
1083 &     etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),&
1084 &     1,iscf_mod,istep,kpt_rbz,maxfor,&
1085 &     mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
1086 &     nstep,occ_rbz,0,prtfor,0,&
1087 &     quit,res2,resid,residm,response,&
1088 &     tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode)
1089 !     !debug: print the information about residuals at "-q"
1090 !     if(.not.kramers_deg) then
1091 !       call scprqt(choice,cpus,deltae_mq,diffor,dtset,eigen0,&
1092 !&       etotal_mq,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),&
1093 !&       1,iscf_mod,istep,kpt_rbz,maxfor,&
1094 !&       mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
1095 !&       nstep,occ_rbz,0,prtfor,0,&
1096 !&       quit,res2_mq,resid_mq,residm_mq,response,&
1097 !&       tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode)
1098 !     endif
1099      call timab(152,2,tsec)
1100 
1101 !    If criteria in scprqt say to quit, then exit the loop over istep
1102      quit_sum=quit
1103      call xmpi_sum(quit_sum,spaceComm,ierr)
1104      if (quit_sum>0) exit
1105 
1106      ! TODO
1107      ! Better error handling is the SCF cycle goes bananas:
1108      !   Write a BIG warning in the output file and save the wavefunctions.
1109      !   so that we can restart.
1110      if(iscf_mod/=-3)then
1111 !      Note that nvresid1 and vtrial1 are called vresid and vtrial inside this routine
1112        call dfpt_newvtr(cplex,dbl_nnsclo,dielar,dtset,etotal,pawfgr%fintocoa,&
1113 &       initialized,iscf_mod,ispmix,istep,mix,pawfgr%coatofin,&
1114 &       mpi_enreg,my_natom,nfftf,nfftmix,ngfftf,ngfftmix,npawmix,pawrhoij1,&
1115 &       qphon,rhor1,rprimd,psps%usepaw,nvresid1,vtrial1)
1116        if (.not.kramers_deg) then
1117        !same problem as with density reconstruction, TODO proper fft parallelization...
1118          do ifft=1,nfftf
1119            vtrial1_mq(2*ifft-1,1)=+vtrial1(2*ifft-1,1)
1120            vtrial1_mq(2*ifft-1,2)=+vtrial1(2*ifft-1,2)
1121            vtrial1_mq(2*ifft  ,1)=-vtrial1(2*ifft  ,1)
1122            vtrial1_mq(2*ifft  ,2)=-vtrial1(2*ifft  ,2)
1123            vtrial1_mq(2*ifft-1,3)= vtrial1(2*ifft  ,4) !Re[V^12]
1124            vtrial1_mq(2*ifft  ,3)= vtrial1(2*ifft-1,4) !Im[V^12],see definition of v(:,4) cplex=2 case
1125            vtrial1_mq(2*ifft  ,4)= vtrial1(2*ifft-1,3) !Re[V^21]=Re[V^12]
1126            vtrial1_mq(2*ifft-1,4)= vtrial1(2*ifft  ,3) !Re[V^21]=Re[V^12]
1127          end do
1128        end if
1129        initialized=1
1130      end if
1131    end if
1132 
1133 !  ######################################################################
1134 !  END MINIMIZATION ITERATIONS
1135 !  Note that there are different "exit" instructions within the loop
1136 !  ######################################################################
1137  end do ! istep
1138 
1139  ! Avoid pending requests if itime == ntime.
1140  call xmpi_wait(quitsum_request,ierr)
1141  if (timelimit_exit == 1) istep = istep - 1
1142 
1143 !SP : Here read the _DDB file and extract the Born effective charge and
1144 !     dielectric constant.
1145 ! The idea is to supress the divergence due to a residual Born effective charge
1146 ! by renormalizing the v_hart1. For this, the difference between the ionic
1147 ! Z_kappa and the Born effective charge divided by the dielectric constant is used.
1148 ! ---------------------------------------------------------------------------------
1149  if ((dtset%getddb .ne. 0 .or. dtset%irdddb .ne.0) .and. qzero .ne. 1) then
1150    ABI_ALLOCATE(rhor2,(cplex*nfftf,nspden))
1151    ABI_ALLOCATE(resid2,(dtset%mband*nkpt_rbz*nspden))
1152    ABI_ALLOCATE(rhog2,(2,nfftf))
1153    ABI_ALLOCATE(vtrial2,(cplex*nfftf,nspden))
1154 
1155    Z_kappa = nint(psps%ziontypat(dtset%typat(ipert))) ! Charge ionic from the psp
1156    qred2cart = two_pi*gprimd
1157    q_cart = MATMUL(qred2cart,qphon)
1158    q_cart = q_cart/SQRT(dot_product(q_cart,q_cart))
1159    diel_q = dot_product(MATMUL(dielt,q_cart),q_cart)
1160    zeff_bar = SUM(zeff(:,:,:),DIM=3)/dtset%natom
1161    zeff_red = MATMUL(zeff_bar(:,:),rprimd(:,idir))/two_pi
1162    qphon_norm = SQRT(dot_product(qphon,qphon))
1163    q_cart = MATMUL(qred2cart,qphon)
1164    born_bar = dot_product(q_cart,zeff_red(:))
1165    zeff_red = MATMUL(zeff(:,:,ipert),rprimd(:,idir))/two_pi
1166    born = dot_product(q_cart,zeff_red(:))
1167 
1168 ! To avoid problem of divergence (0/0) we add a small value to qphon
1169    qphon2 = qphon + tol6
1170    renorm = (1-(qphon2(idir)*Z_kappa-(born-born_bar)/diel_q)/(qphon2(idir)*Z_kappa-born/diel_q))
1171 
1172    vtrial2(:,1) = vtrial1(:,1) -renorm*vhartr1
1173 
1174    call dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,&
1175 &   dbl_nnsclo,dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,dtset%qptn,edocc,&
1176 &   eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,enl0,enl1,fermie1,gh0c1_set,gh1c_set,&
1177 &   gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,&
1178 &   mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,&
1179 &   nhat1,nkpt_rbz,npwarr,npwar1,res3,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid2,&
1180 &   occkq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,&
1181 &   pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid2,residm2,rhog2,&
1182 &   rhor2,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,usecprj,useylmgr1,ddk_f,&
1183 &   vtrial,vtrial2,wtk_rbz,xred,ylm,ylm1,ylmgr1,1)
1184 
1185    write(msg,'(a)') ' '//NEW_LINE('A')//'&
1186 &   ---------------------------------'
1187    call wrtout(ab_out,msg,'COLL')
1188    write(msg,'(a,a)')'  The charge sum rule is activated'//NEW_LINE('A')//'&
1189 &   ---------------------------------'
1190    call wrtout(ab_out,msg,'COLL')
1191    write(msg,'(a,i4)') ' Z_ion (psp):',Z_kappa
1192    call wrtout(ab_out,msg,'COLL')
1193    write(msg,'(a,f12.8)') ' Residual Born effective charge: ',born
1194    call wrtout(ab_out,msg,'COLL')
1195    write(msg,'(a,f12.8)') ' Renormalisation: ',renorm
1196    call wrtout(ab_out,msg,'COLL')
1197    if (renorm > 0.01 ) then
1198      write(msg,'(a,a)')'   WARNING: The renormalisation seems large (> 0.01).'//NEW_LINE('A')//'&
1199 &     You might consider increasing the k-point grid.'
1200      MSG_WARNING(msg)
1201      call wrtout(ab_out,msg,'COLL')
1202    end if
1203    write(msg,'(a)') ' '
1204    call wrtout(ab_out,msg,'COLL')
1205 
1206    ABI_DEALLOCATE(nvresid2)
1207    ABI_DEALLOCATE(rhor2)
1208    ABI_DEALLOCATE(resid2)
1209    ABI_DEALLOCATE(rhog2)
1210    ABI_DEALLOCATE(vtrial2)
1211  end if
1212 
1213  if (iscf_mod>0.or.iscf_mod==-3)  then
1214    ABI_DEALLOCATE(nvresid1)
1215    if (.not.kramers_deg) then
1216      ABI_DEALLOCATE(nvresid1_mq)
1217    end if
1218  end if
1219 
1220 !######################################################################
1221 !Additional steps after SC iterations
1222 !----------------------------------------------------------------------
1223 
1224  call timab(160,1,tsec)
1225 
1226 !Compute Dynamic magnetic charges (dmc) in case of rfphon,
1227 !and magnetic susceptibility in case of rfmagn from first order density
1228 !(results to be comapred to dmc from d2e)
1229 !SPr deb
1230 !if (ipert<=dtset%natom.and.dtset%nspden>=2) then
1231 !
1232 !  mpi_comm_sphgrid=mpi_enreg%comm_fft
1233 !  call mean_fftr(rhor1(:,1),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid)
1234 !  write(*,*) '   Mean 1st order density: ', mean_rhor1
1235 !  call mean_fftr(rhor1(:,2),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid)
1236 !  if (dtset%nspden==2) then
1237 !    write(*,*) '        1st order m_z    : ', mean_rhor1
1238 !  else !nspden==4
1239 !    write(*,*) '        1st order m_x    : ', mean_rhor1
1240 !    call mean_fftr(rhor1(:,3),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid)
1241 !    write(*,*) '        1st order m_y    : ', mean_rhor1
1242 !    call mean_fftr(rhor1(:,4),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid)
1243 !    write(*,*) '        1st order m_z    : ', mean_rhor1
1244 !  endif
1245 !
1246 !endif
1247 
1248 
1249 !Eventually close the dot file, before calling dfpt_nstdy
1250  if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<=1.0d-7.and.&
1251 & (dtset%berryopt/=4 .and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
1252 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.&
1253 & ipert==dtset%natom+10.or.ipert==dtset%natom+11) then
1254    call wfk_close(ddk_f(1))
1255  end if
1256  if ((ipert==dtset%natom+10 .and. idir>3) .or. ipert==dtset%natom+11) then
1257    call wfk_close(ddk_f(2))
1258  end if
1259  if (ipert==dtset%natom+11) then
1260    call wfk_close(ddk_f(3))
1261    if(idir>3) call wfk_close(ddk_f(4))
1262  end if
1263 
1264 !Deallocate the no more needed arrays
1265  if (iscf_mod>0.and.nstep>0) then
1266    call ab7_mixing_deallocate(mix)
1267  end if
1268  if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then
1269    ABI_DEALLOCATE(dielinv)
1270    ABI_DEALLOCATE(susmat)
1271  end if
1272  if(allocated(rhorfermi))  then
1273    ABI_DEALLOCATE(rhorfermi)
1274  end if
1275  if(allocated(rhorfermi_mq)) then
1276    ABI_DEALLOCATE(rhorfermi_mq)
1277  end if
1278  if(allocated(nhatfermi))  then
1279    ABI_DEALLOCATE(nhatfermi)
1280  end if
1281  if(allocated(pawrhoijfermi))  then
1282    call pawrhoij_free(pawrhoijfermi)
1283    ABI_DATATYPE_DEALLOCATE(pawrhoijfermi)
1284  end if
1285  if(psps%usepaw==1) then
1286    if (mk1mem/=0.and.usecprj==1) then
1287      call pawcprj_free(cprj1)
1288    end if
1289    do iatom=1,my_natom
1290      if (pawfgrtab(iatom)%nhatfr_allocated>0)  then
1291        ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr)
1292      end if
1293      pawfgrtab(iatom)%nhatfr_allocated=0
1294    end do
1295    if (nstep>0.and.iscf_mod>0) then
1296      do iatom=1,my_natom
1297        pawrhoij1(iatom)%lmnmix_sz=0
1298        pawrhoij1(iatom)%use_rhoijres=0
1299        ABI_DEALLOCATE(pawrhoij1(iatom)%kpawmix)
1300        ABI_DEALLOCATE(pawrhoij1(iatom)%rhoijres)
1301      end do
1302    end if
1303  end if ! PAW
1304  ABI_DATATYPE_DEALLOCATE(cprj1)
1305  ABI_DEALLOCATE(nhat1gr)
1306 
1307  call timab(160,2,tsec)
1308  call timab(150,1,tsec)
1309 
1310  if (psps%usepaw==0.and.dtset%userie/=919.and. &
1311 & (ipert==dtset%natom+3.or.ipert==dtset%natom+4)) then
1312    call status(0,dtfil%filstat,iexit,level,'enter dfpt_nselt  ')
1313    call dfpt_nselt(blkflg,cg,cg1,cplex,&
1314 &   d2bbb,d2lo,d2nl,ecut,dtset%ecutsm,dtset%effmass_free,&
1315 &   gmet,gprimd,gsqcut,idir,&
1316 &   ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,dtset%mband,mgfftf,&
1317 &   mkmem,mk1mem,mpert,mpi_enreg,psps%mpsang,mpw,mpw1,&
1318 &   dtset%natom,nband_rbz,nfftf,ngfftf,&
1319 &   nkpt_rbz,nkxc,dtset%nloalg,&
1320 &   npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,&
1321 &   nsym1,dtset%ntypat,occ_rbz,&
1322 &   dtset%paral_kgb,ph1d,dtset%prtbbb,psps,dtset%qptn,rhog,&
1323 &   rhor,rhor1,rmet,rprimd,symrc1,dtset%typat,ucvol,&
1324 &   wtk_rbz,xred,ylm,ylm1,ylmgr,ylmgr1)
1325  end if
1326 
1327 !Use of NSTPAW3 for NCPP (instead of DFPT_NSELT/DFPT_NSTDY) can be forced with userie=919
1328 !MT oct. 2015: this works perfectly on all automatic tests
1329  if(ipert<=dtset%natom+4)then
1330    if (psps%usepaw==1.or.dtset%userie==919) then
1331      call status(0,dtfil%filstat,iexit,level,'enter dfpt_nstpaw ')
1332      call dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,&
1333 &     eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,&
1334 &     kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,ncpgr,&
1335 &     nfftf,ngfftf,nhat,nhat1,nkpt_rbz,nkxc,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,&
1336 &     nsym1,n3xccc,occkq,occ_rbz,paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,&
1337 &     pawrhoij,pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,&
1338 &     symrl1,ucvol,usecprj,psps%usepaw,usexcnhat,useylmgr1,vhartr1,vpsp1,vtrial,vtrial1,vxc,wtk_rbz,&
1339 &     xccc3d1,xred,ylm,ylm1,ylmgr1)
1340    else
1341      call status(0,dtfil%filstat,iexit,level,'enter dfpt_nstdy  ')
1342      if (dtset%nspden==4) then
1343        call dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,gmet,&
1344 &       gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,mpert,mpi_enreg,&
1345 &       mpw,mpw1,nattyp,nband_rbz,nfftf,ngfftf,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,&
1346 &       dtset%nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,symrc1,ucvol,&
1347 &       wtk_rbz,xred,ylm,ylm1,rhor=rhor,vxc=vxc)
1348      else
1349        call dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,gmet,&
1350 &       gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,mpert,mpi_enreg,&
1351 &       mpw,mpw1,nattyp,nband_rbz,nfftf,ngfftf,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,&
1352 &       dtset%nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,symrc1,ucvol,&
1353 &       wtk_rbz,xred,ylm,ylm1)
1354      end if
1355    end if
1356  end if
1357 
1358  call timab(150,2,tsec)
1359  call timab(160,1,tsec)
1360 
1361 !calculate Born effective charge and store it in d2lo
1362  if ((dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1363 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17).and.&
1364 & ipert<=dtset%natom) then
1365    call dfptff_bec(cg,cg1,dtefield,dtset%natom,d2lo,idir,ipert,dtset%mband,mkmem,&
1366 &   mpw,mpw1,mpert,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd)
1367    blkflg(:,dtset%natom+2,:,1:dtset%natom)=1
1368  end if
1369 
1370 !calculate dielectric tensor and store it in d2lo
1371  if ((dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1372 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17).and.&
1373 & ipert==dtset%natom+2) then
1374    call dfptff_die(cg,cg1,dtefield,d2lo,idir,ipert,dtset%mband,mkmem,&
1375 &   mpw,mpw1,mpert,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd)
1376    blkflg(:,dtset%natom+2,:,dtset%natom+2)=1
1377  end if
1378 
1379 !If SCF convergence was not reached (for nstep>0),
1380 !print a warning to the output file (non-dummy arguments: nstep,
1381 !residm, diffor - infos from tollist have been saved inside )
1382 !Set also the value of conv_retcode
1383  choice=3
1384  call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
1385 & etotal,favg,fcart,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),&
1386 & 1,iscf_mod,istep,kpt_rbz,maxfor,&
1387 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
1388 & nstep,occ_rbz,0,prtfor,0,&
1389 & quit,res2,resid,residm,response,&
1390 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode)
1391 
1392 !Update the content of the header (evolving variables)
1393  bantot_rbz = sum(nband_rbz(1:nkpt_rbz*dtset%nsppol))
1394  call hdr_update(hdr,bantot_rbz,etotal,fermie,&
1395 & residm,rprimd,occ_rbz,pawrhoij1,xred,dtset%amu_orig(:,1),&
1396 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
1397 
1398 !Optionally provide output of charge density and/or potential in real space,
1399 !as well as analysis of geometrical factors (bond lengths and bond angles).
1400 !Warnings :
1401 !- core charge is excluded from the charge density;
1402 !- the potential is the INPUT vtrial.
1403 
1404  if(ipert==dtset%natom+5.or.ipert<=dtset%natom)then
1405    prtopt=1
1406    if(ipert==dtset%natom+5) then
1407      prtopt=idir+1;
1408      call calcdensph(gmet,mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,&
1409 &     dtset%ntypat,ab_out,dtset%ratsph,rhor1,rprimd,dtset%typat,ucvol,xred,&
1410 &     prtopt,cplex,intgden=intgden,dentot=dentot)
1411      !debug: write out the vtk first-order density components
1412 !    call appdig(pertcase,dtfil%fnameabo_den,fi1o_vtk)
1413 !    call printmagvtk(mpi_enreg,cplex,nspden,nfftf,ngfftf,rhor1,rprimd,adjustl(adjustr(fi1o_vtk)//"_PQ"))
1414 !    call printmagvtk(mpi_enreg,cplex,nspden,nfftf,ngfftf,rhor1,rprimd,adjustl(adjustr(fi1o_vtk)//"_MQ"))
1415      !SPr: add calculation of the contributions to susceptibility from all attomic spheres
1416    end if
1417  end if
1418 
1419  if (iwrite_fftdatar(mpi_enreg)) then
1420    if (dtset%prtden>0) then
1421      rdwrpaw=0
1422      call appdig(pertcase,dtfil%fnameabo_den,fi1o)
1423      ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij
1424      call fftdatar_write_from_hdr("first_order_density",fi1o,dtset%iomode,hdr,&
1425      ngfftf,cplex,nfftf,dtset%nspden,rhor1,mpi_enreg)
1426    end if
1427 
1428    ! first order potentials are always written because the eph code requires them
1429    ! the files are small (much much smaller that 1WFK, actually we should avoid writing 1WFK)
1430    rdwrpaw=0
1431    call appdig(pertcase,dtfil%fnameabo_pot,fi1o)
1432    ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij
1433    call fftdatar_write_from_hdr("first_order_potential",fi1o,dtset%iomode,hdr,&
1434    ngfftf,cplex,nfftf,dtset%nspden,vtrial1,mpi_enreg)
1435 
1436 ! output files for perturbed potential components: vhartr1,vpsp1,vxc
1437 ! NB: only 1 spin for these
1438    if (dtset%prtvha > 0) then
1439      rdwrpaw=0
1440      ABI_ALLOCATE(vhartr1_tmp, (cplex*nfftf, dtset%nspden))
1441      vhartr1_tmp = zero
1442      vhartr1_tmp(:,1) = vhartr1(:)
1443      call appdig(pertcase,dtfil%fnameabo_vha,fi1o)
1444      ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij
1445      call fftdatar_write_from_hdr("first_order_vhartree",fi1o,dtset%iomode,hdr,&
1446      ngfftf,cplex,nfftf,dtset%nspden,vhartr1_tmp,mpi_enreg)
1447      ABI_DEALLOCATE(vhartr1_tmp)
1448    end if
1449 
1450 
1451 ! vpsp1 needs to be copied to a temp array - intent(inout) in fftdatar_write_from_hdr though I do not know why
1452 !   if (dtset%prtvpsp > 0) then
1453 !     rdwrpaw=0
1454 !     call appdig(pertcase,dtfil%fnameabo_vpsp,fi1o)
1455 !     ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij
1456 !     call fftdatar_write_from_hdr("first_order_vpsp",fi1o,dtset%iomode,hdr,&
1457 !       ngfftf,cplex,nfftf,1,vpsp1,mpi_enreg)
1458 !   end if
1459 
1460    if (dtset%prtvxc > 0) then
1461      rdwrpaw=0
1462      call appdig(pertcase,dtfil%fnameabo_vxc,fi1o)
1463      ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij
1464      call fftdatar_write_from_hdr("first_order_vxc",fi1o,dtset%iomode,hdr,&
1465      ngfftf,cplex,nfftf,dtset%nspden,vxc1,mpi_enreg)
1466    end if
1467 
1468 
1469    ! Add rhog1(G=0) to file
1470    if (mpi_enreg%me_g0 == 1) then
1471      if (dtset%iomode == IO_MODE_ETSF) then
1472 #ifdef HAVE_NETCDF
1473        NCF_CHECK(nctk_open_modify(ncid, nctk_ncify(fi1o), xmpi_comm_self))
1474        ncerr = nctk_def_one_array(ncid, nctkarr_t('rhog1_g0', "dp", "two"), varid=varid)
1475        NCF_CHECK(ncerr)
1476        NCF_CHECK(nctk_set_datamode(ncid))
1477        NCF_CHECK(nf90_put_var(ncid, varid, rhog1(:,1)))
1478        NCF_CHECK(nf90_close(ncid))
1479 #endif
1480      else
1481        ! Handle Fortran files.
1482        if (open_file(fi1o, msg, newunit=ncid, form='unformatted', status='old', action="readwrite") /= 0) then
1483          MSG_ERROR(msg)
1484        end if
1485        if (fort_denpot_skip(ncid, msg) /= 0) MSG_ERROR(msg)
1486        write(ncid) rhog1(:,1)
1487        close(ncid)
1488      end if
1489    end if
1490 
1491  end if ! iwrite_fftdatar(mpi_enreg)
1492 
1493 !All procs waiting here...
1494  if(mpi_enreg%paral_kgb==1)then
1495    call timab(61,1,tsec)
1496    call xmpi_barrier(spaceComm)
1497    call timab(61,2,tsec)
1498  end if
1499 
1500 !Deallocate arrays
1501  ABI_DEALLOCATE(fcart)
1502  ABI_DEALLOCATE(vtrial1)
1503  if (.not.kramers_deg) then
1504    ABI_DEALLOCATE(vhartr1_mq)
1505    ABI_DEALLOCATE(vxc1_mq)
1506    ABI_DEALLOCATE(vtrial1_pq)
1507    ABI_DEALLOCATE(vtrial1_mq)
1508  end if
1509  ABI_DEALLOCATE(vhartr1)
1510  ABI_DEALLOCATE(vxc1)
1511  ABI_DEALLOCATE(pwindall)
1512  ABI_DEALLOCATE(qmat)
1513  if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1514 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
1515    call destroy_efield(dtefield)
1516    if(allocated(mpi_enreg%kpt_loc2ibz_sp))  then
1517      ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp)
1518    end if
1519  end if
1520 
1521  if(psps%usepaw==1) then
1522    call paw_an_free(paw_an1)
1523    call paw_ij_free(paw_ij1)
1524  end if
1525  ABI_DATATYPE_DEALLOCATE(paw_an1)
1526  ABI_DATATYPE_DEALLOCATE(paw_ij1)
1527  ABI_DEALLOCATE(nhat1)
1528 
1529  call status(0,dtfil%filstat,iexit,level,'exit')
1530 
1531  call timab(160,2,tsec)
1532  call timab(120,2,tsec)
1533 
1534  DBG_EXIT("COLL")
1535 
1536 end subroutine dfpt_scfcv

ABINIT/dfpt_wfkfermi [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_wfkfermi

FUNCTION

 This routine computes the partial Fermi-level density at a given k-point,
 and the fixed contribution to the 1st-order Fermi energy (nonlocal and kinetic)

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
  cplex=1 if rhoaug is real, 2 if rhoaug is complex
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
              projected with non-local projectors: cprj=<p_i|Cnk>
  cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q
              projected with non-local projectors: cprjq=<p_i|Cnk+q>
  dtfil <type(datafiles_type)>=variables related to files
  eig0_k(nband_k)=GS eigenvalues at k (hartree)
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  ibg=shift to be applied on the location of data in the array cprj
  ibgq=shift to be applied on the location of data in the array cprjq
  icg=shift to be applied on the location of data in the array cg
  icgq=shift to be applied on the location of data in the array cgq
  idir=direction of the current perturbation
  ikpt=number of the k-point
  ipert=type of the perturbation
  isppol=1 for unpolarized, 2 for spin-polarized
  kptopt=option for the generation of k points
  mband=maximum number of bands
  mcgq=second dimension of the cgq array
  mcprjq=second dimension of the cprjq array
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw or wfs at k
  nband_k=number of bands at this k point for that spin polarization
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  npw_k=number of plane waves at this k point
  npw1_k=number of plane waves at this k+q point
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  occ_k(nband_k)=occupation number for each band (usually 2) for each k.
  prtvol=control print volume and debugging output
  rf_hamkq <type(gs_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
  rhoaug(cplex*n4,n5,n6)= density in electrons/bohr**3,
   on the augmented fft grid. (cumulative, so input as well as output)
  rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
   if this ratio has been attributed to the band n (second argument), zero otherwise
  wtk_k=weight assigned to the k point.

OUTPUT

  eig1_k(2*nband_k**2)=first-order eigenvalues (hartree)
  fe1fixed_k(nband_k)=contribution to 1st-order Fermi energy
      from changes of occupation from all bands at this k point.
  fe1norm_k(nband_k)=contribution to normalization for above
  rhoaug(cplex*n4,n5,n6)= Fermi-level density in electrons/bohr**3,
   on the augmented fft grid. (cumulative, so input as well as output).
  ==== if (gs_hamkq%usepaw==1) ====
    pawrhoijfermi(natom) <type(pawrhoij_type)>= paw rhoij occupancies
       at Fermi level (cumulative, so input as well as output)

PARENTS

      dfpt_rhofermi

CHILDREN

      dfpt_accrho,dotprod_g,getgh1c,pawcprj_alloc,pawcprj_axpby,pawcprj_copy
      pawcprj_free,pawcprj_get,status,timab,wrtout

SOURCE

4368 subroutine dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,&
4369 &          dtfil,eig0_k,eig1_k,fe1fixed_k,fe1norm_k,gs_hamkq,&
4370 &          ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,&
4371 &          kptopt,mband,mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,&
4372 &          npw_k,npw1_k,nspinor,nsppol,occ_k,pawrhoijfermi,prtvol,&
4373 &          rf_hamkq,rhoaug,rocceig,wtk_k)
4374 
4375  use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type
4376  use m_getgh1c,     only : getgh1c
4377  use m_dfpt_mkrho,  only : dfpt_accrho
4378 
4379 !This section has been created automatically by the script Abilint (TD).
4380 !Do not modify the following lines by hand.
4381 #undef ABI_FUNC
4382 #define ABI_FUNC 'dfpt_wfkfermi'
4383 !End of the abilint section
4384 
4385  implicit none
4386 
4387 !Arguments ------------------------------------
4388 !scalars
4389  integer,intent(in) :: cplex,ibg,ibgq,icg,icgq,idir,ikpt
4390  integer,intent(in) :: ipert,isppol,kptopt,mband,mcgq,mcprjq,mkmem,mpw,ncpgr
4391  integer,intent(in) :: npw1_k,nspinor,nsppol,prtvol
4392  integer,intent(inout) :: nband_k,npw_k
4393  real(dp),intent(in) :: wtk_k
4394  type(MPI_type),intent(in) :: mpi_enreg
4395  type(datafiles_type),intent(in) :: dtfil
4396  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
4397  type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
4398 !arrays
4399  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),cgq(2,mcgq)
4400  real(dp),intent(in) :: eig0_k(nband_k),occ_k(nband_k),rocceig(nband_k,nband_k)
4401  real(dp),intent(inout) :: rhoaug(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc)
4402  real(dp),intent(inout) :: eig1_k(2*nband_k**2)
4403  real(dp),intent(out) :: fe1fixed_k(nband_k)
4404  real(dp),intent(out) :: fe1norm_k(nband_k)
4405  type(pawcprj_type),intent(in) :: cprj(gs_hamkq%natom,nspinor*mband*mkmem*nsppol*gs_hamkq%usecprj)
4406  type(pawcprj_type),intent(in) :: cprjq(gs_hamkq%natom,mcprjq)
4407  type(pawrhoij_type),intent(inout) :: pawrhoijfermi(gs_hamkq%natom*gs_hamkq%usepaw)
4408 
4409 !Local variables-------------------------------
4410 !scalars
4411  integer,parameter :: level=18
4412  integer :: berryopt,counter,iband,iexit,ii,indx,iorder_cprj
4413  integer :: ipw,me,nkpt_max,optlocal,optnl,opt_accrho,opt_corr
4414  integer :: opt_gvnl1,sij_opt,tim_fourwf,tim_getgh1c,usevnl
4415  real(dp) :: dotr,lambda,wtband
4416  character(len=500) :: msg
4417 !arrays
4418  real(dp) :: dum_grad_berry(1,1),dum_gvnl1(1,1),dum_gs1(1,1),tsec(2)
4419  real(dp),allocatable :: cwave0(:,:),cwaveq(:,:),gh1(:,:)
4420  type(pawcprj_type),allocatable :: cwaveprj0(:,:),cwaveprjq(:,:),cwaveprj_tmp(:,:)
4421 
4422 ! *********************************************************************
4423 
4424  DBG_ENTER('COLL')
4425 
4426 !Check arguments validity
4427  if (ipert>gs_hamkq%natom.and.ipert/=gs_hamkq%natom+3.and.ipert/=gs_hamkq%natom+4.and.ipert/=gs_hamkq%natom+5) then !SPr rfmagn deb
4428    msg='  wrong ipert argument !'
4429    MSG_BUG(msg)
4430  end if
4431  if (cplex/=1) then
4432    MSG_BUG('wrong cplex/=1 argument !')
4433  end if
4434 
4435 !Debugging statements
4436  call status(0,dtfil%filstat,iexit,level,'enter dfpt_wfkfermi')
4437  if(prtvol==-level)then
4438    write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,'dfpt_wfkfermi : enter'
4439    call wrtout(std_out,msg,'PERS')
4440  end if
4441  nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1
4442 
4443  if(prtvol>2 .or. ikpt<=nkpt_max)then
4444    write(msg, '(a,a,i5,2x,a,3f9.5,2x,a)' ) ch10,&
4445 &   ' Non-SCF iterations; k pt #',ikpt,'k=',gs_hamkq%kpt_k(:),' band residuals:'
4446    call wrtout(std_out,msg,'PERS')
4447  end if
4448 
4449 !Retrieve parallelism data
4450  me=mpi_enreg%me_kpt
4451 !Initializations and allocations
4452 
4453  ABI_ALLOCATE(gh1,(2,npw1_k*nspinor))
4454  ABI_ALLOCATE(cwave0,(2,npw_k*nspinor))
4455  ABI_ALLOCATE(cwaveq,(2,npw1_k*nspinor))
4456  iorder_cprj=0 ; eig1_k(:)=zero
4457  if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
4458    ABI_DATATYPE_ALLOCATE(cwaveprj0,(gs_hamkq%natom,nspinor))
4459    ABI_DATATYPE_ALLOCATE(cwaveprjq,(gs_hamkq%natom,nspinor))
4460    call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj)
4461    call pawcprj_alloc(cwaveprjq,0,gs_hamkq%dimcprj)
4462  else
4463    ABI_DATATYPE_ALLOCATE(cwaveprj0,(0,0))
4464    ABI_DATATYPE_ALLOCATE(cwaveprjq,(0,0))
4465  end if
4466 !Arguments of getgh1c routine (want only (NL+kin) frozen H(1))
4467  berryopt=0;usevnl=0;sij_opt=-gs_hamkq%usepaw;tim_getgh1c=3
4468  optlocal=0;optnl=1;opt_gvnl1=0
4469  if(ipert==gs_hamkq%natom+5) optnl=0;    ! no 1st order NL in H(1), also no kin, but this will be taken into account later
4470 !if(ipert==gs_hamkq%natom+5) optlocal=0; ! 1st order LOCAL potential present
4471 
4472 !Arguments of the dfpt_accrho routine
4473  tim_fourwf=5 ; opt_accrho=1 ; opt_corr=0
4474 !Null potentially unassigned output variables
4475  fe1fixed_k(:)=zero; fe1norm_k(:)=zero
4476 
4477 !Read the npw and kg records of wf files
4478  call status(0,dtfil%filstat,iexit,level,'before WffRead')
4479 
4480  call timab(139,1,tsec)
4481 
4482 
4483 !Loop over bands
4484  do iband=1,nband_k
4485    counter=100*iband+1
4486 
4487 !  Skip bands not treated by current proc
4488    if(mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me) cycle
4489 
4490    if(prtvol>=10)then
4491      call status(counter,dtfil%filstat,iexit,level,'loop iband    ')
4492    end if
4493 
4494 !  Select occupied bands
4495    if(abs(occ_k(iband))>tol8.and.abs(rocceig(iband,iband))>tol8)then
4496 
4497      wtband=rocceig(iband,iband)/occ_k(iband)
4498 !    Get ground-state wavefunctions at k
4499      do ipw=1,npw_k*nspinor
4500        cwave0(1,ipw)=cg(1,ipw+(iband-1)*npw_k*nspinor+icg)
4501        cwave0(2,ipw)=cg(2,ipw+(iband-1)*npw_k*nspinor+icg)
4502      end do
4503 
4504      if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
4505 !      Read PAW ground state projected WF (cprj)
4506        call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,gs_hamkq%natom,iband,ibg,ikpt,iorder_cprj,&
4507 &       isppol,mband,mkmem,gs_hamkq%natom,1,nband_k,nspinor,nsppol,dtfil%unpaw,&
4508 &       mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,&
4509 &       icpgr=idir,ncpgr=ncpgr)
4510      end if
4511 
4512 !    Read ground-state wavefunctions at k+q
4513      indx=npw1_k*nspinor*(iband-1)+icgq
4514      cwaveq(:,1:npw_k*nspinor)=wtband*cgq(:,1+indx:npw_k*nspinor+indx)
4515      if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
4516 !      Read PAW ground state projected WF (cprj)
4517        indx=nspinor*(iband-1)+ibgq
4518        call pawcprj_copy(cprjq(:,1+indx:nspinor+indx),cwaveprjq)
4519        call pawcprj_axpby(zero,wtband,cwaveprj_tmp,cwaveprjq)
4520      end if
4521 
4522      if(prtvol>=10)then
4523        call status(0,dtfil%filstat,iexit,level,'after wf read ')
4524      end if
4525 
4526 
4527 !    Apply H^(1)-Esp.S^(1) to Psi^(0) (H(^1)=only (NL+kin) frozen part)
4528      lambda=eig0_k(iband)
4529      call getgh1c(berryopt,cwave0,cwaveprj0,gh1,dum_grad_berry,dum_gs1,gs_hamkq,dum_gvnl1,&
4530 &     idir,ipert,lambda,mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamkq,sij_opt,&
4531 &     tim_getgh1c,usevnl)
4532 !    Compute Eig1=<Psi^(0)|H^(1)-Eps.S^(1)|Psi(0)>
4533      call dotprod_g(dotr,lambda,gs_hamkq%istwf_k,npw_k*nspinor,1,cwave0,gh1,mpi_enreg%me_g0, &
4534 &     mpi_enreg%comm_spinorfft)
4535      indx=2*iband-1+(iband-1)*2*nband_k
4536      eig1_k(indx)=dotr
4537 !    Compute the fixed contribution to the 1st-order Fermi energy
4538      fe1fixed_k(iband)=two*wtband*eig1_k(indx)
4539      fe1norm_k(iband) =two*wtband
4540 
4541 !    Accumulate contribution to density and PAW occupation matrix
4542 
4543      call dfpt_accrho(counter,cplex,cwave0,cwaveq,cwaveq,cwaveprj0,cwaveprjq,dotr,&
4544 &     dtfil%filstat,gs_hamkq,iband,0,0,isppol,kptopt,mpi_enreg,gs_hamkq%natom,nband_k,ncpgr,&
4545 &     npw_k,npw1_k,nspinor,occ_k,opt_accrho,pawrhoijfermi,prtvol,rhoaug,tim_fourwf,&
4546 &     opt_corr,wtk_k)
4547 
4548    end if ! End of non-zero occupation and rocceig
4549 
4550  end do ! End loop over bands
4551 
4552  call timab(139,2,tsec)
4553  call timab(130,1,tsec)
4554  call status(0,dtfil%filstat,iexit,level,'after loops   ')
4555 
4556  ABI_DEALLOCATE(cwave0)
4557  ABI_DEALLOCATE(cwaveq)
4558  ABI_DEALLOCATE(gh1)
4559  if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
4560    call pawcprj_free(cwaveprj0)
4561    call pawcprj_free(cwaveprjq)
4562  end if
4563  ABI_DATATYPE_DEALLOCATE(cwaveprj0)
4564  ABI_DATATYPE_DEALLOCATE(cwaveprjq)
4565 
4566 !Structured debugging : if prtvol=-level, stop here.
4567  if(prtvol==-level)then
4568    write(msg,'(a,a1,a,i2,a)')' fermie3 : exit prtvol=-',level,', debugging mode => stop '
4569    MSG_ERROR(msg)
4570  end if
4571 
4572  call status(0,dtfil%filstat,iexit,level,'exit dfpt_wfkfermi')
4573 
4574  call timab(130,2,tsec)
4575 
4576  DBG_EXIT('COLL')
4577 
4578 end subroutine dfpt_wfkfermi

ABINIT/m_dfpt_scfcv [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_scfcv

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XG, DRH, MB, XW, MT, SPr, XW, MV, MM, AR)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_dfpt_scfcv
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_ab7_mixing
33  use m_efield
34  use m_errors
35  use m_abicore
36  use m_wfk
37  use m_xmpi
38  use m_nctk
39  use m_hdr
40 #ifdef HAVE_NETCDF
41  use netcdf
42 #endif
43 
44  use m_dtfil,    only : status
45  use m_cgtools,  only : mean_fftr, overlap_g, dotprod_vn, dotprod_vn, dotprod_g
46  use m_fstrings, only : int2char4, sjoin
47  use m_geometry, only : metric, stresssym
48  use m_time,     only : abi_wtime, sec2str, timab
49  use m_io_tools, only : open_file, file_exists, get_unit, iomode_from_fname
50  use m_exit,     only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
51  use m_mpinfo,   only : iwrite_fftdatar, proc_distrb_cycle
52  use m_kg,       only : getcut, mkkin, kpgstr, mkkpg
53  use m_fft,      only : fftpac, fourdp
54  use m_ioarr,    only : ioarr, fftdatar_write_from_hdr, fort_denpot_skip
55  use m_pawang,   only : pawang_type
56  use m_pawrad,   only : pawrad_type
57  use m_pawtab,   only : pawtab_type
58  use m_paw_an,   only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
59  use m_paw_ij,   only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
60  use m_pawfgrtab,only : pawfgrtab_type
61  use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_io, pawrhoij_get_nspden
62  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_copy, pawcprj_axpby, pawcprj_free, pawcprj_getdim
63  use m_pawdij,   only : pawdij, pawdijfr, symdij
64  use m_pawfgr,   only : pawfgr_type
65  use m_paw_denpot,  only : pawdenpot
66  use m_paw_dfpt,    only : pawdfptenergy
67  use m_paw_nhat,    only : pawmknhat,pawnhatfr
68  use m_rf2,         only : rf2_getidirs
69  use m_dens,        only : calcdensph
70  use m_dfpt_fef,    only : dfptff_initberry, qmatrix, dfptff_edie, dfptff_ebp, dfptff_die, dfptff_bec
71  use m_dfpt_vtorho, only : dfpt_vtorho
72  use m_paral_atom,  only : get_my_atmtab, free_my_atmtab
73  use m_common,      only : scprqt
74  use m_prcref,      only : moddiel
75  use m_dfpt_rhotov, only : dfpt_rhotov
76  use m_dfpt_mkvxc,    only : dfpt_mkvxc, dfpt_mkvxc_noncoll
77  use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
78  use m_mklocl,     only : dfpt_vlocal, vlocalstr
79  use m_dfpt_nstwf,   only : dfpt_nstpaw, dfpt_nstwf
80 
81  implicit none
82 
83  private

ABINIT/newfermie1 [ Functions ]

[ Top ] [ Functions ]

NAME

 newfermie1

FUNCTION

 This routine computes the derivative of the fermi energy wrt
 the active perturbation for use in evaluating the edocc term
 and active subspace contribution to the first-order wavefunctions
 in the case of metals. This is presently used only for the
 strain and magnetic field perturbations, and only for Q = 0.

INPUTS

  cplex: if 1, real space 1-order functions on FFT grid are REAL,
    if 2, COMPLEX
  fe1fixed=fixed contribution to the first-order Fermi energy
  ipert=index of perturbation
  istep=index of the number of steps in the routine scfcv
  ixc= choice of exchange-correlation scheme
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms
  nfft=(effective) number of FFT grid points (for this processor)
  nfftot= total number of FFT grid points
  nhatfermi(nfft,nspden)=fermi-level compensation charge density (PAW only)
  nspden=number of spin-density components
  ntypat=number of atom types
  occopt=option for occupancies
  paw_an(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
  paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh
  paw_ij1(natom) <type(paw_ij_type)>=(1st-order) paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawnzlm=-- PAW only -- option for the computation of non-zero
          lm moments of the on-sites densities
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij1(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies
  pawrhoijfermi(natom) <type(pawrhoij_type)>=paw rhoij occupancies at Fermi level
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  xclevel= XC functional level
  prtvol=control print volume and debugging output
  rhorfermi(nfft,nspden)=fermi-level electronic density
  ucvol=unit cell volume in bohr**3
  usepaw=1 if PAW is activated
  usexcnhat= -PAW only- flag controling use of compensation density in Vxc
  vtrial1(cplex*nfft,nspden)=1-st order potential
  vxc1(cplex*nfft,nspden)=1-st order XC potential

OUTPUT

  (see side effects)

SIDE EFFECTS

  fermie1=derivative of fermi energy wrt perturbation
   at input  : old value
   at output : updated value

PARENTS

      dfpt_scfcv

CHILDREN

      dotprod_vn,free_my_atmtab,get_my_atmtab,pawdfptenergy,wrtout

SOURCE

1746 subroutine newfermie1(cplex,fermie1,fe1fixed,ipert,istep,ixc,my_natom,natom,nfft,nfftot,&
1747 &                     nhatfermi,nspden,ntypat,occopt,paw_an,paw_an1,paw_ij1,pawang,pawnzlm,pawrad,&
1748 &                     pawrhoij1,pawrhoijfermi,pawtab,pawxcdev,prtvol,rhorfermi,&
1749 &                     ucvol,usepaw,usexcnhat,vtrial1,vxc1,xclevel,&
1750 &                     mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1751 
1752 
1753 !This section has been created automatically by the script Abilint (TD).
1754 !Do not modify the following lines by hand.
1755 #undef ABI_FUNC
1756 #define ABI_FUNC 'newfermie1'
1757 !End of the abilint section
1758 
1759  implicit none
1760 
1761 !Arguments -------------------------------
1762 !scalars
1763  integer,intent(in) :: cplex,ipert,istep,ixc,my_natom,natom,nfft,nfftot,nspden,ntypat
1764  integer,intent(in) :: occopt,pawnzlm,pawxcdev,prtvol,usepaw,usexcnhat,xclevel
1765  integer,optional,intent(in) :: comm_atom
1766  real(dp),intent(in) :: fe1fixed,ucvol
1767  real(dp),intent(inout) :: fermie1
1768  type(pawang_type),intent(in) :: pawang
1769 !arrays
1770  integer,optional,target,intent(in) :: mpi_atmtab(:)
1771  real(dp),intent(in) :: rhorfermi(nfft,nspden),vtrial1(cplex*nfft,nspden)
1772  real(dp),intent(in) :: nhatfermi(:,:),vxc1(:,:)
1773  type(paw_an_type),intent(in) :: paw_an(my_natom*usepaw)
1774  type(paw_an_type),intent(inout) :: paw_an1(my_natom*usepaw)
1775  type(paw_ij_type),intent(inout) :: paw_ij1(my_natom*usepaw)
1776  type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw)
1777  type(pawrhoij_type),intent(in) :: pawrhoij1(my_natom*usepaw),pawrhoijfermi(my_natom*usepaw)
1778  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
1779 
1780 !Local variables-------------------------------
1781 !scalars
1782  integer :: ipert0,my_comm_atom,nzlmopt,nzlmopt_fermi,option,pawprtvol
1783  logical :: my_atmtab_allocated,paral_atom
1784  real(dp) :: doti,fe1_scf,fe1_tmp,fermie1_new,fermie1rs
1785  character(len=500) :: msg
1786 !arrays
1787  integer, pointer :: my_atmtab(:)
1788  real(dp) :: fe1_paw(2)
1789  real(dp), allocatable :: rhor_nonhat(:,:),vtrial1_novxc(:,:)
1790 
1791 ! *********************************************************************
1792 
1793 !Tests
1794  if (cplex==2) then
1795    msg='Not compatible with cplex=2!'
1796    MSG_BUG(msg)
1797  end if
1798  if (usepaw==1.and.usexcnhat==0.and.(size(nhatfermi)<=0.or.size(vxc1)<=0)) then
1799    msg='Should have nhatfermi and vxc1 allocated with usexcnhat=0!'
1800    MSG_BUG(msg)
1801  end if
1802 
1803 !Set up parallelism over atoms
1804  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1805  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1806  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1807  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1808 
1809  if(occopt>=3 .and. occopt <=8) then
1810 
1811 !  The product of the current trial potential and the so-called Fermi level
1812 !  density is integrated to give the local potential contributions to the
1813 !  first-order Fermi level.
1814    option=1
1815    if (usepaw==1.and.usexcnhat==0) then
1816      ABI_ALLOCATE(rhor_nonhat,(nfft,nspden))
1817      ABI_ALLOCATE(vtrial1_novxc,(nfft,nspden))
1818      rhor_nonhat(1:nfft,1:nspden)=rhorfermi(1:nfft,1:nspden)-nhatfermi(1:nfft,1:nspden)
1819      vtrial1_novxc(1:nfft,1:nspden)=vtrial1(1:nfft,1:nspden)-vxc1(1:nfft,1:nspden)
1820      call dotprod_vn(cplex,rhor_nonhat,fe1_scf,doti,nfft,nfftot,&
1821 &     nspden,option,vtrial1,ucvol)
1822      call dotprod_vn(cplex,nhatfermi,fe1_tmp,doti,nfft,nfftot,&
1823 &     nspden,option,vtrial1_novxc,ucvol)
1824      fe1_scf=fe1_scf+fe1_tmp
1825      ABI_DEALLOCATE(rhor_nonhat)
1826      ABI_DEALLOCATE(vtrial1_novxc)
1827    else
1828      call dotprod_vn(cplex,rhorfermi,fe1_scf,doti,nfft,nfftot,&
1829 &     nspden,option,vtrial1,ucvol)
1830    end if
1831 
1832    fe1_paw(:)=zero
1833 !  PAW on-site contribution (use Fermi level occupation matrix)
1834    if (usepaw==1) then
1835      ipert0=0;pawprtvol=0
1836      nzlmopt=0;if (istep>1) nzlmopt=pawnzlm
1837      if (istep==1.and.pawnzlm>0) nzlmopt=-1
1838      nzlmopt_fermi=0;if (pawnzlm>0) nzlmopt_fermi=-1
1839      call pawdfptenergy(fe1_paw,ipert,ipert0,ixc,my_natom,natom,ntypat,nzlmopt,&
1840 &     nzlmopt_fermi,paw_an,paw_an1,paw_ij1,pawang,pawprtvol,pawrad,&
1841 &     pawrhoij1,pawrhoijfermi,pawtab,pawxcdev,xclevel,&
1842 &     mpi_atmtab=my_atmtab, comm_atom=my_comm_atom)
1843    end if
1844 
1845 !  The fixed contributions consisting of non-local potential and kinetic terms
1846 !  are added
1847    fermie1_new=fe1fixed+fe1_scf+fe1_paw(1)
1848    fermie1rs=(fermie1-fermie1_new)**2
1849    fermie1=fermie1_new
1850 
1851    if(prtvol>=10)then
1852      write(msg, '(a,i5,2es18.8)' ) ' fermie1, residual squared',istep,fermie1,fermie1rs
1853      call wrtout(std_out,msg,'COLL')
1854    end if
1855 
1856  else
1857    fermie1=zero
1858  end if
1859 
1860 !Destroy atom table used for parallelism
1861  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1862 
1863 end subroutine newfermie1