TABLE OF CONTENTS


ABINIT/dfpt_nstpaw [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nstpaw

FUNCTION

 Initially designed for PAW approach, but works also for NCPP.
 This routine compute the non-stationary expression for the
 second derivative of the total energy, for a whole row of
 mixed derivatives (including diagonal terms contributing
 to non-stationnary 2nd-order total energy).
 Compared with NC-pseudopotentials, PAW contributions include:
  - changes of the overlap between 0-order wave-functions,
  - on-site contributions.

COPYRIGHT

 Copyright (C) 2010-2018 ABINIT group (MT, AM)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k
  cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
  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
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors
  cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors
  docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
  doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz 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)
  eigen1(2*mband*mband*nkpt_rbz*nsppol)=1st-order eigenvalues at k,q (hartree)
  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
  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.
  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(nfftf,nkxc)=exchange and correlation kernel
  mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90)
  mkmem =number of k points treated by this node.
  mkqmem =number of k+q points treated by this node (GS data).
  mk1mem =number of k points treated by this node (RF data)
  mpert =maximum number of ipert
  mpi_enreg=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
  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
  ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid
  nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
  nhat1(cplex*nfftf,nspden*usepaw)=1st-order compensation charge density (PAW)
  nkpt_rbz=number of k points in the reduced BZ for this perturbation
  nkxc=second dimension of the kxc array
  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 i perturbation
  n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, cplex*nfftf
  occkq(mband*nkpt_rbz*nsppol)=occupation number for each band at each k+q point of the reduced BZ
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced BZ
  paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS
  paw_an1(natom) <type(paw_an_type)>=1st-order paw arrays given on angular mesh for the perturbation (j1)
  paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
  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
  pawang1 <type(pawang_type)>=pawang datastructure containing only the 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
  pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies 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 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
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
  rhor1(cplex*nfftf,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)
  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
  ucvol=unit cell volume in bohr**3.
  usecprj= 1 if cprj, cprjq arrays are stored in memory
  usepaw= 0 for non paw calculation; =1 for paw calculation
  usexcnhat= -PAW only- flag controling use of compensation density in Vxc
  useylmgr1= 1 if ylmgr1 array is allocated
  vhartr1(cplex*nfft)=1-order Hartree potential
  vpsp1(cplex*nfftf)=first-order derivative of the ionic potential
  vtrial(nfftf,nspden)=GS potential (Hartree).
  vtrial1(cplex*nfftf,nspden)= RF 1st-order potential (Hartree).
  vxc(nfftf,nspden)=XC GS potential
  wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
  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
  ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q

OUTPUT

  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
  d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
  d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
  d2ovl(2,3,mpert,3,mpert*usepaw)=overlap contributions to the 2DTEs (PAW only)
  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]]

NOTES

   We perform here the computation of
     delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
     see PRB 78, 035105 (2008), Eq. (42) [[cite:Audouze2008]]

PARENTS

      dfpt_scfcv

CHILDREN

      appdig,destroy_hamiltonian,destroy_mpi_enreg,destroy_rf_hamiltonian
      dfpt_accrho,dfpt_atm2fft,dfpt_mkcore,dfpt_mkvxc,dfpt_mkvxc_noncoll
      dfpt_mkvxcstr,dfpt_sygra,dfpt_vlocal,dotprod_g,dotprod_vn,fftpac
      getcprj,getdc1,getgh1c,hartrestr,init_hamiltonian,init_rf_hamiltonian
      initmpi_seq,initylmg,kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian
      load_kprime_hamiltonian,load_spin_hamiltonian,mkffnl,mkkin,mkkpg,occeig
      paw_an_reset_flags,paw_ij_free,paw_ij_init,paw_ij_nullify
      paw_ij_reset_flags,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawcprj_get
      pawdfptenergy,pawdij2e1kb,pawdijfr,pawmkrho,pawnhatfr,pawrhoij_alloc
      pawrhoij_free,pawrhoij_init_unpacked,pawrhoij_mpisum_unpacked,projbd
      stresssym,symrhg,timab,vlocalstr,wfk_close,wfk_open_read,wfk_read_bks
      wrtout,xmpi_barrier,xmpi_sum

