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)

NOTES

     delta_u^(j1)=-1/2 Sum_{j}[<u0_k+q_j|S^(j1)|u0_k_i>.|u0_k+q_j>]

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

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