SOURCE

 186 subroutine dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,&
 187 &                  eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,&
 188 &                  kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,&
 189 &                  mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nhat1,&
 190 &                  nkpt_rbz,nkxc,npwarr,npwar1,nspden,nspinor,nsppol,nsym1,n3xccc,occkq,occ_rbz,&
 191 &                  paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,&
 192 &                  pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,&
 193 &                  ucvol,usecprj,usepaw,usexcnhat,useylmgr1,vhartr1,vpsp1,vtrial,vtrial1,vxc,&
 194 &                  wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr1)
 195 
 196  use defs_basis
 197  use defs_datatypes
 198  use defs_abitypes
 199  use m_xmpi
 200  use m_errors
 201  use m_abicore
 202  use m_wfk
 203  use m_hamiltonian
 204  use m_cgtools
 205  use m_nctk
 206 
 207  use m_time,     only : timab
 208  use m_io_tools, only : file_exists
 209  use m_geometry, only : stresssym
 210  use m_dynmat,   only : dfpt_sygra
 211  use m_mpinfo,   only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle
 212  use m_hdr,      only : hdr_skip
 213  use m_occ,      only : occeig
 214  use m_pawang,   only : pawang_type
 215  use m_pawrad,   only : pawrad_type
 216  use m_pawtab,   only : pawtab_type
 217  use m_paw_an,   only : paw_an_type, paw_an_reset_flags
 218  use m_paw_ij,   only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
 219  use m_pawfgrtab,only : pawfgrtab_type
 220  use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, pawrhoij_nullify, &
 221 &                       pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_get_nspden
 222  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get, pawcprj_copy
 223  use m_pawdij,   only : pawdijfr
 224  use m_pawfgr,   only : pawfgr_type
 225  use m_paw_mkrho,only : pawmkrho
 226  use m_paw_nhat, only : pawnhatfr
 227  use m_paw_dfpt, only : pawdfptenergy
 228  use m_kg,       only : mkkin, kpgstr, mkkpg
 229  use m_fft,      only : fftpac
 230  use m_spacepar, only : hartrestr, symrhg
 231  use m_initylmg, only : initylmg
 232  use m_mkffnl,   only : mkffnl
 233  use m_getgh1c,  only : getgh1c, getdc1
 234  use m_dfpt_mkrho, only : dfpt_accrho
 235  use m_atm2fft,    only : dfpt_atm2fft
 236  use m_mkcore,     only : dfpt_mkcore
 237  use m_dfpt_mkvxc,    only : dfpt_mkvxc, dfpt_mkvxc_noncoll
 238  use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
 239  use m_mklocl,     only : dfpt_vlocal, vlocalstr
 240  use m_cgprj,           only : getcprj
 241 
 242 !This section has been created automatically by the script Abilint (TD).
 243 !Do not modify the following lines by hand.
 244 #undef ABI_FUNC
 245 #define ABI_FUNC 'dfpt_nstpaw'
 246 !End of the abilint section
 247 
 248  implicit none
 249 
 250 !Arguments -------------------------------
 251 !scalars
 252  integer,intent(in) :: cplex,idir,ipert,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpw,mpw1
 253  integer,intent(in) :: ncpgr,nfftf,nkpt_rbz,nkxc,nspden,nspinor,nsppol,nsym1
 254  integer,intent(in) :: n3xccc,usecprj,usepaw,usexcnhat,useylmgr1
 255  real(dp),intent(in) :: gsqcut,ucvol
 256  real(dp),intent(out) :: eovl1
 257  type(datafiles_type),intent(in) :: dtfil
 258  type(dataset_type),intent(in) :: dtset
 259  type(MPI_type),intent(in) :: mpi_enreg
 260  type(pawang_type),intent(in) :: pawang,pawang1
 261  type(pawfgr_type),intent(in) :: pawfgr
 262  type(pseudopotential_type),intent(in) :: psps
 263 !arrays
 264  integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol)
 265  integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom)
 266  integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))
 267  integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
 268  integer,intent(in) :: ngfftf(18),npwarr(nkpt_rbz),npwar1(nkpt_rbz)
 269  integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1)
 270  integer,intent(inout) :: blkflg(3,mpert,3,mpert)
 271  real(dp),intent(in) :: cg(2,mpw*nspinor*dtset%mband*mkmem*nsppol)
 272  real(dp),intent(in) :: cgq(2,mpw1*nspinor*dtset%mband*mkqmem*nsppol)
 273  real(dp),intent(in) :: cg1(2,mpw1*nspinor*dtset%mband*mk1mem*nsppol)
 274  real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*nsppol)
 275  real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*nsppol)
 276  real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*nsppol),eigen0(dtset%mband*nkpt_rbz*nsppol)
 277  real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol)
 278  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kpt_rbz(3,nkpt_rbz)
 279  real(dp),intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden),nhat1(cplex*nfftf,nspden*usepaw)
 280  real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol)
 281  real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*nsppol)
 282  real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))
 283  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
 284  real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*dtset%natom),rhog(2,nfftf)
 285  real(dp),intent(in) :: rhor(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden),rmet(3,3),rprimd(3,3)
 286  real(dp),intent(in) :: vhartr1(cplex*nfftf),vtrial1(cplex*nfftf,nspden),vxc(nfftf,nspden)
 287  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom)
 288  real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 289  real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
 290  real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
 291  real(dp),target,intent(in) :: vpsp1(cplex*nfftf),vtrial(nfftf,nspden),xccc3d1(cplex*n3xccc)
 292  real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert)
 293  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw)
 294  type(pawcprj_type),intent(in) :: cprj(dtset%natom,nspinor*dtset%mband*mkmem*nsppol*usecprj)
 295  type(pawcprj_type),intent(in) :: cprjq(dtset%natom,nspinor*dtset%mband*mkqmem*nsppol*usecprj)
 296  type(paw_an_type),intent(in) :: paw_an(:)
 297  type(paw_an_type),intent(inout) :: paw_an1(:)
 298  type(paw_ij_type),intent(in) :: paw_ij(:)
 299  type(paw_ij_type),intent(inout) :: paw_ij1(:)
 300  type(pawfgrtab_type),intent(inout) :: pawfgrtab(:)
 301  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*usepaw)
 302  type(pawrhoij_type),intent(in) :: pawrhoij(:)
 303  type(pawrhoij_type),intent(in) :: pawrhoij1(:)
 304  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw)
 305 
 306 !Local variables-------------------------------
 307 !scalars
 308  integer,parameter :: tim_fourwf=18,tim_getgh1c=2,tim_projbd=3,formeig1=1
 309  integer :: bd2tot_index,bdtot_index,berryopt,bufsz,choice,cpopt,counter,cplex_rhoij,ddkcase
 310  integer :: dimffnl,dimffnl1,dimffnl1_idir1,dimylmgr1
 311  integer :: ia,iatom,iband,ibg,ibgq,ibg1,icg,icgq,icg1,ider,idir0,idir1,idir_cprj
 312  integer :: ierr,ii,ikg,ikg1,ikpt,ikpt_me,ilmn,iorder_cprj,ipert1
 313  integer :: ispden,isppol,istwf_k,istr,istr1,itypat,jband,jj,kdir1,kpert1,master,mcgq,mcprjq
 314  integer :: mdir1,me,mpert1,my_natom,my_comm_atom,my_nsppol,nband_k,nband_kocc,need_ylmgr1
 315  integer :: nfftot,nkpg,nkpg1,nkpt_me,npw_,npw_k,npw1_k,nspden_rhoij
 316  integer :: nvh1,nvxc1,nzlmopt_ipert,nzlmopt_ipert1,optlocal,optnl
 317  integer :: option,opt_gvnl1,sij_opt,spaceworld,usevnl,use_rhoijim,wfcorr,ik_ddk
 318  real(dp) :: arg,doti,dotr,dot1i,dot1r,dot2i,dot2r,dot3i,dot3r,elfd_fact,invocc,lambda,wtk_k
 319  logical :: force_recompute,has_dcwf,has_dcwf2,has_drho,has_ddk_file
 320  logical :: is_metal,is_metal_or_qne0,need_ddk_file,need_pawij10
 321  logical :: need_wfk,need_wf1,paral_atom,qne0,t_exist
 322  character(len=500) :: msg
 323  character(len=fnlen) :: fiwfddk(3)
 324  type(gs_hamiltonian_type) :: gs_hamkq
 325  type(rf_hamiltonian_type) :: rf_hamkq
 326  type(MPI_type) :: mpi_enreg_seq
 327 !arrays
 328  integer :: ddkfil(3),my_spintab(2),nband_tmp(1),npwar1_tmp(1)
 329  integer,allocatable :: jpert1(:),jdir1(:),kg1_k(:,:),kg_k(:,:)
 330  integer,pointer :: my_atmtab(:)
 331  real(dp) :: dum1(1,1),dum2(1,1),dum3(1,1),epawnst(2),kpoint(3),kpq(3)
 332  real(dp) :: sumelfd(2),symfact(3),tsec(2),ylmgr_dum(1,3,1)
 333  real(dp),allocatable :: buffer(:),ch1c(:,:,:,:),cs1c(:,:,:,:)
 334  real(dp),allocatable :: cwave0(:,:),cwavef(:,:),dcwavef(:,:)
 335  real(dp),allocatable :: doccde_k(:),doccde_kq(:)
 336  real(dp),allocatable :: dnhat1(:,:),drhoaug1(:,:,:,:)
 337  real(dp),allocatable :: drhor1(:,:),drho1wfg(:,:),drho1wfr(:,:,:)
 338  real(dp),allocatable :: d2nl_elfd(:,:),dkinpw(:)
 339  real(dp),allocatable :: d2nl_k(:,:),d2ovl_drho(:,:,:,:,:),d2ovl_k(:,:)
 340  real(dp),allocatable :: eig_k(:),eig_kq(:),eig1_k(:)
 341  real(dp),allocatable,target :: e1kbfr_spin(:,:,:,:,:,:),ffnlk(:,:,:,:),ffnl1(:,:,:,:)
 342  real(dp),allocatable :: gh1(:,:),gs1(:,:),gvnl1(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:)
 343  real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:),rhotmp(:,:),rocceig(:,:)
 344  real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:),vtmp1(:,:),vxc10(:,:)
 345  real(dp),allocatable,target :: work(:,:,:),e1kb_work(:,:,:,:)
 346  real(dp),pointer :: e1kbfr(:,:,:,:,:),e1kb_ptr(:,:,:,:)
 347  real(dp),pointer :: ffnl1_idir1(:,:,:,:),vhartr01(:),vpsp1_idir1(:),xccc3d1_idir1(:)
 348  type(pawcprj_type),allocatable :: dcwaveprj(:,:)
 349  type(pawcprj_type),allocatable,target :: cwaveprj0(:,:)
 350  type(pawcprj_type),pointer :: cwaveprj0_idir1(:,:)
 351  type(paw_ij_type),allocatable :: paw_ij10(:,:)
 352  type(pawrhoij_type),target,allocatable :: pawdrhoij1(:,:)
 353  type(pawrhoij_type),pointer :: pawdrhoij1_unsym(:,:)
 354  type(wfk_t) :: ddks(3)
 355 
 356 ! *********************************************************************
 357 
 358  DBG_ENTER("COLL")
 359 
 360 !Keep track of total time spent in dfpt_nstpaw
 361  call timab(566,1,tsec)
 362 
 363 !Not valid for PrintBandByBand
 364  if (dtset%prtbbb/=0) then
 365    MSG_BUG('not yet valid for prtbbb/=0!')
 366  end if
 367 
 368 !NCPP restrictions
 369  if (usepaw==0) then
 370 !  cprj cannot be used
 371    if (usecprj/=0) then
 372      MSG_BUG('NCPP: usecprj should be 0!')
 373    end if
 374 !  d2ovl cannot be used
 375    if (size(d2ovl)/=0) then
 376      MSG_BUG('NCPP: d2ovl should not be allocated!')
 377    end if
 378  end if
 379 
 380 !PAW restrictions
 381  if (usepaw==1) then
 382 !  Test on FFT grid sizes
 383    if (pawfgr%nfft/=nfftf) then
 384      MSG_BUG('PAW: wrong values for nfft, nfftf!')
 385    end if
 386 !  Test gradients of cprj
 387    if (ipert<=dtset%natom.and.ncpgr/=3) then
 388      MSG_BUG('PAW: wrong value of ncpgr for ipert<=natom!')
 389    end if
 390    if (ipert==dtset%natom+1.and.ncpgr/=1) then
 391      MSG_BUG('PAW: wrong value of ncpgr for ipert=natom+1!')
 392    end if
 393    if (ipert==dtset%natom+2.and.ncpgr/=3) then
 394      MSG_BUG('PAW: wrong value of ncpgr for ipert=natom+2!')
 395    end if
 396    if ((ipert==dtset%natom+3.or.ipert==dtset%natom+4).and.ncpgr/=1) then
 397      MSG_BUG('PAW: wrong value of ncpgr for ipert=natom+3 or 4!')
 398    end if
 399 !  Test on availability of DijHartree and XC on-site potentials
 400    if (mpi_enreg%my_natom>0.and.ipert/=dtset%natom+1)  then
 401      if (paw_ij1(1)%has_dijhartree==0.or.paw_an1(1)%has_vxc==0) then
 402        msg='PAW: paw_ij1%dijhartree and paw=_an1%vxc1 should be allocated !'
 403      end if
 404    end if
 405  end if
 406 
 407 !Set up parallelism
 408  master=0;me=mpi_enreg%me_kpt
 409  spaceworld=mpi_enreg%comm_cell
 410  paral_atom=(mpi_enreg%my_natom/=dtset%natom)
 411  my_comm_atom=mpi_enreg%comm_atom
 412  my_natom=mpi_enreg%my_natom
 413  my_atmtab=>mpi_enreg%my_atmtab
 414  my_spintab=mpi_enreg%my_isppoltab
 415  my_nsppol=count(my_spintab==1)
 416 
 417 !Fake MPI data to be used in sequential calls to parallel routines
 418  call initmpi_seq(mpi_enreg_seq)
 419  mpi_enreg_seq%my_natom=dtset%natom
 420 
 421 !Compute effective number of k-points
 422  nkpt_me=nkpt_rbz
 423  if(xmpi_paral==1)then
 424    nkpt_me=0
 425    do isppol=1,nsppol
 426      do ikpt=1,nkpt_rbz
 427        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
 428        if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) nkpt_me=nkpt_me+1
 429      end do
 430    end do
 431  end if
 432 
 433 !Sizes for WF at k+q
 434  mcgq=mpw1*nspinor*dtset%mband*mkqmem*nsppol
 435  mcprjq=nspinor*dtset%mband*mkqmem*nsppol*usecprj
 436 
 437 !Check ddk files (needed to compute electric field perturbations)
 438  ddkfil(:)=0
 439  do idir1=1,3
 440    ddkcase=idir1+dtset%natom*3
 441    call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(idir1))
 442    t_exist = file_exists(fiwfddk(idir1))
 443 
 444    if (.not. t_exist) then
 445      ! Try netcdf file.
 446      t_exist = file_exists(nctk_ncify(fiwfddk(idir1)))
 447      if (t_exist) then
 448        fiwfddk(idir1) = nctk_ncify(fiwfddk(idir1))
 449        write(msg,"(3a)")"- File: ",trim(fiwfddk(idir1))," does not exist but found netcdf file with similar name."
 450        call wrtout(std_out,msg,'COLL')
 451      end if
 452    end if
 453 
 454    if (t_exist) then
 455      ddkfil(idir1)=20+idir1 ! Note the use of unit numbers 21, 22 and 23
 456    end if
 457  end do
 458  has_ddk_file=(any(ddkfil(:)>0))
 459 
 460 !Define the set of perturbations (j1)=(ipert1,idir1)
 461 !The first perturbation must be (j1)=(j2)=(ipert,idir)
 462 !because we need to compute <g|H^(j2)-Eps.S^(j2)|u0> first.
 463  if (ipert/=dtset%natom+1) then
 464    mpert1=0
 465    ABI_ALLOCATE(jpert1,(mpert))
 466    jpert1 = 0
 467    if (ipert/=dtset%natom+2.or.has_ddk_file) then
 468      mpert1=mpert1+1;jpert1(mpert1)=ipert
 469    end if
 470    do ipert1=1,mpert
 471      if (ipert1/=ipert.and.&
 472 &     (ipert1<=dtset%natom.or.&
 473 &     (ipert1==dtset%natom+2.and.has_ddk_file).or.&
 474 &     ((ipert>dtset%natom.and.ipert/=dtset%natom+5).and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)).or. &
 475 &     ((ipert1==dtset%natom+2).and.has_ddk_file))) then
 476        mpert1=mpert1+1;jpert1(mpert1)=ipert1
 477      end if
 478    end do
 479  else
 480    mpert1=1
 481    ABI_ALLOCATE(jpert1,(mpert1))
 482    jpert1(1)=dtset%natom+1
 483  end if
 484  mdir1=3
 485  ABI_ALLOCATE(jdir1,(mdir1))
 486  jdir1(1:3)= (/ (idir1,idir1=1,3) /)
 487  jdir1(1)=idir;jdir1(idir)=1
 488 
 489 !Index of strain perturbation, if any
 490  istr=idir;if (ipert==dtset%natom+4) istr=idir+3
 491 
 492 !Open ddk WF file(s)
 493  if (has_ddk_file) then
 494    do kdir1=1,mdir1
 495      idir1=jdir1(kdir1)
 496      if (ddkfil(idir1)/=0) then
 497        write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(idir1))
 498        call wrtout(std_out,msg,'COLL')
 499        call wrtout(ab_out,msg,'COLL')
 500        call wfk_open_read(ddks(idir1),fiwfddk(idir1),formeig1,dtset%iomode,ddkfil(idir1),spaceworld)
 501      end if
 502    end do
 503  end if
 504 
 505 !Zero only portion of matrix to be computed here
 506  d2nl(:,:,1:dtset%natom+4,idir,ipert)=zero
 507  if (usepaw==1) d2ovl(:,:,1:dtset%natom+4,idir,ipert)=zero
 508 
 509 !Update list of computed matrix elements
 510  do kpert1=1,mpert1
 511    ipert1=jpert1(kpert1)
 512    do kdir1=1,mdir1
 513      idir1=jdir1(kdir1)
 514      if ((ipert1<=dtset%natom).or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4.or.&
 515 &     (ipert1==dtset%natom+1.and.((ddkfil(idir1)/=0).or.(dtset%rfdir(idir1)/=0.and.idir1<=idir))).or.&
 516 &     (ipert1==dtset%natom+2.and.ddkfil(idir1)/=0).or.&
 517 &     ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.ddkfil(idir1)/=0)) then
 518        blkflg(idir1,ipert1,idir,ipert)=1
 519      end if
 520    end do
 521  end do
 522 
 523 !Initialize most of the (1st-order) Hamiltonian
 524 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
 525 !2) Perform the setup needed for the non-local factors:
 526  call init_hamiltonian(gs_hamkq,psps,pawtab,nspinor,nsppol,nspden,dtset%natom,&
 527 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
 528 & paw_ij=paw_ij,mpi_atmtab=my_atmtab,comm_atom=my_comm_atom,mpi_spintab=mpi_enreg%my_isppoltab,&
 529 & usecprj=usecprj,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
 530 
 531 !Variables common to all perturbations
 532  arg=maxval(occ_rbz)-minval(occ_rbz)
 533  qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14)
 534  is_metal=((dtset%occopt>=3.and.dtset%occopt<=8).or.(abs(arg)>tol8))
 535  is_metal_or_qne0=((is_metal).or.(qne0))
 536  ABI_ALLOCATE(ch1c,(2,dtset%mband,dtset%mband,nkpt_me))
 537  ch1c(:,:,:,:)=zero
 538  nzlmopt_ipert=0;nzlmopt_ipert1=0
 539  if (usepaw==1) then
 540    ABI_ALLOCATE(d2ovl_drho,(2,3,mpert,3,mpert))
 541    d2ovl_drho=zero
 542    if (dtset%pawnzlm/=0) then
 543      nzlmopt_ipert=1;if (dtset%nstep<2) nzlmopt_ipert=-1
 544      nzlmopt_ipert1=-1
 545    end if
 546    if (is_metal_or_qne0) then
 547      ABI_ALLOCATE(cs1c,(2,dtset%mband,dtset%mband,nkpt_me))
 548      cs1c(:,:,:,:)=zero
 549    end if
 550  end if
 551 
 552 !Force the recomputation of on-site potentials and DijHartree
 553  if (usepaw==1) then
 554    call paw_an_reset_flags(paw_an1)
 555    call paw_ij_reset_flags(paw_ij1,dijhartree=.true.)
 556  end if
 557 
 558 !LOOP OVER PERTURBATION TYPES (j1)
 559  do kpert1=1,mpert1
 560    ipert1=jpert1(kpert1)
 561 
 562 !  Flag for use of DDK file
 563    need_ddk_file=(has_ddk_file.and.(ipert1==dtset%natom+1.or.ipert1==dtset%natom+2))
 564 
 565 !  Factor to be applied for electric Field (Eff. charges and piezo. tensor are "minus" d2E)
 566    elfd_fact=one
 567    if ((ipert <=dtset%natom.or.ipert ==dtset%natom+3.or.ipert ==dtset%natom+4).and. &
 568 &   (ipert1==dtset%natom+2)) elfd_fact=-one
 569    if ((ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and. &
 570 &   (ipert ==dtset%natom+2)) elfd_fact=-one
 571 
 572 !  We want to compute delta_u^(j1))=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
 573 !  see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)
 574    has_dcwf=.false.;has_dcwf2=.false.;has_drho=.false.
 575    if (usepaw==1) then
 576      has_dcwf =(ipert1/=dtset%natom+2)
 577      has_dcwf2=(ipert /=dtset%natom+2)
 578      has_drho =(has_dcwf.and.ipert1/=dtset%natom+1)
 579    end if
 580 
 581 !  Select which WF are needed
 582    need_wfk=(.true.)
 583    need_wf1=(.true.)
 584 
 585 !  Initialize data for NL 1st-order (j1) hamiltonian
 586    call init_rf_hamiltonian(cplex,gs_hamkq,ipert1,rf_hamkq,mpi_spintab=[0,0])
 587 
 588 !  The following contributions are needed only for non-DDK perturbation:
 589 !  - Frozen part of 1st-order Dij
 590 !  - Contribution from local potential to dynamical matrix (due to Vxc^(j1)(tild_nc)+VH^(j1)(tild_nZc))
 591    if (ipert/=dtset%natom+1.and.ipert1/=dtset%natom+1) then
 592 
 593 !    Allocations
 594      force_recompute=(usepaw==0) ! This is dangerous...
 595      nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
 596      nvxc1=0;if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvxc1=nspden
 597      nvh1=0;if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) nvh1=1
 598      ABI_ALLOCATE(vxc10,(cplex*nfftf,nvxc1))
 599      ABI_ALLOCATE(vhartr01,(cplex*nfftf*nvh1))
 600      need_pawij10=(usepaw==1)
 601      if (need_pawij10) then
 602        ABI_DATATYPE_ALLOCATE(paw_ij10,(my_natom,mdir1))
 603        ABI_ALLOCATE(e1kbfr_spin,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,cplex,mdir1,my_nsppol))
 604      else
 605        ABI_DATATYPE_ALLOCATE(paw_ij10,(0,0))
 606      end if
 607 
 608 !    LOOP OVER PERTURBATION DIRECTIONS
 609      do kdir1=1,mdir1
 610        idir1=jdir1(kdir1)
 611        istr1=idir1;if (ipert1==dtset%natom+4) istr1=idir1+3
 612 
 613 !      Get first-order local potential and first-order pseudo core density
 614        if (ipert==ipert1.and.idir==idir1.and.(.not.force_recompute)) then
 615          vpsp1_idir1 => vpsp1
 616          xccc3d1_idir1 => xccc3d1
 617        else
 618          ABI_ALLOCATE(vpsp1_idir1,(cplex*nfftf))
 619          ABI_ALLOCATE(xccc3d1_idir1,(cplex*n3xccc))
 620          if (usepaw==1) then
 621            call dfpt_atm2fft(gs_hamkq%atindx,cplex,gmet,gprimd,gsqcut,istr1,ipert1,&
 622 &           mgfftf,psps%mqgrid_vl,dtset%natom,1,nfftf,ngfftf,dtset%ntypat,&
 623 &           ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
 624 &           atmrhor1=xccc3d1_idir1,atmvlocr1=vpsp1_idir1,optv_in=1,optn_in=n3xccc/nfftf,optn2_in=1,&
 625 &           vspl=psps%vlspl)
 626          else
 627            if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
 628              call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
 629 &             nattyp,nfftf,ngfftf,dtset%ntypat,dtset%paral_kgb,ph1df,psps%qgrid_vl,ucvol,&
 630 &             psps%vlspl,vpsp1_idir1)
 631            else
 632              call dfpt_vlocal(gs_hamkq%atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_vl,&
 633 &             dtset%natom,nattyp,nfftf,ngfftf,dtset%ntypat,ngfftf(1),ngfftf(2),ngfftf(3),dtset%paral_kgb,&
 634 &             ph1df,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_idir1,xred)
 635            end if
 636            if(psps%n1xccc/=0)then
 637              call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,ngfftf(1),psps%n1xccc,&
 638 &             ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,&
 639 &             psps%xcccrc,psps%xccc1d,xccc3d1_idir1,xred)
 640            end if
 641          end if
 642        end if
 643 
 644 !      Compute 1st-order non-local factors (Dij^(j1)_fr)
 645        if (need_pawij10) then
 646          call paw_ij_nullify(paw_ij10(:,idir1))
 647          call paw_ij_init(paw_ij10(:,idir1),cplex,nspinor,dtset%nsppol,dtset%nspden,&
 648 &         0,dtset%natom,dtset%ntypat,dtset%typat,pawtab,has_dijfr=1,&
 649 &         mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
 650          if (ipert/=ipert1.or.idir/=idir1.or.force_recompute) then
 651            option=0
 652            call pawdijfr(cplex,gprimd,idir1,ipert1,my_natom,dtset%natom,nfftf,ngfftf,&
 653 &           nspden,nsppol,dtset%ntypat,option,paw_ij10(:,idir1),pawang,pawfgrtab,pawrad,&
 654 &           pawtab,dtset%qptn,rprimd,ucvol,vpsp1_idir1,vtrial,vxc,xred,&
 655 &           mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
 656          else
 657            do iatom=1,my_natom
 658              paw_ij10(iatom,idir1)%has_dijfr=paw_ij1(iatom)%has_dijfr
 659              if (paw_ij1(iatom)%has_dijfr==2) paw_ij10(iatom,idir1)%dijfr=paw_ij1(iatom)%dijfr
 660            end do
 661          end if
 662        end if
 663 
 664 !      Get first-order exchange-correlation potential (core-correction contribution only)
 665        if (nvxc1>0) then
 666          if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
 667            option=0
 668            call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,&
 669 &           nhat,nhat1,nkxc,nspden,n3xccc,option,dtset%paral_kgb,dtset%qptn,rhor,rhor1,rprimd,&
 670 &           usepaw,usexcnhat,vxc10,xccc3d1_idir1)
 671          else
 672 !          Non-collinear magnetism (should the second nkxc be nkxc_cur ?)
 673            if (nspden==4) then
 674              option=0
 675              call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum1,0,dum2,0,dum3,0,nkxc,&
 676 &             nspden,n3xccc,1,option,dtset%paral_kgb,dtset%qptn,dum1,dum1,rprimd,0,vxc,&
 677 &             vxc10,xccc3d1_idir1)
 678            else
 679              call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfftf,ngfftf,dum2,0,dum3,0,nkxc,&
 680 &             nspden,n3xccc,0,dtset%paral_kgb,dtset%qptn,dum1,rprimd,0,vxc10,xccc3d1_idir1)
 681            end if
 682          end if
 683        end if
 684 
 685 !      Get first-order Hartree potential (metric tensor contribution only)
 686        if (nvh1>0) then
 687          call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,dtset%natom,&
 688 &         nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr01)
 689        end if
 690 
 691 !      Get Hartree + xc + local contributions to dynamical matrix or elastic tensor
 692        if (ipert1<=dtset%natom.or.ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
 693          if (usepaw==0) then
 694 !          vxc1 is integrated with the total 1st-order density (rhor1 )
 695 !          vpsp1 is integrated with the 1st-order pseudo density (rhor1)
 696            call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol)
 697            call dotprod_vn(cplex,rhor1,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol)
 698          else if (usexcnhat/=0) then
 699 !          vxc1 is integrated with the total 1st-order density (rhor1 including nhat1)
 700 !          vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1)
 701            ABI_ALLOCATE(rhotmp,(cplex*nfftf,1))
 702            rhotmp(:,1)=rhor1(:,1)-nhat1(:,1)
 703            call dotprod_vn(cplex,rhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol)
 704            call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol)
 705            ABI_DEALLOCATE(rhotmp)
 706          else
 707 !          vxc1 is integrated with the 1st-order pseudo density (rhor1 without nhat1)
 708 !          vpsp1 is integrated with the 1st-order pseudo density (rhor1 without nhat1)
 709            ABI_ALLOCATE(rhotmp,(cplex*nfftf,nspden))
 710            rhotmp(:,:)=rhor1(:,:)-nhat1(:,:)
 711            call dotprod_vn(cplex,rhotmp,dot1r,dot1i,nfftf,nfftot,nspden,2,vxc10,ucvol)
 712            call dotprod_vn(cplex,rhotmp,dot2r,dot2i,nfftf,nfftot,1,2,vpsp1_idir1,ucvol)
 713            ABI_DEALLOCATE(rhotmp)
 714          end if
 715          if (nvh1>0) then
 716            call dotprod_vn(cplex,rhor1,dot3r,dot3i,nfftf,nfftot,1,2,vhartr01,ucvol)
 717          else
 718            dot3r=zero ; dot3i=zero
 719          end if
 720 !        Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2
 721          dotr=dot1r+dot2r+dot3r;doti=dot1i+dot2i+dot3i
 722 !        In case ipert = natom+2, these lines compute the local part
 723 !        of the Born effective charges from phonon and electric
 724 !        field type perturbations, see eq. 43 of X. Gonze and C. Lee, PRB 55, 10355 (1997) [[cite:Gonze1997a]]
 725 !        The minus sign is due to the fact that the effective charges
 726 !        are minus the second derivatives of the energy
 727          if (ipert/=dtset%natom+1) then
 728            d2lo(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr
 729            d2lo(2,idir1,ipert1,idir,ipert)=elfd_fact*doti
 730          end if
 731        end if ! ipert1<=natom
 732 
 733        if (ipert/=ipert1.or.idir/=idir1.or.force_recompute)  then
 734          ABI_DEALLOCATE(vpsp1_idir1)
 735          ABI_DEALLOCATE(xccc3d1_idir1)
 736        end if
 737 
 738 !      End loop on directions
 739      end do
 740 
 741 !    Free memory
 742      ABI_DEALLOCATE(vxc10)
 743      ABI_DEALLOCATE(vhartr01)
 744 
 745    else ! ddk perturbation
 746      d2lo(1:2,1:mdir1,ipert1,idir,ipert)=zero
 747      need_pawij10=.false.
 748      ABI_DATATYPE_ALLOCATE(paw_ij10,(0,0))
 749    end if
 750 
 751 !  Prepare RF PAW files for reading and writing if mkmem, mkqmem or mk1mem==0
 752    iorder_cprj=0
 753 
 754 !  Allocate arrays used to accumulate density change due to overlap
 755    if (has_drho) then
 756      ABI_ALLOCATE(drhoaug1,(cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),mdir1))
 757      ABI_ALLOCATE(drho1wfr,(cplex*dtset%nfft,dtset%nspden,mdir1))
 758      ABI_DATATYPE_ALLOCATE(pawdrhoij1,(my_natom,mdir1))
 759      if (paral_atom) then
 760        ABI_DATATYPE_ALLOCATE(pawdrhoij1_unsym,(dtset%natom,mdir1))
 761      else
 762        pawdrhoij1_unsym => pawdrhoij1
 763      end if
 764      do kdir1=1,mdir1
 765        idir1=jdir1(kdir1)
 766        drho1wfr(:,:,idir1)=zero
 767        cplex_rhoij=max(cplex,dtset%pawcpxocc)
 768        nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
 769        use_rhoijim = 1
 770 !      use_rhoijim = 0 ; if (sum(dtset%qptn(1:3)**2)>1.d-14) use_rhoijim = 1
 771        call pawrhoij_alloc(pawdrhoij1(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,&
 772 &       dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=1,use_rhoij_=0,use_rhoijim=use_rhoijim,&
 773 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=my_atmtab)
 774        if (paral_atom) then
 775          call pawrhoij_alloc(pawdrhoij1_unsym(:,idir1),cplex_rhoij,nspden_rhoij,dtset%nspinor,&
 776 &         dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1,use_rhoijim=use_rhoijim)
 777        else
 778          call pawrhoij_init_unpacked(pawdrhoij1_unsym(:,idir1))
 779        end if
 780      end do
 781    end if
 782 
 783 !  Initialize shifts for global arrays
 784    bdtot_index=0
 785    bd2tot_index=0
 786    ibg=0;icg=0
 787    ibg1=0;icg1=0
 788    ibgq=0;icgq=0
 789 
 790 !  Has to get 1st-order non-local factors before the loop over spins
 791 !  because this needs a communication over comm_atom (=comm_spinkpt)
 792    if (need_pawij10) then
 793      if (my_nsppol<nsppol) then
 794        ABI_ALLOCATE(e1kb_work,(rf_hamkq%dime1kb1,rf_hamkq%dime1kb2,nspinor**2,cplex))
 795      end if
 796      ii=0
 797      do isppol=1,nsppol
 798        if (my_spintab(isppol)==1) ii=ii+1
 799        if (my_spintab(isppol)/=1) e1kb_ptr => e1kb_work
 800        do kdir1=1,mdir1
 801          idir1=jdir1(kdir1)
 802          if (my_spintab(isppol)==1) e1kb_ptr => e1kbfr_spin(:,:,:,:,idir1,ii)
 803          call pawdij2e1kb(paw_ij10(:,idir1),isppol,my_comm_atom,e1kbfr=e1kb_ptr,mpi_atmtab=my_atmtab)
 804        end do
 805      end do
 806      if (my_nsppol<nsppol) then
 807        ABI_DEALLOCATE(e1kb_work)
 808      end if
 809    end if
 810 
 811 !  LOOP OVER SPINS
 812    do isppol=1,nsppol
 813 
 814 !    either this has to be inside the sppol loop or the size of ch1c and cs1c has to be nkpt_me*nsppol
 815      ikpt_me=0
 816 
 817 !    Rewind (k+G) data if needed
 818      ikg=0;ikg1=0
 819 
 820 !    Continue to initialize the GS/RF Hamiltonian
 821      call load_spin_hamiltonian(gs_hamkq,isppol,with_nonlocal=.true.)
 822      if (need_pawij10) then
 823        ii=min(isppol,size(e1kbfr_spin,6))
 824        if (ii>0) e1kbfr => e1kbfr_spin(:,:,:,:,:,ii)
 825      end if
 826 
 827 !    Initialize accumulation of density
 828      if (has_drho) drhoaug1(:,:,:,:)=zero
 829 
 830 !    LOOP OVER K-POINTS
 831      do ikpt=1,nkpt_rbz
 832 
 833 !      Load dimensions for this k-point
 834        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
 835        istwf_k=istwfk_rbz(ikpt)
 836        npw_k=npwarr(ikpt)
 837        npw1_k=npwar1(ikpt)
 838 
 839 !      Skip loop if this k-point is not to be treated by this proc
 840        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
 841          bdtot_index=bdtot_index+nband_k
 842          bd2tot_index=bd2tot_index+2*nband_k**2
 843          cycle ! Skip the rest of the k-point loop
 844        end if
 845 
 846 !      Allocate/initialize local arrays and scalars for this k-point
 847        ABI_ALLOCATE(d2nl_k,(2,3))
 848        ABI_ALLOCATE(d2ovl_k,(2,3))
 849        ABI_ALLOCATE(eig_k,(nband_k))
 850        ABI_ALLOCATE(eig_kq,(nband_k))
 851        ABI_ALLOCATE(eig1_k,(2*nband_k**2))
 852        ABI_ALLOCATE(occ_k,(nband_k))
 853        d2nl_k(:,:)=zero;d2ovl_k(:,:)=zero
 854        eig_k (:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
 855        eig_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
 856        eig1_k(:)=eigen1(1+bd2tot_index:2*nband_k**2+bd2tot_index)
 857        occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
 858        nband_kocc=count(abs(occ_k(:))>tol8)
 859        kpoint(:)=kpt_rbz(:,ikpt)
 860        kpq(:)=kpoint(:);if (ipert1<dtset%natom+3) kpq(:)=kpq(:)+dtset%qptn(1:3)
 861        wtk_k=wtk_rbz(ikpt)
 862        need_ylmgr1=0;dimylmgr1=0
 863        nkpg=0;nkpg1=0
 864        ikpt_me=ikpt_me+1
 865        if (is_metal) then
 866 !        For each pair of active bands (m,n), generates the ratios
 867 !        rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
 868          ABI_ALLOCATE(doccde_k,(nband_k))
 869          ABI_ALLOCATE(doccde_kq,(nband_k))
 870          ABI_ALLOCATE(occ_kq,(nband_k))
 871          ABI_ALLOCATE(rocceig,(nband_k,nband_k))
 872          doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index)
 873          doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index)
 874          occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index)
 875          call occeig(doccde_k,doccde_kq,eig_k,eig_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig)
 876        end if
 877 
 878 !      Take care of the npw and kg records in WF and DDK files
 879        if (need_ddk_file) then
 880          do kdir1=1,mdir1
 881            idir1=jdir1(kdir1)
 882            if (ddkfil(idir1)/=0)then
 883              !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt)
 884              ik_ddk = indkpt1(ikpt)
 885              npw_ = ddks(idir1)%hdr%npwarr(ik_ddk)
 886              if (npw_/=npw_k) then
 887                write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')&
 888 &               'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,&
 889 &               'the number of plane waves in the ddk file is equal to', npw_,ch10,&
 890 &               'while it should be ',npw_k
 891                MSG_ERROR(msg)
 892              end if
 893 
 894            end if
 895          end do
 896        end if
 897 
 898 !      Allocate arrays used for NL form factors
 899        ABI_ALLOCATE(kg_k,(3,npw_k))
 900        ABI_ALLOCATE(kg1_k,(3,npw1_k))
 901        ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
 902        ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
 903        if (psps%useylm==1.and.(need_ddk_file.or.ipert1==dtset%natom+1)) need_ylmgr1=1
 904        if (psps%useylm==1.and.(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)) need_ylmgr1=1
 905        dimylmgr1=max(useylmgr1,need_ylmgr1)
 906        ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*dimylmgr1))
 907 
 908 !      Get plane-wave vectors and related data at k
 909        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
 910        if (psps%useylm==1) then
 911          do jj=1,psps%mpsang*psps%mpsang
 912            ylm_k(1:npw_k,jj)=ylm(1+ikg:npw_k+ikg,jj)
 913          end do
 914        end if
 915 
 916 !      Get plane-wave vectors and related data at k+q
 917        kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
 918        if (psps%useylm==1) then
 919          do jj=1,psps%mpsang*psps%mpsang
 920            ylm1_k(1:npw1_k,jj)=ylm1(1+ikg1:npw1_k+ikg1,jj)
 921          end do
 922          if (need_ylmgr1==1.and.useylmgr1/=0) then
 923            do jj=1,psps%mpsang*psps%mpsang
 924              do ia=1,3
 925                ylmgr1_k(1:npw1_k,ia,jj)=ylmgr1(1+ikg1:npw1_k+ikg1,ia,jj)
 926              end do
 927            end do
 928          end if
 929        end if
 930 
 931 !      If Ylm gradients at k+q are needed and not in memory, compute them
 932        if (need_ylmgr1==1.and.useylmgr1==0) then
 933          option=-1;npwar1_tmp(1)=npw1_k;nband_tmp(1)=nband_k
 934          !Subtlety: initylmg is called in sequential mode
 935          call initylmg(gprimd,kg1_k,kpq,1,mpi_enreg_seq,psps%mpsang,&
 936 &         npw1_k,nband_tmp,1,npwar1_tmp,nsppol,option,rprimd,&
 937 &         ylm1_k,ylmgr1_k)
 938        end if
 939 
 940 !      Compute (k+G) vectors
 941        nkpg=0;if(ipert1<=dtset%natom) nkpg=3*dtset%nloalg(3)
 942        ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 943        if (nkpg>0) then
 944          call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
 945        end if
 946 
 947 !      Compute (k+q+G) vectors
 948        nkpg1=0;if(ipert1<=dtset%natom.or.need_ylmgr1==1) nkpg1=3*dtset%nloalg(3)
 949        ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1))
 950        if (nkpg1>0) then
 951          call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k)
 952        end if
 953 
 954 !      Allocate kinetic contributions
 955        ABI_ALLOCATE(dkinpw,(npw_k))
 956        ABI_ALLOCATE(kinpw1,(npw1_k))
 957        dkinpw=zero
 958 !      Compute (1/2) (2 Pi)**2 (k+q+G)**2:
 959 !       call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k)
 960        call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
 961 
 962 !      Compute nonlocal form factors ffnl at (k+G), for all atoms
 963        ider=0;idir0=0
 964        dimffnl=0;if (ipert1<=dtset%natom) dimffnl=1
 965        ABI_ALLOCATE(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
 966        if (ipert1<=dtset%natom) then
 967          call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,&
 968 &         ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,&
 969 &         psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,&
 970 &         psps%useylm,ylm_k,ylmgr_dum)
 971        end if
 972 
 973 !      Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
 974        ider=1;if (ipert1<=dtset%natom) ider=0
 975        if(ipert1==dtset%natom+3.or.ipert1==dtset%natom+4)then
 976          dimffnl1=1;if (ider>=1) dimffnl1=2+5*psps%useylm
 977          idir0=0;if (ider>0.and.psps%useylm==1) idir0=-7
 978        else
 979          dimffnl1=1;if (ider>=1) dimffnl1=2+2*psps%useylm
 980          idir0=0;if (ider>0.and.psps%useylm==1) idir0=4
 981        end if
 982        ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
 983        call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,&
 984 &       ider,idir0,psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,&
 985 &       psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,usepaw,&
 986 &       psps%useylm,ylm1_k,ylmgr1_k)
 987 
 988 !      Extract non-local form factors for H^(j1)
 989        if (ipert1<=dtset%natom) then
 990          ffnl1_idir1 => ffnl1(:,:,:,:)
 991          dimffnl1_idir1=dimffnl1
 992        else
 993          dimffnl1_idir1=1+ider
 994          ABI_ALLOCATE(ffnl1_idir1,(npw1_k,dimffnl1_idir1,psps%lmnmax,psps%ntypat))
 995          ii=1;if (psps%useylm==0) ii=1+ider
 996          do itypat=1,psps%ntypat
 997            do ilmn=1,psps%lmnmax
 998              ffnl1_idir1(1:npw1_k,1:ii,ilmn,itypat)=ffnl1(1:npw1_k,1:ii,ilmn,itypat)
 999            end do
1000          end do
1001        end if
1002 
1003 !      Load k-dependent part in the Hamiltonian datastructure
1004        ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamkq%matblk))
1005        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,&
1006 &       ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
1007        if (size(ffnlk)>0) then
1008          call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnlk)
1009        else
1010          call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1)
1011        end if
1012 
1013 !      Load k+q-dependent part in the Hamiltonian datastructure
1014 !          Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
1015        call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
1016 &       kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1_idir1,&
1017 &       compute_gbound=.true.)
1018        if (qne0) then
1019          ABI_ALLOCATE(ph3d1,(2,npw1_k,gs_hamkq%matblk))
1020          call load_kprime_hamiltonian(gs_hamkq,ph3d_kp=ph3d1,compute_ph3d=.true.)
1021        end if
1022 
1023 !      Load k-dependent part in the 1st-order Hamiltonian datastructure
1024        call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw)
1025 
1026 !      Allocate memory space for one band
1027        if (need_wfk)  then
1028          ABI_ALLOCATE(cwave0,(2,npw_k*nspinor))
1029        end if
1030        if (need_wf1)  then
1031          ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor))
1032        end if
1033        ABI_ALLOCATE(gh1,(2,npw1_k*nspinor))
1034        nullify(cwaveprj0_idir1)
1035        if (usecprj==1) then
1036          ABI_DATATYPE_ALLOCATE(cwaveprj0,(dtset%natom,nspinor))
1037          call pawcprj_alloc(cwaveprj0,ncpgr,gs_hamkq%dimcprj)
1038        end if
1039        if (has_dcwf) then
1040          ABI_ALLOCATE(gs1,(2,npw1_k*nspinor))
1041        else
1042          ABI_ALLOCATE(gs1,(0,0))
1043        end if
1044 
1045 !      LOOP OVER BANDS
1046        do iband=1,nband_k
1047 
1048 !        Skip band if not to be treated by this proc
1049          if (xmpi_paral==1) then
1050            if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle
1051          end if
1052 
1053 !        Extract GS wavefunctions
1054          if (need_wfk) then
1055            cwave0(:,:)=cg(:,1+(iband-1)*npw_k*nspinor+icg:iband*npw_k*nspinor+icg)
1056            if (usecprj==1) then
1057              call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,dtset%natom,iband,ibg,ikpt,iorder_cprj,&
1058 &             isppol,dtset%mband,mkmem,dtset%natom,1,nband_k,nspinor,nsppol,dtfil%unpaw,&
1059 &             mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1060            end if
1061          end if
1062 
1063 !        Extract 1st-order wavefunctions
1064          if (need_wf1) then
1065            cwavef(:,:)=cg1(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)
1066          end if
1067 
1068 !        LOOP OVER PERTURBATION DIRECTIONS
1069          do kdir1=1,mdir1
1070            idir1=jdir1(kdir1)
1071            istr1=idir1;if(ipert1==dtset%natom+4) istr1=idir1+3
1072 
1073 !          Not able to compute if ipert1=(Elect. field) and no ddk WF file
1074            if (ipert1==dtset%natom+2.and.ddkfil(idir1)==0) cycle
1075 
1076 !          Extract 1st-order NL form factors derivatives for this idir1
1077            if (dimffnl1_idir1>=2.and.psps%useylm==1.and.ipert1>dtset%natom) then
1078              do itypat=1,psps%ntypat
1079                do ilmn=1,psps%lmnmax
1080                  ffnl1_idir1(1:npw1_k,2,ilmn,itypat)=ffnl1(1:npw1_k,1+istr1,ilmn,itypat)
1081                end do
1082              end do
1083            end if
1084 
1085 !          Extract ground state projected WF and derivatives in idir1 direction
1086            if (usecprj==1) then
1087              cpopt=1
1088 
1089 !            === Atomic displ. perturbation
1090              if (ipert<=dtset%natom) then
1091                ABI_DATATYPE_ALLOCATE(cwaveprj0_idir1,(dtset%natom,nspinor))
1092                call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj)
1093                if (ipert1<=dtset%natom) then
1094                  call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1)
1095                else
1096                  if (ipert1==dtset%natom+2) then
1097                    idir_cprj=idir1;choice=5
1098                  end if
1099                  if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1100                    idir_cprj=istr1;choice=3
1101                  end if
1102                  call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1)
1103                  call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,&
1104 &                 gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,&
1105 &                 gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,&
1106 &                 gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,&
1107 &                 gs_hamkq%nloalg,gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,&
1108 &                 gs_hamkq%ph1d,gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm)
1109                end if
1110 
1111 !            === Wave-vector perturbation
1112              else if (ipert==dtset%natom+1) then
1113                cwaveprj0_idir1 => cwaveprj0
1114 
1115 !            == Electric field perturbation
1116              else if (ipert==dtset%natom+2) then
1117                ABI_DATATYPE_ALLOCATE(cwaveprj0_idir1,(dtset%natom,nspinor))
1118                call pawcprj_alloc(cwaveprj0_idir1,1,gs_hamkq%dimcprj)
1119                if (ipert1==dtset%natom+2) then
1120                  call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=idir1)
1121                else
1122                  if (ipert1<=dtset%natom) then
1123                    idir_cprj=idir1;choice=2
1124                  end if
1125                  if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1126                    idir_cprj=istr1;choice=3
1127                  end if
1128                  call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1)
1129                  call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,&
1130 &                 gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,&
1131 &                 gs_hamkq%kg_kp,gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,&
1132 &                 gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,&
1133 &                 gs_hamkq%npw_kp,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,&
1134 &                 gs_hamkq%ph3d_kp,gs_hamkq%ucvol,gs_hamkq%useylm)
1135                end if
1136 
1137 !            === Strain perturbation
1138              else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
1139                if ((ipert1==dtset%natom+3.or.ipert1==dtset%natom+4).and.(istr==istr1)) then
1140                  cwaveprj0_idir1 => cwaveprj0
1141                else
1142                  ABI_DATATYPE_ALLOCATE(cwaveprj0_idir1,(dtset%natom,nspinor))
1143                  call pawcprj_alloc(cwaveprj0_idir1,ncpgr,gs_hamkq%dimcprj)
1144                  if (ipert1<=dtset%natom) then
1145                    idir_cprj=idir1;choice=2
1146                  end if
1147                  if (ipert1==dtset%natom+2) then
1148                    idir_cprj=idir1;choice=5
1149                  end if
1150                  if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1151                    idir_cprj=istr1;choice=3
1152                  end if
1153                  call pawcprj_copy(cwaveprj0,cwaveprj0_idir1,icpgr=-1)
1154                  call getcprj(choice,cpopt,cwave0,cwaveprj0_idir1,&
1155 &                 gs_hamkq%ffnl_kp,idir_cprj,gs_hamkq%indlmn,gs_hamkq%istwf_kp,gs_hamkq%kg_kp,&
1156 &                 gs_hamkq%kpg_kp,gs_hamkq%kpt_kp,gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,&
1157 &                 gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,gs_hamkq%nloalg,gs_hamkq%npw_kp,&
1158 &                 gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkpxred,gs_hamkq%ph1d,gs_hamkq%ph3d_kp,&
1159 &                 gs_hamkq%ucvol,gs_hamkq%useylm)
1160                end if
1161              end if ! ipert
1162 
1163            else ! usecprj=0: cwaveprj0_idir1 is not used
1164              cwaveprj0_idir1 => cwaveprj0
1165            end if
1166 
1167 !          Eventually compute 1st-order kinetic operator
1168            if (ipert1==dtset%natom+1) then
1169              call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0)
1170            else if (ipert1==dtset%natom+3.or.ipert1==dtset%natom+4) then
1171              call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr1,kg_k,kpoint,npw_k)
1172            end if
1173 
1174 !          Finalize initialization of 1st-order NL hamiltonian
1175            if (need_pawij10) rf_hamkq%e1kbfr => e1kbfr(:,:,:,:,idir1)
1176 
1177 !          Read DDK wave function (if ipert1=electric field)
1178            if (ipert1==dtset%natom+2) then
1179              usevnl=1
1180              ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor*usevnl))
1181              if (need_ddk_file) then
1182                if (ddkfil(idir1)/=0) then
1183                  !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt)
1184                  ik_ddk = indkpt1(ikpt)
1185                  call wfk_read_bks(ddks(idir1), iband, ik_ddk, isppol, xmpio_single, cg_bks=gvnl1)
1186                else
1187                  gvnl1=zero
1188                end if
1189                if (ipert1==dtset%natom+2) then
1190                  do ii=1,npw1_k*nspinor ! Multiply ddk by +i (to be consistent with getgh1c)
1191                    arg=gvnl1(1,ii)
1192                    gvnl1(1,ii)=-gvnl1(2,ii)
1193                    gvnl1(2,ii)=arg
1194                  end do
1195                end if
1196              else
1197                gvnl1=zero
1198              end if
1199            else
1200              usevnl=0
1201              ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor*usevnl))
1202            end if
1203 
1204 !          Get |H^(j2)-Eps_k_i.S^(j2)|u0_k_i> (VHxc-dependent part not taken into account) and S^(j2)|u0>
1205            lambda=eig_k(iband);berryopt=0;optlocal=0
1206            optnl=0;if (ipert1/=dtset%natom+1.or.idir==idir1) optnl=1
1207            opt_gvnl1=0;if (ipert1==dtset%natom+2) opt_gvnl1=2
1208            sij_opt=-1;if (has_dcwf) sij_opt=1
1209            if (usepaw==0) sij_opt=0
1210            call getgh1c(berryopt,cwave0,cwaveprj0_idir1,gh1,dum1,gs1,gs_hamkq,gvnl1,idir1,ipert1,&
1211 &           lambda,mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
1212            if (sij_opt==1.and.optnl==1) gh1=gh1-lambda*gs1
1213            ABI_DEALLOCATE(gvnl1)
1214 
1215 !          If needed, compute here <delta_u^(j1)_k_i|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>
1216 !          with delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
1217 !          (see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42))
1218 !          This can be rewritten as:
1219 !          -1/2.<u0_k_i|S^(j1)| Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j>
1220 !          The sum over j can be computed with a single call to projbd routine
1221 !          At first call (when j1=j2), ch1c=<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i> is stored
1222 !          For the next calls, it is re-used.
1223            if (has_dcwf.or.(ipert==ipert1.and.idir==idir1.and.usepaw==1)) then
1224 !            note: gvnl1 used as temporary space
1225              ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor))
1226              if (ipert==ipert1.and.idir==idir1) then
1227                option=0;gvnl1=gh1
1228              else
1229                option=1;gvnl1=zero
1230              end if
1231 !            Compute -Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j>
1232              call projbd(cgq,gvnl1,-1,icgq,0,istwf_k,mcgq,0,nband_k,npw1_k,nspinor,&
1233 &             dum1,ch1c(:,1:nband_k,iband,ikpt_me),option,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1234              if (has_dcwf) then
1235                if (ipert==ipert1.and.idir==idir1) gvnl1=gvnl1-gh1
1236                dotr=zero;doti=zero
1237                if (abs(occ_k(iband))>tol8) then
1238 !                Compute: -<u0_k_i|S^(j1)| Sum_{j}[<u0_k+q_j|H^(j2)-Eps_k_i.S^(j2)|u0_k_i>.|u0_k+q_j>
1239                  call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gs1,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1240 !                Add contribution to DDB
1241 !                Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2
1242 !                (-1) factor already present
1243                  d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband)*elfd_fact*dotr
1244                  d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband)*elfd_fact*doti
1245                end if
1246              end if
1247              ABI_DEALLOCATE(gvnl1)
1248            end if
1249 
1250 !          If needed, compute here <delta_u^(j1)_k_i|H-Eps_k_i.S|u^(j2)_k_i>
1251 !          This is equal to <delta_u^(j1)_k_i|H-Eps_k_i.S|delta_u^(j2)_k_i>  (I)
1252 !                          +<delta_u^(j1)_k_i|H-Eps_k_i.S|u^paral^(j2)_k_i>  (II)
1253 !          (u^paral^(j2)_k_i is the part of u^(j2)_k_i parallel to active space : metals)
1254 !          (I) can be rewritten as:
1255 !          Sum_j{ 1/4.<u0_k_i|S^(j1)|u0_k+q_j>.<u0_k+q_j|S^(j2)|u0_k_i>.(Eps_k+q_j-Eps_k_i) }
1256 !          (II) can be rewritten as:
1257 !          Sum_j{1/2.(occ_kq_j-occ_k_i).Eps1_k,q_ij.<u0_k_i|S^(j1)|u0_k+q_j> }
1258 !          where Eps1_k,q_ij=<u0_k+q_j|H^(j2)-1/2(Eps_k+q_j-Eps_k_i)S^(j2)|u0_k_i>
1259 !          At first call (when j1=j2), cs1c=<u0_k_i|S^(j1)|u0_k+q_j> is stored
1260 !          For the next calls, it is re-used.
1261            if (has_dcwf.and.is_metal_or_qne0) then
1262              ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor))
1263              dotr=zero;doti=zero
1264              invocc=zero ; if (abs(occ_k(iband))>tol8) invocc=two/occ_k(iband)
1265              do jband=1,nband_k
1266 !              Computation of cs1c=<u0_k_i|S^(j1)|u0_k+q_j>
1267                if ((ipert==ipert1.and.idir==idir1).or.(abs(occ_k(iband))>tol8)) then
1268                  gvnl1(:,1:npw1_k*nspinor)=cgq(:,1+npw1_k*nspinor*(jband-1)+icgq:npw1_k*nspinor*jband+icgq)
1269                  call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*nspinor,2,gs1,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1270                  if (ipert==ipert1.and.idir==idir1.and.has_dcwf2) then
1271                    cs1c(1,jband,iband,ikpt_me)=dot1r
1272                    cs1c(2,jband,iband,ikpt_me)=dot1i
1273                  end if
1274                end if
1275                if (abs(occ_k(iband))>tol8) then
1276 !                Computation of term (I)
1277                  if (has_dcwf2) then
1278                    arg=eig_kq(jband)-eig_k(iband)
1279                    dot2r=cs1c(1,jband,iband,ikpt_me)
1280                    dot2i=cs1c(2,jband,iband,ikpt_me)
1281                    dotr=dotr+(dot1r*dot2r+dot1i*dot2i)*arg
1282                    doti=doti+(dot1i*dot2r-dot1r*dot2i)*arg
1283                  end if
1284 !                Computation of term (II)
1285                  if (is_metal) then
1286                    if (abs(rocceig(jband,iband))>tol8) then
1287                      ii=2*jband-1+(iband-1)*2*nband_k
1288                      arg=invocc*rocceig(jband,iband)*(eig_k(iband)-eig_kq(jband))
1289                      dot2r=eig1_k(ii);dot2i=eig1_k(ii+1)
1290                      dotr=dotr+arg*(dot1r*dot2r-dot1i*dot2i)
1291                      doti=doti+arg*(dot1r*dot2i+dot2r*dot1i)
1292                    end if
1293                  end if
1294                end if
1295              end do
1296              dotr=quarter*dotr;doti=quarter*doti
1297 !            Note: factor 2 (from d2E/dj1dj2=2E^(j1j2))
1298              d2ovl_k(1,idir1)=d2ovl_k(1,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*dotr
1299              d2ovl_k(2,idir1)=d2ovl_k(2,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*doti
1300              ABI_DEALLOCATE(gvnl1)
1301            end if
1302 
1303 !          Build the matrix element <u0_k_i|H^(j1)-Eps_k_i.S^(j1)|u^(j2)_k,q_i>
1304 !          and add contribution to DDB
1305            if (ipert1/=dtset%natom+1) then
1306              if (abs(occ_k(iband))>tol8) then
1307                call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gh1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1308 !              Case ipert1=natom+2 (electric field):
1309 !              gh1 contains H^(j1)|u0_k_i> (VHxc constant) which corresponds
1310 !              to i.d/dk in Eq. (38) of Gonze, PRB 55, 10355 (1997) [[cite:Gonze1997a]].
1311 !              * if ipert==natom+2, we apply directly Eq. (38)
1312 !              * if ipert/=natom+2, Born effective charges are minus D2E
1313                d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*dotr
1314                d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband)*two*elfd_fact*doti
1315              end if
1316 
1317 !            Or compute localisation tensor (ddk)
1318 !            See M. Veithen thesis Eq(2.5)
1319 !            MT jan-2010: this is probably not correctly implemented for PAW !!!
1320 !            missing terms due to S^(1) and S^(2)
1321            else
1322 !            note: gh1 used as temporary space (to store idir ddk WF)
1323              if (idir==idir1) then
1324                gh1=cwavef
1325              else
1326                if (need_ddk_file.and.ddkfil(idir1)/=0) then
1327                  !ik_ddk = wfk_findk(ddks(idir1), kpt_rbz(:,ikpt)
1328                  ik_ddk = indkpt1(ikpt)
1329                  call wfk_read_bks(ddks(idir1), iband, ik_ddk, isppol, xmpio_single, cg_bks=gh1)
1330                else
1331                  gh1=zero
1332                end if
1333              end if
1334              if (abs(occ_k(iband))>tol8) then
1335                call dotprod_g(dotr,doti,istwf_k,npw1_k*nspinor,2,gh1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1336                call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*nspinor,2,cwave0,gh1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1337                call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1338                dotr=dotr-(dot1r*dot2r-dot1i*dot2i)
1339                doti=doti-(dot1r*dot2i+dot1i*dot2r)
1340                d2nl_k(1,idir1)=d2nl_k(1,idir1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two)
1341                d2nl_k(2,idir1)=d2nl_k(2,idir1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two)
1342              end if
1343            end if
1344 
1345 !          Accumulate here 1st-order density change due to overlap operator changes (if any)
1346            if (has_drho) then
1347              if (abs(occ_k(iband))>tol8) then
1348 !              Compute here delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]
1349 !              (see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42))
1350                ABI_ALLOCATE(dcwavef,(2,npw1_k*nspinor))
1351                ABI_DATATYPE_ALLOCATE(dcwaveprj,(dtset%natom,nspinor))
1352                call pawcprj_alloc(dcwaveprj,0,gs_hamkq%dimcprj)
1353                call getdc1(cgq,cprjq,dcwavef,dcwaveprj,ibgq,icgq,istwf_k,mcgq,&
1354 &               mcprjq,mpi_enreg,dtset%natom,nband_k,npw1_k,nspinor,1,gs1)
1355 
1356 !              Accumulate 1st-order density due to delta_u^(j1)
1357                counter=500*iband;option=1;wfcorr=0
1358                call dfpt_accrho(counter,cplex,cwave0,dcwavef,dcwavef,cwaveprj0_idir1,dcwaveprj,&
1359 &               lambda,dtfil%filstat,gs_hamkq,iband,idir1,ipert1,isppol,dtset%kptopt,&
1360 &               mpi_enreg,dtset%natom,nband_k,1,npw_k,npw1_k,nspinor,occ_k,option,&
1361 &               pawdrhoij1_unsym(:,idir1),dtset%prtvol,drhoaug1(:,:,:,idir1),tim_fourwf,wfcorr,wtk_k)
1362                call pawcprj_free(dcwaveprj)
1363                ABI_DEALLOCATE(dcwavef)
1364                ABI_DATATYPE_DEALLOCATE(dcwaveprj)
1365              end if
1366            end if
1367 
1368            if((usecprj==1).and..not.(associated(cwaveprj0_idir1,cwaveprj0)))then
1369              call pawcprj_free(cwaveprj0_idir1)
1370              ABI_DATATYPE_DEALLOCATE(cwaveprj0_idir1)
1371            end if
1372 
1373 !          End of loops
1374          end do   ! idir1
1375        end do     ! iband
1376 
1377 !      Accumulate contribution of this k-point
1378        d2nl (:,:,ipert1,idir,ipert)=d2nl (:,:,ipert1,idir,ipert)+d2nl_k (:,:)
1379        if (usepaw==1) d2ovl(:,:,ipert1,idir,ipert)=d2ovl(:,:,ipert1,idir,ipert)+d2ovl_k(:,:)
1380 
1381 !      Deallocations of arrays used for this k-point
1382        ABI_DEALLOCATE(gh1)
1383        ABI_DEALLOCATE(gs1)
1384        if (need_wfk)  then
1385          ABI_DEALLOCATE(cwave0)
1386        end if
1387        if (need_wf1)  then
1388          ABI_DEALLOCATE(cwavef)
1389        end if
1390        ABI_DEALLOCATE(kg_k)
1391        ABI_DEALLOCATE(kg1_k)
1392        ABI_DEALLOCATE(ylm_k)
1393        ABI_DEALLOCATE(ylm1_k)
1394        ABI_DEALLOCATE(ylmgr1_k)
1395        ABI_DEALLOCATE(kpg_k)
1396        ABI_DEALLOCATE(kpg1_k)
1397        ABI_DEALLOCATE(d2nl_k)
1398        ABI_DEALLOCATE(d2ovl_k)
1399        ABI_DEALLOCATE(eig_k)
1400        ABI_DEALLOCATE(eig_kq)
1401        ABI_DEALLOCATE(eig1_k)
1402        ABI_DEALLOCATE(occ_k)
1403        if (is_metal)  then
1404          ABI_DEALLOCATE(doccde_k)
1405          ABI_DEALLOCATE(doccde_kq)
1406          ABI_DEALLOCATE(occ_kq)
1407          ABI_DEALLOCATE(rocceig)
1408        end if
1409        ABI_DEALLOCATE(dkinpw)
1410        ABI_DEALLOCATE(kinpw1)
1411        ABI_DEALLOCATE(ph3d)
1412        if (allocated(ph3d1)) then
1413          ABI_DEALLOCATE(ph3d1)
1414        end if
1415        ABI_DEALLOCATE(ffnlk)
1416        ABI_DEALLOCATE(ffnl1)
1417        if (ipert1>dtset%natom)  then
1418          ABI_DEALLOCATE(ffnl1_idir1)
1419        end if
1420        nullify(ffnl1_idir1)
1421        if (usecprj==1) then
1422          call pawcprj_free(cwaveprj0)
1423          ABI_DATATYPE_DEALLOCATE(cwaveprj0)
1424        end if
1425        nullify(cwaveprj0_idir1)
1426 !      Shift arrays
1427        bdtot_index=bdtot_index+nband_k
1428        bd2tot_index=bd2tot_index+2*nband_k**2
1429        if (mkmem/=0) then
1430          ibg=ibg+nspinor*nband_k
1431          icg=icg+npw_k*nspinor*nband_k
1432          ikg=ikg+npw_k
1433        end if
1434        if (mkqmem/=0) then
1435          ibgq=ibgq+nspinor*nband_k
1436          icgq=icgq+npw1_k*nspinor*nband_k
1437        end if
1438        if (mk1mem/=0) then
1439          ibg1=ibg1+nspinor*nband_k
1440          icg1=icg1+npw1_k*nspinor*nband_k
1441          ikg1=ikg1+npw1_k
1442        end if
1443 
1444      end do ! End loop over K-POINTS
1445 
1446 !    Transfer 1st-order density change due to overlap; also take into account the spin.
1447      if(has_drho) then
1448        do kdir1=1,mdir1
1449          idir1=jdir1(kdir1)
1450          call fftpac(isppol,mpi_enreg,nspden,cplex*dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),&
1451 &         cplex*dtset%ngfft(4),dtset%ngfft(5),dtset%ngfft(6),&
1452 &         dtset%ngfft,drho1wfr(:,:,idir1),drhoaug1(:,:,:,idir1),1)
1453        end do
1454      end if
1455 
1456    end do ! End loop over SPINS
1457 
1458 !  Free memory used for this type of perturbation
1459    call destroy_rf_hamiltonian(rf_hamkq)
1460    if (has_drho)  then
1461      ABI_DEALLOCATE(drhoaug1)
1462    end if
1463    if (need_pawij10) then
1464      do kdir1=1,mdir1
1465        idir1=jdir1(kdir1)
1466        call paw_ij_free(paw_ij10(:,idir1))
1467      end do
1468      ABI_DEALLOCATE(e1kbfr_spin)
1469    end if
1470    ABI_DATATYPE_DEALLOCATE(paw_ij10)
1471 
1472 !  In case of parallelism, sum 1st-order density and occupation matrix over processors
1473    if (has_drho.and.xmpi_paral==1) then
1474 
1475 !    Accumulate 1st-order density
1476      call timab(48,1,tsec)
1477      bufsz=cplex*dtset%nfft*nspden*mdir1
1478      ABI_ALLOCATE(buffer,(bufsz))
1479      buffer(1:bufsz)=reshape(drho1wfr,(/bufsz/))
1480      call xmpi_sum(buffer,bufsz,spaceworld,ierr)
1481      drho1wfr(:,:,:)=reshape(buffer(1:bufsz),(/cplex*dtset%nfft,nspden,mdir1/))
1482      ABI_DEALLOCATE(buffer)
1483      call timab(48,2,tsec)
1484 
1485 !    Accumulate 1st-order PAW occupancies
1486      if (usepaw==1) then
1487        call pawrhoij_mpisum_unpacked(pawdrhoij1_unsym,spaceworld)
1488      end if
1489 
1490    end if
1491 
1492 !  Compute second part of overlap contribution (due to VHxc^(j2)(tild_n+hat_n))
1493    if (has_drho) then
1494 
1495      ABI_ALLOCATE(drhor1,(cplex*nfftf,nspden))
1496      ABI_ALLOCATE(dnhat1,(cplex*nfftf,nspden))
1497 
1498 !    LOOP OVER PERTURBATION DIRECTIONS
1499      do kdir1=1,mdir1
1500        idir1=jdir1(kdir1)
1501 
1502 !      Build and symmetrize 1st-order density change due to change of overlap
1503        ABI_ALLOCATE(drho1wfg,(2,dtset%nfft))
1504        call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,&
1505 &       nspden,nsppol,nsym1,dtset%paral_kgb,phnons1,drho1wfg,drho1wfr(:,:,idir1),&
1506 &       rprimd,symaf1,symrl1)
1507        if (dtset%pawstgylm/=0) then
1508          option=0
1509          call pawnhatfr(option,idir1,ipert1,my_natom,dtset%natom,nspden,dtset%ntypat,&
1510 &         pawang,pawfgrtab,pawrhoij,pawtab,rprimd,&
1511 &         mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
1512        end if
1513        call pawmkrho(1,arg,cplex,gprimd,idir1,indsy1,ipert1,&
1514 &       mpi_enreg,my_natom,dtset%natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,&
1515 &       pawfgr,pawfgrtab,-10001,pawdrhoij1(:,idir1),pawdrhoij1_unsym(:,idir1),pawtab,&
1516 &       dtset%qptn,drho1wfg,drho1wfr(:,:,idir1),drhor1,rprimd,symaf1,symrc1,dtset%typat,&
1517 &       ucvol,dtset%usewvl,xred,pawang_sym=pawang1,pawnhat=dnhat1,pawrhoij0=pawrhoij)
1518        ABI_DEALLOCATE(drho1wfg)
1519 
1520 !      Compute plane-wave contribution to overlap contribution
1521 !      This is subtle as it is a mix of Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
1522 !      Details:
1523 !      The VH(tild_nZc)^(1) term of Eq(79) is:
1524 !      <VH(tild_nZc)^(j2)|delta_tild_rho^(j1)>            = <vpsp1|drhor1-dnhat1>
1525 !      The first term of Eq(80) is:
1526 !      <VHxc^(j2)|delta_tild_rho^(j1)+delta_hat_rho^(j1)> = <vtrial1-vpsp1|drhor1>
1527 !      The addition of these two terms gives:
1528 !      <vtrial1|drhor1>-<vpsp1|dnhat1>
1529 !      And this is more subtle when usexcnhat=0
1530        call dotprod_vn(cplex,drhor1,dot1r,dot1i,nfftf,nfftot,nspden,2,vtrial1,ucvol)
1531        if (usexcnhat/=0) then
1532          call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,1   ,2,vpsp1,ucvol)
1533        else
1534          ABI_ALLOCATE(vtmp1,(cplex*nfftf,nspden))
1535          do ispden=1,nspden
1536            vtmp1(:,ispden)=vtrial1(:,ispden)-vhartr1(:)
1537          end do
1538          call dotprod_vn(cplex,dnhat1,dot2r,dot2i,nfftf,nfftot,nspden,2,vtmp1,ucvol)
1539          ABI_DEALLOCATE(vtmp1)
1540        end if
1541        dotr=dot1r-dot2r;doti=dot1i-dot2i
1542 
1543 !      Compute on-site contributions to overlap contribution
1544 !      (two last terms of Eq(80) of PRB 78, 035105 (2008)) [[cite:Audouze2008]]
1545 !      (note: Dij^(j2) and Vxc^(j2) are computed for ipert at first call)
1546        call pawdfptenergy(epawnst,ipert,ipert1,dtset%ixc,my_natom,dtset%natom,dtset%ntypat,&
1547 &       nzlmopt_ipert,nzlmopt_ipert1,paw_an,paw_an1,paw_ij1,pawang,dtset%pawprtvol,&
1548 &       pawrad,pawrhoij1,pawdrhoij1(:,idir1),pawtab,dtset%pawxcdev,dtset%xclevel,&
1549 &       mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
1550 
1551 !      Accumulate in 2nd-order matrix:
1552 !      Note: factor 2 (from d2E/dj1dj2=2E^(j1j2)) eliminated by factor 1/2
1553 !      has to take the complex conjugate because we want here Int[VHxc^(j1)^*.delta_rho^(j2)]
1554        dotr=dotr+epawnst(1);doti=-(doti+epawnst(2))
1555        d2ovl_drho(1,idir1,ipert1,idir,ipert)=elfd_fact*dotr
1556        d2ovl_drho(2,idir1,ipert1,idir,ipert)=elfd_fact*doti
1557 
1558      end do ! End loop over perturbation directions
1559 
1560 !    Free no more needed memory
1561      ABI_DEALLOCATE(drhor1)
1562      ABI_DEALLOCATE(dnhat1)
1563      ABI_DEALLOCATE(drho1wfr)
1564      do iatom=1,my_natom
1565        if (pawfgrtab(iatom)%nhatfr_allocated>0)  then
1566          ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr)
1567        end if
1568        pawfgrtab(iatom)%nhatfr_allocated=0
1569      end do
1570      if (paral_atom) then
1571        do kdir1=1,mdir1
1572          idir1=jdir1(kdir1)
1573          call pawrhoij_free(pawdrhoij1_unsym(:,idir1))
1574        end do
1575        ABI_DATATYPE_DEALLOCATE(pawdrhoij1_unsym)
1576      end if
1577      do kdir1=1,mdir1
1578        idir1=jdir1(kdir1)
1579        call pawrhoij_free(pawdrhoij1(:,idir1))
1580      end do
1581      ABI_DATATYPE_DEALLOCATE(pawdrhoij1)
1582    end if ! has_drho
1583 
1584 !  End loop over perturbations (j1)
1585  end do
1586 
1587 !Final deallocations
1588  ABI_DEALLOCATE(ch1c)
1589  if (usepaw==1.and.is_metal_or_qne0) then
1590    ABI_DEALLOCATE(cs1c)
1591  end if
1592  call destroy_hamiltonian(gs_hamkq)
1593 
1594 !In case of parallelism, sum over processors
1595  if (xmpi_paral==1)then
1596    call timab(161,1,tsec)
1597    call xmpi_barrier(spaceworld)
1598    call timab(161,2,tsec)
1599    ABI_ALLOCATE(buffer,(6*mpert*(1+usepaw)))
1600    buffer(1:6*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/6*mpert/))
1601    if (usepaw==1) buffer(6*mpert+1:6*mpert+6*mpert)=reshape(d2ovl(:,:,:,idir,ipert),(/6*mpert/))
1602    call timab(48,1,tsec)
1603    call xmpi_sum(buffer,6*mpert*(1+usepaw),spaceworld,ierr)
1604    call timab(48,2,tsec)
1605    d2nl (:,:,:,idir,ipert)=reshape(buffer(1:6*mpert),(/2,3,mpert/))
1606    if (usepaw==1) d2ovl(:,:,:,idir,ipert)=reshape(buffer(6*mpert+1:6*mpert+6*mpert),(/2,3,mpert/))
1607    ABI_DEALLOCATE(buffer)
1608  end if
1609 
1610 !Build complete d2ovl matrix
1611  if (usepaw==1) d2ovl(:,:,:,idir,ipert)=d2ovl(:,:,:,idir,ipert)+d2ovl_drho(:,:,:,idir,ipert)
1612 
1613  if (usepaw==1) then
1614    ABI_DEALLOCATE(d2ovl_drho)
1615  end if
1616 
1617 !Close the ddk WF files
1618  if (has_ddk_file) then
1619    do kdir1=1,mdir1
1620      idir1=jdir1(kdir1)
1621      if (ddkfil(idir1)/=0)then
1622        call wfk_close(ddks(idir1))
1623      end if
1624    end do
1625  end if
1626  ABI_DEALLOCATE(jpert1)
1627  ABI_DEALLOCATE(jdir1)
1628 
1629 !Symmetrize the phonons contributions, as was needed for the forces in a GS calculation
1630  ABI_ALLOCATE(work,(2,3,dtset%natom))
1631  do ipert1=1,dtset%natom
1632    do idir1=1,3
1633      work(:,idir1,ipert1)=d2nl(:,idir1,ipert1,idir,ipert)
1634    end do
1635  end do
1636  call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1)
1637  if (usepaw==1) then
1638    do ipert1=1,dtset%natom
1639      do idir1=1,3
1640        work(:,idir1,ipert1)=d2ovl(:,idir1,ipert1,idir,ipert)
1641      end do
1642    end do
1643    call dfpt_sygra(dtset%natom,d2ovl(:,:,:,idir,ipert),work,indsy1,ipert,nsym1,dtset%qptn,symrc1)
1644  end if
1645  ABI_DEALLOCATE(work)
1646 
1647 !In the case of the strain perturbation time-reversal symmetry will always
1648 !be true so imaginary part of d2nl will be must be set to zero here since
1649 !the symmetry-reduced kpt set will leave a non-zero imaginary part.
1650  if(ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
1651    d2nl(2,:,:,idir,ipert)=zero
1652    if (usepaw==1) d2ovl(2,:,:,idir,ipert)=zero
1653  else
1654    d2nl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero
1655    if (usepaw==1) d2ovl(2,:,dtset%natom+3:dtset%natom+4,idir,ipert)=zero
1656  end if
1657 
1658 
1659 !Symmetrize the strain perturbation contributions, as was needed for the stresses in a GS calculation
1660 !if (ipert==dtset%natom+3.or.ipert==dtset%natom+4)then
1661  if (nsym1>1) then
1662    ABI_ALLOCATE(work,(6,1,1))
1663    ii=0
1664    do ipert1=dtset%natom+3,dtset%natom+4
1665      do idir1=1,3
1666        ii=ii+1
1667        work(ii,1,1)=d2nl(1,idir1,ipert1,idir,ipert)
1668      end do
1669    end do
1670    call stresssym(gprimd,nsym1,work(:,1,1),symrc1)
1671    ii=0
1672    do ipert1=dtset%natom+3,dtset%natom+4
1673      do idir1=1,3
1674        ii=ii+1
1675        d2nl(1,idir1,ipert1,idir,ipert)=work(ii,1,1)
1676      end do
1677    end do
1678    if (usepaw==1) then
1679      ii=0
1680      do ipert1=dtset%natom+3,dtset%natom+4
1681        do idir1=1,3
1682          ii=ii+1
1683          work(ii,1,1)=d2ovl(1,idir1,ipert1,idir,ipert)
1684        end do
1685      end do
1686      call stresssym(gprimd,nsym1,work(:,1,1),symrc1)
1687      ii=0
1688      do ipert1=dtset%natom+3,dtset%natom+4
1689        do idir1=1,3
1690          ii=ii+1
1691          d2ovl(1,idir1,ipert1,idir,ipert)=work(ii,1,1)
1692        end do
1693      end do
1694      ABI_DEALLOCATE(work)
1695    end if
1696  end if
1697 !end if
1698 
1699 !Must also symmetrize the electric field perturbation response !
1700 !Note: d2ovl is not symetrized because it is zero for electric field perturbation
1701  if (has_ddk_file) then
1702    ABI_ALLOCATE(d2nl_elfd,(2,3))
1703 !  There should not be any imaginary part, but stay general (for debugging)
1704    d2nl_elfd (:,:)=d2nl(:,:,dtset%natom+2,idir,ipert)
1705    do ii=1,3
1706      sumelfd(:)=zero
1707      do ia=1,nsym1
1708        do jj=1,3
1709          if(symrl1(ii,jj,ia)/=0)then
1710            if(ddkfil(jj)==0)then
1711              blkflg(ii,dtset%natom+2,idir,ipert)=0
1712            end if
1713          end if
1714        end do
1715        symfact(1)=dble(symrl1(ii,1,ia))
1716        symfact(2)=dble(symrl1(ii,2,ia))
1717        symfact(3)=dble(symrl1(ii,3,ia))
1718        sumelfd(:)=sumelfd(:)+symfact(1)*d2nl_elfd(:,1) &
1719 &       +symfact(2)*d2nl_elfd(:,2)+symfact(3)*d2nl_elfd(:,3)
1720      end do
1721      d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1)
1722    end do
1723    ABI_DEALLOCATE(d2nl_elfd)
1724  end if
1725 
1726 !Overlap: store the diagonal part of the matrix in the
1727 !         2nd-order energy non-stationnary expression
1728  eovl1=zero;if (usepaw==1) eovl1=d2ovl(1,idir,ipert,idir,ipert)
1729 
1730  call destroy_mpi_enreg(mpi_enreg_seq)
1731  call timab(566,2,tsec)
1732 
1733  DBG_EXIT("COLL")
1734 
1735 end subroutine dfpt_nstpaw

ABINIT/dfpt_nstwf [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nstwf

FUNCTION

 This routine computes the non-local contribution to the
 2DTE matrix elements, in the non-stationary formulation
 Only for norm-conserving pseudopotentials (no PAW)

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (XG,AR,MB,MVer,MT, MVeithen)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  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.
  ddkfil(3)=unit numbers for the three possible ddk files for ipert1
       equal to 0 if no dot file is available for this direction
  dtset <type(dataset_type)>=all input variables for this dataset
  eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
  eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  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
  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
  istwf_k=parameter that describes 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
  kpt(3)=reduced coordinates of k point
  kpq(3)=reduced coordinates of k+q point
  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).
  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
  nsppol=1 for unpolarized, 2 for spin-polarized
  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)
  ddks(3)<wfk_t>=struct info for for the three possible DDK files for ipert1
  wtk_k=weight assigned to the k point.
  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

  d2bbb_k(2,3,mband,mband*prtbbb)=band by band decomposition of the second
   order derivatives, for the present k point, and perturbation idir, ipert
  d2nl_k(2,3,mpert)=non-local contributions to
   non-stationary 2DTE, for the present k point, and perturbation idir, ipert

TODO

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

PARENTS

      dfpt_nstdy

CHILDREN

      destroy_rf_hamiltonian,dotprod_g,gaugetransfo,getgh1c
      init_rf_hamiltonian,load_k_hamiltonian,load_k_rf_hamiltonian
      load_kprime_hamiltonian,mkffnl,mkkpg,timab,wfk_read_bks

SOURCE

1811 subroutine dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,&
1812 &                 icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpt,kpq,&
1813 &                 mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,&
1814 &                 occ_k,psps,rmet,ddks,wtk_k,ylm,ylm1)
1815 
1816 
1817  use defs_basis
1818  use defs_datatypes
1819  use defs_abitypes
1820  use m_abicore
1821  use m_cgtools
1822  use m_hamiltonian
1823  use m_errors
1824  use m_wfk
1825  use m_xmpi
1826 
1827  use m_time,    only : timab
1828  use m_pawcprj, only : pawcprj_type
1829  use m_kg,      only : mkkpg
1830  use m_mkffnl,  only : mkffnl
1831  use m_getgh1c, only : getgh1c
1832 
1833 !This section has been created automatically by the script Abilint (TD).
1834 !Do not modify the following lines by hand.
1835 #undef ABI_FUNC
1836 #define ABI_FUNC 'dfpt_nstwf'
1837 !End of the abilint section
1838 
1839  implicit none
1840 
1841 !Arguments ------------------------------------
1842 !scalars
1843  integer,intent(in) :: icg,icg1,idir,ikpt,ipert,isppol,istwf_k
1844  integer,intent(in) :: mkmem,mk1mem,mpert,mpw,mpw1,nsppol
1845  integer,intent(inout) :: nband_k,npw1_k,npw_k
1846  real(dp),intent(in) :: wtk_k
1847  type(MPI_type),intent(in) :: mpi_enreg
1848  type(dataset_type),intent(in) :: dtset
1849  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
1850  type(pseudopotential_type),intent(in) :: psps
1851 !arrays
1852  integer,intent(in) :: ddkfil(3),kg1_k(3,npw1_k)
1853  integer,intent(in) :: kg_k(3,npw_k)
1854  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
1855  real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*nsppol)
1856  real(dp),intent(in) :: eig_k(dtset%mband*nsppol),kpt(3),kpq(3),occ_k(nband_k),rmet(3,3)
1857  real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
1858  real(dp),intent(in) :: ylm1(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)
1859  real(dp),intent(inout) :: eig1_k(2*nsppol*dtset%mband**2)
1860  real(dp),intent(out) :: d2bbb_k(2,3,dtset%mband,dtset%mband*dtset%prtbbb)
1861  real(dp),intent(inout) :: d2nl_k(2,3,mpert)
1862  type(wfk_t),intent(inout) :: ddks(3)
1863 
1864 !Local variables-------------------------------
1865 !scalars
1866  integer :: berryopt,dimffnl,dimffnl1,dimph3d
1867  integer :: iband,ider,idir1,ipert1,ipw,jband,nband_kocc,nkpg,nkpg1 !ierr,ii
1868  integer :: npw_disk,nsp,optlocal,optnl,opt_gvnl1,sij_opt,tim_getgh1c,usevnl
1869  logical :: ddk
1870  real(dp) :: aa,dot1i,dot1r,dot2i,dot2r,dot_ndiagi,dot_ndiagr,doti,dotr,lambda
1871  character(len=500) :: msg
1872  type(rf_hamiltonian_type) :: rf_hamkq
1873 !arrays
1874  integer :: ik_ddks(3)
1875  real(dp) :: dum_grad_berry(1,1),dum_gvnl1(1,1),dum_gs1(1,1),dum_ylmgr(1,3,1),tsec(2)
1876  real(dp),allocatable :: cg_k(:,:),cwave0(:,:),cwavef(:,:),cwavef_da(:,:)
1877  real(dp),allocatable :: cwavef_db(:,:),dkinpw(:),eig2_k(:),ffnl1(:,:,:,:),ffnlk(:,:,:,:)
1878  real(dp),allocatable :: gvnl1(:,:),kinpw1(:),kpg1_k(:,:),kpg_k(:,:),ph3d(:,:,:)
1879  type(pawcprj_type),allocatable :: dum_cwaveprj(:,:)
1880 
1881 ! *********************************************************************
1882 
1883  DBG_ENTER("COLL")
1884 
1885 !Not valid for PAW
1886  if (psps%usepaw==1) then
1887    msg='  This routine cannot be used for PAW (use pawnst3 instead) !'
1888    MSG_BUG(msg)
1889  end if
1890 
1891 !Keep track of total time spent in dfpt_nstwf
1892  call timab(102,1,tsec)
1893  tim_getgh1c=2
1894 
1895 !Miscelaneous inits
1896  ABI_DATATYPE_ALLOCATE(dum_cwaveprj,(0,0))
1897  ddk=(ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11)
1898 
1899 !Additional allocations
1900  if (.not.ddk) then
1901    ABI_ALLOCATE(dkinpw,(npw_k))
1902    ABI_ALLOCATE(kinpw1,(npw1_k))
1903    kinpw1=zero;dkinpw=zero
1904  else
1905    ABI_ALLOCATE(dkinpw,(0))
1906    ABI_ALLOCATE(kinpw1,(0))
1907  end if
1908  ABI_ALLOCATE(gvnl1,(2,npw1_k*dtset%nspinor))
1909  ABI_ALLOCATE(eig2_k,(2*nsppol*dtset%mband**2))
1910  ABI_ALLOCATE(cwave0,(2,npw_k*dtset%nspinor))
1911  ABI_ALLOCATE(cwavef,(2,npw1_k*dtset%nspinor))
1912 
1913 !Compute (k+G) vectors
1914  nkpg=0;if (.not.ddk) nkpg=3*gs_hamkq%nloalg(3)
1915  ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
1916  if (nkpg>0) then
1917    call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
1918  end if
1919 
1920 !Compute (k+q+G) vectors
1921  nkpg1=0;if (.not.ddk) nkpg1=3*gs_hamkq%nloalg(3)
1922  ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1))
1923  if (nkpg1>0) then
1924    call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k)
1925  end if
1926 
1927 !Compute nonlocal form factors ffnl at (k+G)
1928  dimffnl=0;if (.not.ddk) dimffnl=1
1929  ABI_ALLOCATE(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
1930  if (.not.ddk) then
1931    ider=0
1932    call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamkq%gmet,&
1933 &   gs_hamkq%gprimd,ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,&
1934 &   psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,&
1935 &   psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,dum_ylmgr)
1936  end if
1937 
1938 !Compute nonlocal form factors ffnl1 at (k+q+G)
1939  dimffnl1=0;if (.not.ddk) dimffnl1=1
1940  ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
1941  if (.not.ddk) then
1942    ider=0
1943    call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,&
1944 &   gs_hamkq%gprimd,ider,ider,psps%indlmn,kg1_k,kpg1_k,kpq,&
1945 &   psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,npw1_k,psps%ntypat,&
1946 &   psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1,dum_ylmgr)
1947  end if
1948 
1949 !Load k-dependent part in the Hamiltonian datastructure
1950  call load_k_hamiltonian(gs_hamkq,kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,&
1951 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk)
1952 
1953 !Load k+q-dependent part in the Hamiltonian datastructure
1954 !    Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
1955  dimph3d=0;if (.not.ddk) dimph3d=gs_hamkq%matblk
1956  ABI_ALLOCATE(ph3d,(2,npw1_k,dimph3d))
1957  call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
1958 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,&
1959 & ph3d_kp=ph3d,compute_ph3d=(.not.ddk))
1960 
1961 !Load k-dependent part in the 1st-order Hamiltonian datastructure
1962  call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw)
1963 
1964 !Take care of the npw and kg records
1965 !NOTE : one should be able to modify the rwwf routine to take care
1966 !of the band parallelism, which is not the case yet ...
1967  ik_ddks = 0
1968  do idir1=1,3
1969    if (ddkfil(idir1)/=0)then
1970 !    Read npw record
1971      nsp=dtset%nspinor
1972      ik_ddks(idir1) = wfk_findk(ddks(idir1), kpt)
1973      ABI_CHECK(ik_ddks(idir1) /= -1, "Cannot find kpt")
1974      npw_disk = ddks(idir1)%hdr%npwarr(ik_ddks(idir1))
1975      if (npw_k /= npw_disk) then
1976        write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')&
1977 &       'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',idir,ch10,&
1978 &       'the number of plane waves in the ddk file is equal to', npw_disk,ch10,&
1979 &       'while it should be ',npw_k
1980        MSG_BUG(msg)
1981      end if
1982    end if
1983  end do
1984 
1985  if (ipert==dtset%natom+1) then
1986    nband_kocc = 0
1987    do iband = 1,nband_k
1988      if (abs(occ_k(iband)) > tol8) nband_kocc = nband_kocc + 1
1989      nband_kocc = max (nband_kocc, 1)
1990    end do
1991  end if
1992 
1993  if(dtset%prtbbb==1)then
1994    ABI_ALLOCATE(cwavef_da,(2,npw1_k*dtset%nspinor))
1995    ABI_ALLOCATE(cwavef_db,(2,npw1_k*dtset%nspinor))
1996    ABI_ALLOCATE(cg_k,(2,npw_k*dtset%nspinor*nband_k))
1997    if ((ipert == dtset%natom + 1).or.(ipert <= dtset%natom).or. &
1998 &   (ipert == dtset%natom + 2).or.(ipert == dtset%natom + 5)) then
1999      cg_k(:,:) = cg(:,1+icg:icg+nband_k*npw_k*dtset%nspinor)
2000    end if
2001    d2bbb_k(:,:,:,:) = zero
2002  end if
2003 
2004 !Loop over bands
2005  do iband=1,nband_k
2006 
2007    if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
2008 
2009 !  Read ground-state wavefunctions
2010    if (dtset%prtbbb==0 .or. ipert==dtset%natom+2) then
2011      cwave0(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg)
2012    else    ! prtbbb==1 and ipert<=natom , already in cg_k
2013      cwave0(:,:)=cg_k(:,1+(iband-1)*npw_k*dtset%nspinor:iband*npw_k*dtset%nspinor)
2014    end if
2015 
2016 !  Get first-order wavefunctions
2017    cwavef(:,:)=cg1(:,1+(iband-1)*npw1_k*dtset%nspinor+icg1:iband*npw1_k*dtset%nspinor+icg1)
2018 
2019 !  In case non ddk perturbation
2020    if (ipert /= dtset%natom + 1) then
2021 
2022      do ipert1=1,mpert
2023 
2024        if( ipert1<=dtset%natom .or. ipert1==dtset%natom+2 )then
2025 
2026 !        Initialize data for NL 1st-order hamiltonian
2027          call init_rf_hamiltonian(1,gs_hamkq,ipert1,rf_hamkq)
2028 
2029          if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)) &
2030 &         .and.(ipert1 == dtset%natom+2).and. dtset%prtbbb==1) then
2031            call gaugetransfo(cg_k,cwavef,cwavef_db,eig_k,eig1_k,iband,nband_k, &
2032 &           dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k)
2033            cwavef(:,:) = cwavef_db(:,:)
2034          end if
2035 
2036 !        Define the direction along which to move the atom :
2037 !        the polarisation (ipert1,idir1) is refered as j1.
2038          do idir1=1,3
2039            if (ipert1<=dtset%natom.or.(ipert1==dtset%natom+2.and.ddkfil(idir1)/=0)) then
2040 
2041 !            Get |Vnon-locj^(1)|u0> :
2042 !            First-order non-local, applied to zero-order wavefunction
2043 !            This routine gives MINUS the non-local contribution
2044 
2045 !            ==== Atomic displ. perturbation
2046              if( ipert1<=dtset%natom )then
2047                lambda=eig_k((isppol-1)*nband_k+iband)
2048                berryopt=1;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0
2049                call getgh1c(berryopt,cwave0,dum_cwaveprj,gvnl1,dum_grad_berry,&
2050 &               dum_gs1,gs_hamkq,dum_gvnl1,idir1,ipert1,lambda,mpi_enreg,optlocal,&
2051 &               optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
2052 
2053 !              ==== Electric field perturbation
2054              else if( ipert1==dtset%natom+2 )then
2055                ! TODO: Several tests fail here ifdef HAVE_MPI_IO_DEFAULT
2056                ! The problem is somehow related to the use of MPI-IO file views!.
2057                call wfk_read_bks(ddks(idir1), iband, ik_ddks(idir1), isppol, xmpio_single, cg_bks=gvnl1, &
2058                eig1_bks=eig2_k(1+(iband-1)*2*nband_k:))
2059                  !eig1_bks=eig2_k(1+(iband-1)*2*nband_k:2*iband*nband_k))
2060                !write(777,*)"eig2_k, gvnl1 for band: ",iband,", ikpt: ",ikpt
2061                !do ii=1,2*nband_k
2062                !  write(777,*)eig2_k(ii+(iband-1))
2063                !end do
2064                !write(777,*)gvnl1
2065 
2066 !              In case of band-by-band,
2067 !              construct the first-order wavefunctions in the diagonal gauge
2068                if (((ipert <= dtset%natom).or.(ipert == dtset%natom + 2)).and.(dtset%prtbbb==1)) then
2069                  call gaugetransfo(cg_k,gvnl1,cwavef_da,eig_k,eig2_k,iband,nband_k, &
2070 &                 dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k)
2071                  gvnl1(:,:) = cwavef_da(:,:)
2072                end if
2073 !              Multiplication by -i
2074                do ipw=1,npw1_k*dtset%nspinor
2075                  aa=gvnl1(1,ipw)
2076                  gvnl1(1,ipw)=gvnl1(2,ipw)
2077                  gvnl1(2,ipw)=-aa
2078                end do
2079              end if
2080 
2081 !            MVeithen 021212 :
2082 !            1) Case ipert1 = natom + 2 and ipert = natom + 2:
2083 !            the second derivative of the energy with respect to an electric
2084 !            field is computed from Eq. (38) of X. Gonze, PRB 55 ,10355 (1997) [[cite:Gonze1997a]].
2085 !            The evaluation of this formula needs the operator $i \frac{d}{dk}.
2086 !            2) Case ipert1 = natom + 2 and ipert < natom:
2087 !            the computation of the Born effective charge tensor uses
2088 !            the operator $-i \frac{d}{dk}.
2089              if (ipert==dtset%natom+2) gvnl1(:,:) = -gvnl1(:,:)
2090 
2091 !            <G|Vnl1|Cnk> is contained in gvnl1
2092 !            construct the matrix element (<uj2|vj1|u0>)complex conjug and add it to the 2nd-order matrix
2093              call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,cwavef,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2094              d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*two*dotr
2095              d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*two*doti
2096 
2097 !            Band by band decomposition of the Born effective charges
2098 !            calculated from a phonon perturbation
2099              if(dtset%prtbbb==1)then
2100                d2bbb_k(1,idir1,iband,iband) = wtk_k*occ_k(iband)*two*dotr
2101                d2bbb_k(2,idir1,iband,iband) = -one*wtk_k*occ_k(iband)*two*doti
2102              end if
2103 
2104            end if
2105          end do
2106 
2107          call destroy_rf_hamiltonian(rf_hamkq)
2108        end if
2109      end do
2110    end if     ! ipert /= natom +1
2111 
2112 !  Compute the localization tensor
2113 
2114    if (ipert==dtset%natom+1) then
2115 
2116      ipert1=dtset%natom+1
2117      if(dtset%prtbbb==1)then
2118        call gaugetransfo(cg_k,cwavef,cwavef_db,eig_k,eig1_k,iband,nband_k, &
2119 &       dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k)
2120        cwavef(:,:) = cwavef_db(:,:)
2121      end if
2122 
2123      do idir1 = 1,3
2124        eig2_k(:) = zero
2125        gvnl1(:,:) = zero
2126        if (idir == idir1) then
2127          gvnl1(:,:) = cwavef(:,:)
2128          eig2_k(:) = eig1_k(:)
2129        else
2130          if (ddkfil(idir1) /= 0) then
2131            call wfk_read_bks(ddks(idir1), iband, ik_ddks(idir1), isppol, xmpio_single, cg_bks=gvnl1, &
2132            eig1_bks=eig2_k(1+(iband-1)*2*nband_k:))
2133              !eig1_bks=eig2_k(1+(iband-1)*2*nband_k:2*iband*nband_k))
2134            !write(778,*)"eig2_k, gvnl1 for band: ",iband,", ikpt: ",ikpt
2135            !do ii=1,2*nband_k
2136            !  write(778,*)eig2_k(ii+(iband-1))
2137            !end do
2138            !write(778,*)gvnl1
2139 
2140            if(dtset%prtbbb==1)then
2141              call gaugetransfo(cg_k,gvnl1,cwavef_da,eig_k,eig2_k,iband,nband_k, &
2142 &             dtset%mband,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k)
2143              gvnl1(:,:) = cwavef_da(:,:)
2144            end if
2145 
2146          end if    !ddkfil(idir1)
2147        end if    !idir == idir1
2148 
2149 !      <G|du/dqa> is contained in gvnl1 and <G|du/dqb> in cwavef
2150 !      construct the matrix elements <du/dqa|du/dqb> -> dot
2151 !      <u|du/dqa> -> dot1
2152 !      <du/dqb|u> -> dot2
2153 !      and add them to the 2nd-order matrix
2154 
2155        call dotprod_g(dotr,doti,istwf_k,npw1_k*dtset%nspinor,2,gvnl1,cwavef,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2156        d2nl_k(1,idir1,ipert1)=d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*dotr/(nband_kocc*two)
2157        d2nl_k(2,idir1,ipert1)=d2nl_k(2,idir1,ipert1)+wtk_k*occ_k(iband)*doti/(nband_kocc*two)
2158 
2159 
2160 !      XG 020216 : Marek, could you check the next forty lines
2161 !      In the parallel gauge, dot1 and dot2 vanishes
2162        if(dtset%prtbbb==1)then
2163          d2bbb_k(1,idir1,iband,iband)=d2bbb_k(1,idir1,iband,iband)+dotr
2164          d2bbb_k(2,idir1,iband,iband)=d2bbb_k(2,idir1,iband,iband)+doti
2165          dot_ndiagr=zero ; dot_ndiagi=zero
2166          do jband = 1,nband_k              !compute dot1 and dot2
2167            if (abs(occ_k(jband)) > tol8) then
2168              dot1r=zero ; dot1i=zero
2169              dot2r=zero ; dot2i=zero
2170              cwave0(:,:)=cg_k(:,1+(jband-1)*npw_k*dtset%nspinor:jband*npw_k*dtset%nspinor)
2171 
2172              call dotprod_g(dot1r,dot1i,istwf_k,npw1_k*dtset%nspinor,2,cwave0,gvnl1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2173              call dotprod_g(dot2r,dot2i,istwf_k,npw1_k*dtset%nspinor,2,cwavef,cwave0,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
2174 
2175              dot_ndiagr = dot_ndiagr + dot1r*dot2r - dot1i*dot2i
2176              dot_ndiagi = dot_ndiagi + dot1r*dot2i + dot1i*dot2r
2177              d2bbb_k(1,idir1,iband,jband) = d2bbb_k(1,idir1,iband,jband) - &
2178 &             (dot1r*dot2r - dot1i*dot2i)
2179              d2bbb_k(2,idir1,iband,jband) = d2bbb_k(2,idir1,iband,jband) - &
2180 &             (dot1r*dot2i + dot1i*dot2r)
2181            end if  ! occ_k
2182          end do !jband
2183          d2bbb_k(:,idir1,iband,:)=d2bbb_k(:,idir1,iband,:)*wtk_k*occ_k(iband)*half
2184          d2nl_k(1,idir1,ipert1)= &
2185 &         d2nl_k(1,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagr/(nband_kocc*two)
2186          d2nl_k(2,idir1,ipert1)=&
2187 &         d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*dot_ndiagi/(nband_kocc*two)
2188        end if ! prtbbb==1
2189 
2190      end do  ! idir1
2191    end if   ! Compute localization tensor, ipert=natom+1
2192 
2193 !  End loop over bands
2194  end do
2195 
2196 !Final deallocations
2197  ABI_DEALLOCATE(cwave0)
2198  ABI_DEALLOCATE(cwavef)
2199  ABI_DEALLOCATE(eig2_k)
2200  ABI_DEALLOCATE(gvnl1)
2201  ABI_DEALLOCATE(ffnlk)
2202  ABI_DEALLOCATE(ffnl1)
2203  ABI_DEALLOCATE(dkinpw)
2204  ABI_DEALLOCATE(kinpw1)
2205  ABI_DEALLOCATE(kpg_k)
2206  ABI_DEALLOCATE(kpg1_k)
2207  ABI_DEALLOCATE(ph3d)
2208  ABI_DATATYPE_DEALLOCATE(dum_cwaveprj)
2209  if(dtset%prtbbb==1)  then
2210    ABI_DEALLOCATE(cg_k)
2211    ABI_DEALLOCATE(cwavef_da)
2212    ABI_DEALLOCATE(cwavef_db)
2213  end if
2214 
2215  call timab(102,2,tsec)
2216 
2217  DBG_EXIT("COLL")
2218 
2219   contains

ABINIT/gaugetransfo [ Functions ]

[ Top ] [ Functions ]

NAME

 gaugetransfo

FUNCTION

 This routine allows the passage from the parallel-transport gauge
 to the diagonal gauge for the first-order wavefunctions

INPUTS

  cg_k(2,mpw*nspinor*mband*nsppol)=planewave coefficients of wavefunctions
                                   for a particular k point.
  cwavef(2,npw1_k*nspinor)=first order wavefunction for a particular k point
                           in the parallel gauge
  eig_k(mband*nsppol)=GS eigenvalues at k (hartree)
  eig1_k(2*nsppol*mband**2)=matrix of first-order eigenvalues (hartree)
  iband=band index of the 1WF for which the transformation has to be applied
  mband=maximum number of bands
  nband_k=number of bands for this k point
  npw_k=maximum dimensioned size of npw or wfs at k
  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

OUTPUT

  cwavef_d(2,npw1_k*nspinor)=first order wavefunction for a particular k point
                             in the diagonal gauge

PARENTS

      dfpt_nstwf

CHILDREN

SOURCE

2257 subroutine gaugetransfo(cg_k,cwavef,cwavef_d,eig_k,eig1_k,iband,nband_k, &
2258 &                      mband,npw_k,npw1_k,nspinor,nsppol,occ_k)
2259 
2260 
2261 !This section has been created automatically by the script Abilint (TD).
2262 !Do not modify the following lines by hand.
2263 #undef ABI_FUNC
2264 #define ABI_FUNC 'gaugetransfo'
2265 !End of the abilint section
2266 
2267  implicit none
2268 
2269 !Arguments ------------------------------------
2270 !scalars
2271  integer,intent(in) :: iband,mband,nband_k,npw1_k,npw_k,nspinor,nsppol
2272 !arrays
2273  real(dp),intent(in) :: cg_k(2,npw_k*nspinor*nband_k),cwavef(2,npw1_k*nspinor)
2274  real(dp),intent(in) :: eig1_k(2*nsppol*mband**2),eig_k(mband*nsppol)
2275  real(dp),intent(in) :: occ_k(nband_k)
2276  real(dp),intent(out) :: cwavef_d(2,npw1_k*nspinor)
2277 
2278 !Local variables-------------------------------
2279 !tolerance for non degenerated levels
2280 !scalars
2281  integer :: jband
2282  real(dp),parameter :: etol=1.0d-3
2283 !arrays
2284  real(dp) :: cwave0(2,npw1_k*nspinor),eig1(2)
2285 
2286 ! *********************************************************************
2287 
2288 !DEBUG
2289 !write(100,*) 'gaugetransfo: ',iband
2290 !ENDDEBUG
2291 
2292    cwavef_d(:,:) = cwavef(:,:)
2293 
2294    do jband = 1,nband_k !loop over bands
2295 
2296      if ((abs(eig_k(iband)-eig_k(jband)) > etol).and.(abs(occ_k(jband)) > tol8 )) then
2297 
2298        cwave0(:,:) = cg_k(:,1+(jband-1)*npw_k*nspinor:jband*npw_k*nspinor)
2299 
2300        eig1(1) = eig1_k(2*jband-1+(iband-1)*2*nband_k)
2301        eig1(2) = eig1_k(2*jband +(iband-1)*2*nband_k)
2302 
2303        cwavef_d(1,:)=cwavef_d(1,:) &
2304 &       - (eig1(1)*cwave0(1,:)-eig1(2)*cwave0(2,:))/(eig_k(jband)-eig_k(iband))
2305        cwavef_d(2,:)=cwavef_d(2,:) &
2306 &       - (eig1(1)*cwave0(2,:)+eig1(2)*cwave0(1,:))/(eig_k(jband)-eig_k(iband))
2307 
2308      end if
2309 
2310    end do    !loop over bands
2311 
2312   end subroutine gaugetransfo

ABINIT/m_dfpt_nstwf [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_nstwf

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group ()
  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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_dfpt_nstwf
27 
28  implicit none
29 
30  private