TABLE OF CONTENTS


ABINIT/dfptnl_pert [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptnl_pert

FUNCTION

 Compute the linear response part to the 3dte. The main inputs are :
   - GS WFs and Hamiltonian (cg,gs_hamkq)
   - 1st-order WFs for three perturbations i1pert/i1dir,i2pert/i2dir,i3pert/i3dir (cg1,cg2,cg3)
   - 1st-order potentials for i2pert (vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert)
   - 1st-order WFs DDK,DDE and 2nd-order WF DKDE (ddk_f)

COPYRIGHT

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

  atindx(natom)=index table for atoms (see gstate.f)
  cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave
                                          coefficients of wavefunctions
  cg1 = first derivative of cg with respect the perturbation i1pert
  cg2 = first derivative of cg with respect the perturbation i2pert
  cg3 = first derivative of cg with respect the perturbation i3pert
  cplex= if 1, real space 1-order functions on FFT grid are REAL,
          if 2, COMPLEX
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  i1dir,i2dir,i3dir=directions of the corresponding perturbations
  i1pert,i2pert,i3pert = type of perturbation that has to be computed
  kg(3,mpw*mkmem)=reduced planewave coordinates
  mband = maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem = maximum number of k points which can fit in core memory
  mk1mem = maximum number of k points for first-order WF
           which can fit in core memory
  mpert =maximum number of ipert
  mpi_enreg=MPI-parallelisation information
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw   = maximum number of planewaves in basis sphere (large number)
  natom = number of atoms in unit cell
  nattyp(ntypat)= # atoms of each type.
  nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
  nfftotf=total number of real space fine grid points
  ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90)
  nkpt = number of k points
  nk3xc=second dimension of the array k3xc
  nspden = number of spin-density components
  nspinor = number of spinorial components of the wavefunctions
  nsppol = number of channels for spin-polarization (1 or 2)
  nsym1=number of symmetry elements in space group consistent with the perturbation
  npwarr(nkpt) = array holding npw for each k point
  occ(mband*nkpt*nsppol) = occupation number for each band and k
  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
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  pawrhoij0(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
  pawrhoij1_i1pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i1pert)
  pawrhoij1_i2pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i2pert)
  pawrhoij1_i3pert(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data (i3pert)
  paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
  paw_an1_i2pert(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh (i2pert)
  paw_ij1_i2pert(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels (i2pert)
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  psps <type(pseudopotential_type)> = variables related to pseudopotentials
  rho1r1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i1pert)
  rho1r2(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i2pert)
  rho1r3(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3 (i3pert)
  rprimd(3,3) = dimensional primitive translations (bohr)
  symaf1(nsym)=(anti)ferromagnetic part of symmetry operations
  symrc1(3,3,nsym)=symmetries of group in terms of operations on reciprocal space primitive translations
  ucvol=volume of the unit cell
  vtrial(nfftf,nspden)=GS Vtrial(r).
  vhartr1_i2pert(cplex*nfftf,nspden)=firs-order hartree potential
  vtrial1_i2pert(cplex*nfft,nspden)=firs-order local potential
  vxc1_i2pert(cplex*nfft,nspden)=firs-order exchange-correlation potential
  ddk_f = wf files
  xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc (i1pert)
  xccc3d2(cplex*n3xccc)=3D change in core charge density, see n3xccc (i2pert)
  xccc3d3(cplex*n3xccc)=3D change in core charge density, see n3xccc (i3pert)
  xred(3,natom) = reduced atomic coordinates

OUTPUT

  d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
                                    = \sum_{i=1}^9 d3etot_i
  d3etot_1(2,3,mpert,3,mpert,3,mpert) = 1st term of d3etot
  d3etot_2(2,3,mpert,3,mpert,3,mpert) = 2nd term of d3etot
  d3etot_3(2,3,mpert,3,mpert,3,mpert) = 3rd term of d3etot
  d3etot_4(2,3,mpert,3,mpert,3,mpert) = 4th term of d3etot
  d3etot_5(2,3,mpert,3,mpert,3,mpert) = 5th term of d3etot
  d3etot_6(2,3,mpert,3,mpert,3,mpert) = 6th term of d3etot
  d3etot_7(2,3,mpert,3,mpert,3,mpert) = 7th term of d3etot
  d3etot_8(2,3,mpert,3,mpert,3,mpert) = 8th term of d3etot
  d3etot_9(2,3,mpert,3,mpert,3,mpert) = 9th term of d3etot

SIDE EFFECTS

  TO DO!

PARENTS

      dfptnl_loop

CHILDREN

      destroy_hamiltonian,dotprod_g,fftpac,fourwf,init_hamiltonian
      load_k_hamiltonian,mkffnl,mkkpg,nonlop,status,xmpi_sum

SOURCE

 155 subroutine dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,i2dir,i3dir,&
 156 & i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,&
 157 & nspden,nspinor,nsppol,nsym1,npwarr,occ,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,&
 158 & pawrhoij0,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,&
 159 & paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,rprimd,symaf1,symrc1,&
 160 & ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,&
 161 & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9)
 162 
 163  use defs_basis
 164  use defs_datatypes
 165  use defs_abitypes
 166  use m_abicore
 167  use m_wffile
 168  use m_wfk
 169  use m_xmpi
 170  use m_hamiltonian
 171  use m_errors
 172  use m_rf2
 173  use m_kg
 174 
 175  use m_dtfil,      only : status
 176  use m_cgtools,    only : dotprod_g,sqnorm_g,dotprod_vn
 177  use m_pawang,     only : pawang_type
 178  use m_pawfgrtab,  only : pawfgrtab_type
 179  use m_pawrad,     only : pawrad_type
 180  use m_pawtab,     only : pawtab_type
 181  use m_pawcprj,    only : pawcprj_type, pawcprj_alloc, pawcprj_free
 182  use m_paw_ij,     only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify,paw_ij_reset_flags
 183  use m_pawdij,     only : pawdijfr
 184  use m_pawfgr,     only : pawfgr_type
 185  use m_pawrhoij,   only : pawrhoij_type, pawrhoij_alloc , pawrhoij_nullify, pawrhoij_free,&
 186 &                         pawrhoij_init_unpacked, pawrhoij_mpisum_unpacked
 187  use m_paw_an,     only : paw_an_type
 188  use m_paw_mkrho,  only : pawmkrho
 189  use m_paw_nhat,   only : pawnhatfr
 190  use m_paw_dfpt,   only : pawdfptenergy
 191  use m_paw_dfptnl, only : paw_dfptnl_accrhoij,paw_dfptnl_energy
 192  use m_initylmg,   only : initylmg
 193  use m_mkffnl,     only : mkffnl
 194  use m_getghc,     only : getgsc
 195  use m_getgh1c,    only : rf_transgrid_and_pack
 196  use m_mpinfo,     only : proc_distrb_cycle
 197  use m_nonlop,     only : nonlop
 198  use m_fourier_interpol, only : transgrid
 199  use m_cgprj,     only : getcprj
 200 
 201 !This section has been created automatically by the script Abilint (TD).
 202 !Do not modify the following lines by hand.
 203 #undef ABI_FUNC
 204 #define ABI_FUNC 'dfptnl_pert'
 205 !End of the abilint section
 206 
 207  implicit none
 208 
 209 !Arguments ------------------------------------
 210 !scalars
 211  integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft
 212  integer,intent(in) :: mk1mem,mkmem,mpert,mpsang,mpw,natom,nfftf,nfftotf,nkpt,nspden,nsym1
 213  integer,intent(in) :: nk3xc,nspinor,nsppol
 214  real(dp),intent(in) :: ucvol
 215  type(MPI_type),intent(inout) :: mpi_enreg
 216  type(datafiles_type),intent(in) :: dtfil
 217  type(dataset_type),intent(in) :: dtset
 218  type(pseudopotential_type),intent(in) :: psps
 219  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
 220  type(pawang_type),intent(inout) :: pawang,pawang1
 221  type(pawfgr_type),intent(in) :: pawfgr
 222  type(wfk_t),intent(inout) :: ddk_f(5)
 223 
 224 !arrays
 225  integer,intent(in) :: atindx(natom),kg(3,mpw*mkmem),nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt)
 226  integer,intent(in) :: indsy1(4,nsym1,dtset%natom),symaf1(nsym1),symrc1(3,3,nsym1)
 227  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
 228  real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
 229  real(dp),intent(in) :: cg2(2,mpw*nspinor*mband*mk1mem*nsppol)
 230  real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
 231  real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)
 232  real(dp),intent(in) :: k3xc(nfftf,nk3xc)
 233  real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom)
 234  real(dp),intent(in) :: rho1r1(cplex*nfftf,dtset%nspden),rho2r1(cplex*nfftf,dtset%nspden)
 235  real(dp),intent(in) :: rho3r1(cplex*nfftf,dtset%nspden),rprimd(3,3)
 236  real(dp),intent(in) :: vtrial(cplex*nfftf,nspden)
 237  real(dp),intent(in) :: xccc3d1(cplex*nfftf),xccc3d2(cplex*nfftf),xccc3d3(cplex*nfftf),xred(3,natom)
 238  real(dp),intent(in) :: vxc1_i2pert(cplex*nfftf,nspden),vhartr1_i2pert(cplex*nfftf)
 239  real(dp),intent(inout) :: vtrial1_i2pert(cplex*nfftf,nspden),d3etot(2,3,mpert,3,mpert,3,mpert)
 240  real(dp),intent(inout) :: d3etot_1(2,3,mpert,3,mpert,3,mpert)
 241  real(dp),intent(inout) :: d3etot_2(2,3,mpert,3,mpert,3,mpert)
 242  real(dp),intent(inout) :: d3etot_3(2,3,mpert,3,mpert,3,mpert)
 243  real(dp),intent(inout) :: d3etot_4(2,3,mpert,3,mpert,3,mpert)
 244  real(dp),intent(inout) :: d3etot_5(2,3,mpert,3,mpert,3,mpert)
 245  real(dp),intent(inout) :: d3etot_6(2,3,mpert,3,mpert,3,mpert)
 246  real(dp),intent(inout) :: d3etot_7(2,3,mpert,3,mpert,3,mpert)
 247  real(dp),intent(inout) :: d3etot_8(2,3,mpert,3,mpert,3,mpert)
 248  real(dp),intent(inout) :: d3etot_9(2,3,mpert,3,mpert,3,mpert)
 249  type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom*psps%usepaw)
 250  type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
 251  type(pawrhoij_type),intent(in) :: pawrhoij0(natom*psps%usepaw)
 252  type(pawrhoij_type),intent(in),target :: pawrhoij1_i1pert(natom*psps%usepaw)
 253  type(pawrhoij_type),intent(in)        :: pawrhoij1_i2pert(natom*psps%usepaw)
 254  type(pawrhoij_type),intent(in),target :: pawrhoij1_i3pert(natom*psps%usepaw)
 255  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 256  type(paw_an_type),intent(in) :: paw_an0(natom*psps%usepaw)
 257  type(paw_an_type),intent(inout) :: paw_an1_i2pert(natom*psps%usepaw)
 258  type(paw_ij_type),intent(inout) :: paw_ij1_i2pert(natom*psps%usepaw)
 259 
 260 !Local variables-------------------------------
 261 !scalars
 262  logical :: has_cprj_jband,compute_conjugate,compute_rho21
 263  integer,parameter :: level=52,tim_nonlop=0
 264  integer :: bandtot,choice,counter,cplex_cprj,cplex_loc,cplex_rhoij,cpopt,dimffnl1,iband,icg0,ider,ierr,iexit
 265  integer :: idir0,idir_getgh2c,idir_phon,idir_elfd,ipert_phon,ipert_elfd
 266  integer :: ia,iatm,ibg,ii,ikg,ikg1,ikpt,ifft,ifft_re,ifft_im,ilm,isppol,istwf_k,jband
 267  integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nkpg1,nnlout,nsp,nspden_rhoij,npert_phon,npw_k,npw1_k,nzlmopt
 268  integer :: offset_cgi,offset_cgj,offset_eig0,option,paw_opt,debug_mode
 269  integer :: signs,size_wf,size_cprj,spaceComm,typat_ipert_phon,usepaw,useylmgr1
 270  real(dp) :: arg,dot1i,dot1r,dot2i,dot2r,doti,dotr,e3tot,lagi,lagi_paw,lagr,lagr_paw
 271  real(dp) :: rho2ur,rho2ui,rho2dr,rho2di,rho3ur,rho3ui,rho3dr,rho3di
 272  real(dp) :: sumi,sum_psi1H1psi1,sum_psi1H1psi1_i
 273  real(dp) :: sum_lambda1psi1psi1,sum_lambda1psi1psi1_i
 274  real(dp) :: sum_psi0H2psi1a,sum_psi0H2psi1a_i,sum_psi0H2psi1b,sum_psi0H2psi1b_i
 275  real(dp) :: sum_lambda1psi0S1psi1,sum_lambda1psi0S1psi1_i
 276  character(len=1000) :: msg
 277 !arrays
 278  integer,allocatable :: kg_k(:,:),kg1_k(:,:)
 279  real(dp) :: buffer(10),eHxc21_paw(2),eHxc21_nhat(2),exc3(2),exc3_paw(2),kpt(3),eig0_k(mband)
 280  real(dp) :: enlout1(2),enlout2(2)
 281  real(dp) :: rmet(3,3),wtk_k
 282  real(dp),allocatable :: cgi(:,:),cgj(:,:),cg_jband(:,:,:),cwavef1(:,:),cwavef2(:,:),cwavef3(:,:),dkinpw(:)
 283  real(dp),allocatable :: eig1_k_i2pert(:),eig1_k_stored(:)
 284  real(dp),allocatable :: chi_ij(:,:,:,:),cwave_right(:,:),cwave_left(:,:),dudk(:,:),dudkde(:,:),dummy_array(:),dummy_array2(:,:)
 285  real(dp),allocatable :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
 286  real(dp),allocatable :: h_cwave(:,:),iddk(:,:),kinpw1(:),kpg_k(:,:),kpg1_k(:,:),nhat21(:,:),occ_k(:)
 287  real(dp),allocatable :: phkxred(:,:),ph3d(:,:,:),rho1r1_tot(:,:),s_cwave(:,:)
 288  real(dp),allocatable :: vlocal(:,:,:,:),vlocal1_i2pert(:,:,:,:),v_i2pert(:,:),wfraug(:,:,:,:)
 289  real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:)
 290  real(dp),allocatable :: ylm_k(:,:),ylm1_k(:,:),ylmgr1_k(:,:,:)
 291  real(dp),allocatable :: xc_tmp(:,:)
 292  type(pawcprj_type),allocatable :: cwaveprj0(:,:),cwaveprj1(:,:)
 293  type(pawcprj_type),target :: cprj_empty(0,0)
 294  type(pawcprj_type),allocatable,target :: cprj_jband(:,:)
 295  type(pawrhoij_type),allocatable,target  :: pawrhoij21(:)
 296  type(pawrhoij_type),pointer :: pawrhoij21_unsym(:),pawrhoij11(:)
 297  type(paw_ij_type),allocatable :: paw_ij_tmp(:)
 298  type(rf_hamiltonian_type) :: rf_hamkq_i2pert
 299 
 300 !***********************************************************************
 301 
 302  DBG_ENTER("COLL")
 303 
 304  me = mpi_enreg%me
 305  spaceComm=mpi_enreg%comm_cell
 306 
 307  call status(0,dtfil%filstat,iexit,level,'enter         ')
 308 
 309  npert_phon = 0
 310  if(i1pert<=natom) npert_phon = npert_phon + 1
 311  if(i2pert<=natom) npert_phon = npert_phon + 1
 312  if(i3pert<=natom) npert_phon = npert_phon + 1
 313  if (npert_phon>1) then
 314    MSG_ERROR("dfptnl_pert is available with at most one phonon perturbation. Change your input!")
 315  end if
 316 
 317  usepaw = psps%usepaw
 318  size_cprj = nspinor
 319 
 320  call init_rf_hamiltonian(cplex,gs_hamkq,i2pert,rf_hamkq_i2pert,paw_ij1=paw_ij1_i2pert,has_e1kbsc=.true.)
 321 
 322  ABI_ALLOCATE(dummy_array,(0))
 323  ABI_ALLOCATE(dummy_array2,(0,0))
 324 
 325 !Acivate computation of rho^(2:1) and related energy derivatives if needed
 326  compute_rho21 = .false.
 327  if (usepaw==1.and.npert_phon==1.and.(i1pert<=natom.or.i3pert<=natom)) then ! so i2pert==natom+2
 328    compute_rho21 = .true.
 329    if (i1pert<=natom) then
 330      ipert_phon  = i1pert
 331      idir_phon   = i1dir
 332      ipert_elfd = i3pert
 333      idir_elfd  = i3dir
 334      pawrhoij11 => pawrhoij1_i3pert
 335    else if (i3pert<=natom) then
 336      ipert_phon  = i3pert
 337      idir_phon   = i3dir
 338      ipert_elfd = i1pert
 339      idir_elfd  = i1dir
 340      pawrhoij11 => pawrhoij1_i1pert
 341    end if
 342    cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=dtset%nspden
 343    ABI_DATATYPE_ALLOCATE(pawrhoij21,(natom))
 344    call pawrhoij_nullify(pawrhoij21)
 345    call pawrhoij_alloc(pawrhoij21,cplex_rhoij,nspden_rhoij,nspinor,dtset%nsppol,&
 346 &   dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 347    ABI_DATATYPE_ALLOCATE(cwaveprj0,(natom,size_cprj))
 348    ABI_DATATYPE_ALLOCATE(cwaveprj1,(natom,size_cprj))
 349    call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj)
 350    call pawcprj_alloc(cwaveprj1,1,gs_hamkq%dimcprj)
 351 !   if (paral_atom) then
 352 !     ABI_DATATYPE_ALLOCATE(pawrhoij1_unsym,(natom))
 353 !     cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=dtset%nspden
 354 !     call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
 355 !&     dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
 356 !   else
 357    pawrhoij21_unsym => pawrhoij21
 358    call pawrhoij_init_unpacked(pawrhoij21_unsym)
 359 !   end if
 360 !  Compute phkxred :
 361    ABI_ALLOCATE(phkxred,(2,natom))
 362    do ia=1,natom
 363      iatm=min(atindx(ia),natom)
 364      arg=two_pi*(kpt(1)*xred(1,ia)+kpt(2)*xred(2,ia)+kpt(3)*xred(3,ia))
 365      phkxred(1,iatm)=cos(arg);phkxred(2,iatm)=sin(arg)
 366    end do
 367    ABI_DATATYPE_ALLOCATE(paw_ij_tmp,(natom))
 368    call paw_ij_nullify(paw_ij_tmp)
 369    cplex_loc=1;nsp=1 ! Force nsppol/nspden to 1 because Dij^(1) due to electric field is spin-independent
 370    call paw_ij_init(paw_ij_tmp,cplex_loc,dtset%nspinor,nsp,nsp,dtset%pawspnorb,natom,psps%ntypat,&
 371 &   dtset%typat,pawtab,has_dijfr=1,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 372    call paw_ij_reset_flags(paw_ij_tmp,all=.True.)
 373    call pawdijfr(cplex_loc,gs_hamkq%gprimd,idir_elfd,ipert_elfd,natom,natom,nfftf,ngfftf,nsp,nsp,psps%ntypat,&
 374 &   1,paw_ij_tmp,pawang,pawfgrtab,pawrad,pawtab,&
 375 &   (/zero,zero,zero/),rprimd,ucvol,dummy_array2,dummy_array2,dummy_array2,xred,&
 376 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 377    ABI_ALLOCATE(chi_ij,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,dtset%nspinor**2,cplex_loc))
 378    call pawdij2e1kb(paw_ij_tmp,1,mpi_enreg%comm_atom,mpi_enreg%my_atmtab,e1kbfr=chi_ij)
 379    call paw_ij_free(paw_ij_tmp)
 380    ABI_DATATYPE_DEALLOCATE(paw_ij_tmp)
 381  else
 382    ABI_ALLOCATE(chi_ij,(0,0,0,0))
 383    ABI_ALLOCATE(phkxred,(0,0))
 384    ABI_DATATYPE_ALLOCATE(pawrhoij21,(0))
 385    pawrhoij21_unsym => pawrhoij21
 386    ABI_DATATYPE_ALLOCATE(cwaveprj0,(0,0))
 387    ABI_DATATYPE_ALLOCATE(cwaveprj1,(0,0))
 388  end if
 389 
 390  nnlout = 0
 391  bandtot = 0
 392  icg0 = 0
 393  option = 2
 394  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
 395  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
 396 
 397  ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
 398  ABI_ALLOCATE(vlocal1_i2pert,(cplex*n4,n5,n6,gs_hamkq%nvloc))
 399 
 400  ABI_ALLOCATE(wfraug,(2,n4,n5,n6))
 401 
 402  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
 403 
 404  sumi = zero
 405 
 406 !Set up the Ylm for each k point
 407  ABI_ALLOCATE(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
 408  ABI_ALLOCATE(ylmgr,(dtset%mpw*dtset%mkmem,9,psps%mpsang*psps%mpsang*psps%useylm))
 409  if (psps%useylm==1) then
 410    call status(0,dtfil%filstat,iexit,level,'call initylmg ')
 411    option=2
 412    call initylmg(gs_hamkq%gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,&
 413    dtset%nkpt,npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
 414  end if
 415 
 416 !Set up the spherical harmonics (Ylm) at k+q
 417  useylmgr1=0; option=0
 418  if (psps%useylm==1.and. &
 419 & (i2pert==natom+1.or.i2pert==natom+3.or.i2pert==natom+4.or.(usepaw==1.and.i2pert==natom+2))) then
 420    useylmgr1=1; option=1
 421  end if
 422  ABI_ALLOCATE(ylm1,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
 423  ABI_ALLOCATE(ylmgr1,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
 424 !To change the following when q/=0
 425  if (psps%useylm==1) then
 426    call initylmg(gs_hamkq%gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,&
 427    dtset%nkpt,npwarr,dtset%nsppol,option,rprimd,ylm1,ylmgr1)
 428  end if
 429 
 430  debug_mode = 0
 431  if (dtset%nonlinear_info>3.or.dtset%nonlinear_info==2) debug_mode = 1
 432 
 433 !Real parts
 434  sum_psi1H1psi1 =  zero
 435  sum_lambda1psi1psi1 = zero
 436  sum_lambda1psi0S1psi1 = zero
 437  sum_psi0H2psi1a = zero
 438  sum_psi0H2psi1b = zero
 439 !Imaginary parts
 440  sum_psi1H1psi1_i =  zero
 441  sum_lambda1psi1psi1_i = zero
 442  sum_lambda1psi0S1psi1_i = zero
 443  sum_psi0H2psi1a_i = zero
 444  sum_psi0H2psi1b_i = zero
 445 
 446  compute_conjugate = .false.
 447 !We have to compute < u^(ip1) | H^(ip2) | u^(ip3) >
 448 !For some cases, we want to apply H^(ip2) on < u^(ip1) |, not on | u^(ip3) > (see below)
 449  if (i3pert<=natom) then ! As npert_phon<=1, this implies that i1pert=natom+2 and i2pert=natom+2
 450    compute_conjugate = .true.
 451  end if
 452 
 453 !Loop over spins
 454  do isppol = 1, nsppol
 455 
 456 !  Set up local potential vlocal1 with proper dimensioning, from vtrial1
 457 !  Same thing for vlocal from vtrial Also take into account the spin.
 458    call rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,&
 459 &   gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1_i2pert,vlocal,vlocal1_i2pert)
 460 
 461 !  Continue to initialize the Hamiltonian
 462    call load_spin_hamiltonian(gs_hamkq,isppol,vlocal=vlocal,with_nonlocal=.true.)
 463    call load_spin_rf_hamiltonian(rf_hamkq_i2pert,isppol,vlocal1=vlocal1_i2pert,with_nonlocal=.true.)
 464 
 465 !  Loop over k-points
 466 
 467    ikg = 0
 468    ikg1 = 0
 469 
 470    do ikpt = 1, nkpt
 471 
 472      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,mpi_enreg%me)) then
 473        cycle ! Skip the rest of the k-point loop
 474      end if
 475 
 476      counter = 100*ikpt
 477 
 478      nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
 479      npw_k = npwarr(ikpt)
 480      npw1_k = npw_k ! To change for q/=0
 481      istwf_k = dtset%istwfk(ikpt)
 482      ABI_ALLOCATE(occ_k,(nband_k))
 483      occ_k(:) = occ(1+bandtot:nband_k+bandtot)
 484      wtk_k    = dtset%wtk(ikpt)
 485 
 486      size_wf = nspinor*npw_k
 487 
 488      kpt(:) = dtset%kptns(:,ikpt)
 489 
 490      ABI_ALLOCATE(cwavef1,(2,npw_k*nspinor))
 491      ABI_ALLOCATE(cwavef3,(2,npw_k*nspinor))
 492      if (compute_rho21) then
 493        ABI_ALLOCATE(cwavef2,(2,npw_k*nspinor))
 494      end if
 495 
 496      ABI_ALLOCATE(kg_k,(3,npw_k))
 497      ABI_ALLOCATE(kg1_k,(3,npw1_k))
 498      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
 499      ABI_ALLOCATE(ylm1_k,(npw1_k,mpsang*mpsang*psps%useylm))
 500      ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
 501 
 502 !    Get (k+G) wave vectors and associated spherical harmonics
 503      kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg)
 504      if (psps%useylm==1) then
 505        do ilm=1,mpsang*mpsang
 506          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
 507        end do
 508      end if
 509 
 510 !    Get (k+q+G) wave vectors and associated spherical harmonics
 511      kg1_k(:,1:npw1_k)=kg(:,1+ikg1:npw1_k+ikg1) ! To change for q/=0
 512      if (psps%useylm==1) then
 513        do ilm=1,psps%mpsang*psps%mpsang
 514          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
 515        end do
 516        if (useylmgr1==1) then
 517          do ilm=1,psps%mpsang*psps%mpsang
 518            do ii=1,3
 519              ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
 520            end do
 521          end do
 522        end if
 523      end if
 524 
 525 !    Compute (k+G) vectors
 526      nkpg=0;if(i2pert>=1.and.i2pert<=natom) nkpg=3*dtset%nloalg(3)
 527      ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 528      if (nkpg>0) then
 529        call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
 530      end if
 531 
 532 !    Compute (k+q+G) vectors
 533      nkpg1=0;if(i2pert>=1.and.i2pert<=natom) nkpg1=3*dtset%nloalg(3)
 534      ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1))
 535      if (nkpg1>0) then
 536        call mkkpg(kg1_k,kpg1_k,kpt,nkpg1,npw1_k)
 537      end if
 538 
 539 !    ===== Preparation of the non-local contributions
 540 
 541 !    Compute nonlocal form factors ffnl1 at (k+q+G)
 542      !-- Atomic displacement perturbation
 543      if (i2pert<=natom) then
 544        ider=0;idir0=0
 545      !-- Electric field perturbation
 546      else if (i2pert==natom+2) then
 547        if (psps%usepaw==1) then
 548          ider=1;idir0=i2dir
 549        else
 550          ider=0;idir0=0
 551        end if
 552      end if
 553      if (compute_rho21) then ! compute_rho21 implies i2pert==natom+2
 554        ider=1; idir0=4
 555      end if
 556 
 557 !    Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
 558      dimffnl1=1+ider
 559      if (ider==1.and.(idir0==0.or.idir0==4)) dimffnl1=2+2*psps%useylm
 560 !     if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm
 561      ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
 562      call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,&
 563 &     psps%indlmn,kg1_k,kpg1_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
 564 &     npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
 565 
 566 !    Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
 567      if (compute_rho21.and.debug_mode/=0) then
 568        ABI_ALLOCATE(ffnl1_test,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
 569        idir0 = 0 ! for nonlop with signs = 1
 570        call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1_test,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,&
 571 &       psps%indlmn,kg1_k,kpg1_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
 572 &       npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
 573      end if
 574 
 575 !    ===== Preparation of the kinetic contributions
 576 
 577 !    Note that not all these arrays should be allocated in the general case when wtk_k vanishes
 578 
 579 !    Compute (1/2) (2 Pi)**2 (k+q+G)**2:
 580      ABI_ALLOCATE(kinpw1,(npw1_k))
 581      kinpw1(:)=zero
 582      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gs_hamkq%gmet,kg1_k,kinpw1,kpt,npw1_k,0,0)
 583 
 584      ABI_ALLOCATE(dkinpw,(npw_k)) ! 1st derivative (1st direction)
 585      dkinpw(:)=zero
 586 
 587 !===== Load the k/k+q dependent parts of the Hamiltonian
 588 
 589 !  Load k-dependent part in the Hamiltonian datastructure
 590      ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamkq%matblk))
 591      call load_k_hamiltonian(gs_hamkq,kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,&
 592 &     ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
 593      call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1,kpt_k=kpt,npw_k=npw1_k,istwf_k=istwf_k,&
 594 &     kinpw_k=kinpw1,kg_k=kg1_k,kpg_k=kpg1_k,compute_gbound=.true.)
 595 !   end if
 596 
 597 !    Load k-dependent part in the 1st-order Hamiltonian datastructure
 598      call load_k_rf_hamiltonian(rf_hamkq_i2pert,npw_k=npw_k,dkinpw_k=dkinpw)
 599 
 600      ABI_STAT_ALLOCATE(dudk,  (2,nband_k*size_wf), ierr)
 601      ABI_STAT_ALLOCATE(dudkde,(2,nband_k*size_wf), ierr)
 602      ABI_STAT_ALLOCATE(eig1_k_i2pert,(2*nband_k), ierr)
 603      ABI_ALLOCATE(eig1_k_stored,(2*nband_k**2))
 604      ABI_ALLOCATE(cgi,(2,size_wf))
 605      ABI_ALLOCATE(cwave_right,(2,size_wf))
 606      ABI_ALLOCATE(cwave_left,(2,size_wf))
 607 
 608 ! **************************************************************************************************
 609 !      Read dudk and dudkde
 610 ! **************************************************************************************************
 611 
 612      do iband = 1,nband_k
 613        if (occ_k(iband)>tol10) then
 614 
 615 !        Read dude file
 616          call wfk_read_bks(ddk_f(1), iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert)
 617 !        Copy eig1_k_i2pert in "eig1_k_stored"
 618          eig1_k_stored(1+(iband-1)*2*nband_k:2*nband_k+(iband-1)*2*nband_k)=eig1_k_i2pert(:)
 619 
 620          if (i2pert==natom+2) then
 621 !          Read dudk file
 622            call wfk_read_bks(ddk_f(2), iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert)
 623            offset_cgi = (iband-1)*size_wf+icg0
 624            cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi)
 625 !          Copy cwave_right in "dudk"
 626            dudk(:,1+(iband-1)*size_wf:iband*size_wf)=cwave_right(:,:)
 627 
 628 !          Read dudkde file
 629            call wfk_read_bks(ddk_f(3), iband, ikpt, isppol, xmpio_single, cg_bks=cwave_right,eig1_bks=eig1_k_i2pert)
 630            offset_cgi = (iband-1)*size_wf+icg0
 631            cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi)
 632 !          Copy cwave_right in "dudkde"
 633            dudkde(:,1+(iband-1)*size_wf:iband*size_wf)=cwave_right(:,:)
 634          end if
 635 
 636        end if
 637      end do
 638 
 639      ABI_ALLOCATE(cgj,(2,size_wf))
 640      ABI_ALLOCATE(iddk,(2,size_wf))
 641 
 642      offset_eig0 = mband*(ikpt-1+nkpt*(isppol-1))
 643      eig0_k(:) = eigen0(1+offset_eig0:mband+offset_eig0)
 644 
 645      ABI_STAT_ALLOCATE(h_cwave,(2,size_wf), ierr)
 646      ABI_STAT_ALLOCATE(s_cwave,(2,size_wf), ierr)
 647 
 648 !    Allocate work spaces when debug_mode is activated
 649      has_cprj_jband=.false.
 650      if (debug_mode/=0) then ! Only for test purposes
 651        ABI_ALLOCATE(cg_jband,(2,size_wf*nband_k,2))
 652        cg_jband(:,:,1) = cg(:,1+icg0:size_wf*nband_k+icg0)
 653        if (i2pert==natom+2) then ! Note the multiplication by "i"
 654          cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k)
 655          cg_jband(2,:,2) =  dudk(1,1:size_wf*nband_k)
 656        end if
 657        if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
 658          ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,size_cprj*nband_k))
 659          has_cprj_jband=.true.
 660        else
 661          ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,0))
 662        end if
 663      else
 664        ABI_ALLOCATE(cg_jband,(2,0,2))
 665        ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,0))
 666      end if
 667 
 668 !    Loop over bands
 669      do jband = 1,nband_k
 670        !  Skip bands not treated by current proc
 671        if((mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me)) cycle
 672 
 673        if (occ_k(jband)>tol10) then
 674 
 675   !      tol_test = tol8
 676          offset_cgj = (jband-1)*size_wf+icg0
 677          cgj(:,:) = cg(:,1+offset_cgj:size_wf+offset_cgj)
 678 
 679 ! **************************************************************************************************
 680 !        Compute < u^(1) | ( H^(1) - eps^(0) S^(1) ) | u^(1) >
 681 ! **************************************************************************************************
 682 
 683          eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k)
 684          cwavef1(:,:) = cg1(:,1+offset_cgj:size_wf+offset_cgj)
 685          cwavef3(:,:) = cg3(:,1+offset_cgj:size_wf+offset_cgj)
 686          if (i2pert==natom+2) then ! Note the multiplication by i
 687            iddk(1,:) = -dudkde(2,1+(jband-1)*size_wf:jband*size_wf)
 688            iddk(2,:) =  dudkde(1,1+(jband-1)*size_wf:jband*size_wf)
 689          else
 690            iddk(:,:) = zero
 691          end if
 692          cwave_right(:,:) = cwavef3(:,:)
 693          cwave_left(:,:)  = cwavef1(:,:)
 694          if (compute_conjugate) then
 695            cwave_right(:,:) = cwavef1(:,:)
 696            cwave_left(:,:)  = cwavef3(:,:)
 697          end if
 698 
 699   !      Compute : < u^(ip1) | ( H^(ip2) - eps^(0) S^(ip2) ) | u^(ip3) >
 700   !           or : < u^(ip3) | ( H^(ip2) - eps^(0) S^(ip2) ) | u^(ip1) >
 701          call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_right,cprj_empty,h_cwave,s_cwave,eig0_k,eig1_k_i2pert,&
 702 &         jband,gs_hamkq,iddk,i2dir,i2pert,ikpt,isppol,mkmem,mpi_enreg,nband_k,nsppol,&
 703 &         debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf)
 704          call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_left,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 705 
 706          if (usepaw==1.and.i2pert/=natom+2) then ! S^(1) is zero for ipert=natom+2
 707            call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwave_left,s_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 708            dotr = dotr - eig0_k(jband)*dot2r
 709            doti = doti - eig0_k(jband)*dot2i
 710          end if
 711          if (compute_conjugate) doti = -doti
 712 
 713 ! **************************************************************************************************
 714 !        Compute sum_i Lambda_ij^(1) < u_j^(1) | u_i^(1)>
 715 ! **************************************************************************************************
 716 
 717 !         eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k)
 718          lagr = zero ; lagi = zero
 719          lagr_paw = zero ; lagi_paw = zero
 720 
 721          do iband = 1, nband_k
 722            if(occ_k(iband)>tol10) then
 723 
 724              offset_cgi = (iband-1)*size_wf+icg0
 725              cwavef3(:,:) = cg3(:,1+offset_cgi:size_wf+offset_cgi)
 726 
 727              if(debug_mode/=0) then
 728                eig1_k_i2pert(:) = eig1_k_stored(1+(jband-1)*2*nband_k:jband*2*nband_k)
 729              end if
 730 
 731 !            Get Lambda_ij^(1) = < u_i^(0) | H^(1) - eps^(0) S^(1) | u_j^(0) > (see dfpt_cgwf.F90)
 732              dot1r = eig1_k_i2pert(2*iband-1)
 733              dot1i = eig1_k_i2pert(2*iband  )
 734 
 735 !            Compute < u_j^(1) | S^(0) | u_i^(1) >
 736              if (usepaw==1) then
 737                ibg = 0
 738                call getgsc(cwavef3,cprj_empty,gs_hamkq,s_cwave,ibg,0,0,ikpt,isppol,&
 739 &               size_wf,size_cprj,size_wf,mpi_enreg,natom,-1,npw_k,nspinor,select_k=KPRIME_H_KPRIME)
 740              else
 741                s_cwave(:,:) = cwavef3(:,:)
 742              end if
 743 
 744              call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwavef1,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
 745              lagr = lagr + dot1r*dot2r - dot1i*dot2i
 746              lagi = lagi + dot1r*dot2i + dot1i*dot2r
 747 
 748 !            Compute < u_j^(0) | S^(1) | u_i^(1) >
 749              if (usepaw==1.and.i1pert<=natom) then ! S^(1) is zero for ipert=natom+2
 750 
 751                cpopt=-1+5*gs_hamkq%usecprj ; choice=2 ; signs=2 ; paw_opt=3
 752                call nonlop(choice,cpopt,cprj_empty,dummy_array,gs_hamkq,i1dir,(/zero/),mpi_enreg,1,nnlout,&
 753 &               paw_opt,signs,s_cwave,tim_nonlop,cwavef3,dummy_array2,iatom_only=i1pert)
 754                call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cgj,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
 755                lagr_paw = lagr_paw + dot1r*dot2r - dot1i*dot2i
 756                lagi_paw = lagi_paw + dot1r*dot2i + dot1i*dot2r
 757 
 758              end if
 759 
 760 !            Compute < u_j^(1) | S^(1) | u_i^(0) >
 761              if (usepaw==1.and.i3pert<=natom) then ! S^(1) is zero for ipert=natom+2
 762 
 763                cgi(:,:) = cg(:,1+offset_cgi:size_wf+offset_cgi)
 764                cpopt=-1+5*gs_hamkq%usecprj ; choice=2 ; signs=2 ; paw_opt=3
 765                call nonlop(choice,cpopt,cprj_empty,dummy_array,gs_hamkq,i3dir,(/zero/),mpi_enreg,1,nnlout,&
 766 &               paw_opt,signs,s_cwave,tim_nonlop,cgi,dummy_array2,iatom_only=i3pert)
 767                call dotprod_g(dot2r,dot2i,gs_hamkq%istwf_k,size_wf,2,cwavef1,s_cwave,mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
 768                lagr_paw = lagr_paw + dot1r*dot2r - dot1i*dot2i
 769                lagi_paw = lagi_paw + dot1r*dot2i + dot1i*dot2r
 770 
 771              end if
 772 
 773            end if
 774          end do    ! iband
 775 
 776 ! **************************************************************************************************
 777 !        Sum all band_by_band contributions
 778 ! **************************************************************************************************
 779 
 780 !        Real part
 781          sum_psi1H1psi1 = sum_psi1H1psi1 + dtset%wtk(ikpt)*occ_k(jband)*dotr
 782          sum_lambda1psi1psi1 = sum_lambda1psi1psi1 - dtset%wtk(ikpt)*occ_k(jband)*lagr
 783          sum_lambda1psi0S1psi1 = sum_lambda1psi0S1psi1 - dtset%wtk(ikpt)*occ_k(jband)*lagr_paw
 784 
 785 !        Imaginary part
 786          sum_psi1H1psi1_i = sum_psi1H1psi1_i + dtset%wtk(ikpt)*occ_k(jband)*doti
 787          sum_lambda1psi1psi1_i = sum_lambda1psi1psi1_i - dtset%wtk(ikpt)*occ_k(jband)*lagi
 788          sum_lambda1psi0S1psi1_i = sum_lambda1psi0S1psi1_i - dtset%wtk(ikpt)*occ_k(jband)*lagi_paw
 789 
 790 ! **************************************************************************************************
 791 !        If compute_rho21 : accumulate rhoij and compute term with H_KV^(2)
 792 ! **************************************************************************************************
 793 
 794          if (compute_rho21) then
 795 
 796            if (i1pert<=natom) then ! If true, i3pert==natom+2
 797              cwave_right = cg3(:,1+offset_cgj:size_wf+offset_cgj)
 798            else if (i3pert<=natom) then  ! If true, i1pert==natom+2
 799              cwave_right = cg1(:,1+offset_cgj:size_wf+offset_cgj)
 800            end if
 801            choice = 2
 802            cpopt  = 0
 803            call getcprj(choice,cpopt,cgj,cwaveprj0,&
 804 &           ffnl1,idir_phon,psps%indlmn,gs_hamkq%istwf_k,kg_k,kpg_k,kpt,psps%lmnmax,&
 805 &           mgfft,mpi_enreg,natom,nattyp,dtset%ngfft,dtset%nloalg,&
 806 &           npw_k,nspinor,psps%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm)
 807            call getcprj(choice,cpopt,cwave_right,cwaveprj1,&
 808 &           ffnl1,idir_phon,psps%indlmn,gs_hamkq%istwf_k,kg_k,kpg_k,kpt,psps%lmnmax,&
 809 &           mgfft,mpi_enreg,natom,nattyp,dtset%ngfft,dtset%nloalg,&
 810 &           npw_k,nspinor,psps%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm)
 811 
 812            cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1
 813            call paw_dfptnl_accrhoij(atindx,cplex_cprj,cwaveprj0,cwaveprj0,cwaveprj1,cwaveprj1,i1pert,i3pert,isppol,natom,natom,&
 814 &           nspinor,occ_k(jband),pawrhoij21_unsym,wtk_k)
 815 !&           comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
 816 
 817 !          Compute < psi^(0) | H_KV^(pert1pert3) | psi^(pert2) > + < psi^(pert2) | H_KV^(pert1pert3) | psi^(0) >
 818            cwavef2(:,:) = cg2(:,1+offset_cgj:size_wf+offset_cgj)
 819 
 820 !          Read dkk file (for tests only)
 821            if (debug_mode/=0) then
 822              if(idir_elfd==i2dir) then
 823                call wfk_read_bks(ddk_f(2), jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
 824              else
 825                call wfk_read_bks(ddk_f(4), jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
 826              end if
 827              cg_jband(1,1+size_wf*(jband-1):size_wf*jband,2) = -iddk(2,1:size_wf)
 828              cg_jband(2,1+size_wf*(jband-1):size_wf*jband,2) =  iddk(1,1:size_wf)
 829            end if
 830 !          Read dkde file
 831            if(idir_elfd==i2dir) then
 832              call wfk_read_bks(ddk_f(3), jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
 833            else
 834              call wfk_read_bks(ddk_f(5), jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
 835            end if
 836            s_cwave = iddk
 837            iddk(1,:) = -s_cwave(2,:)
 838            iddk(2,:) =  s_cwave(1,:)
 839            call rf2_getidir(idir_phon,idir_elfd,idir_getgh2c)
 840            call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwavef2,cprj_empty,s_cwave,dummy_array2,eig0_k,eig1_k_i2pert,&
 841 &           jband,gs_hamkq,iddk,idir_getgh2c,ipert_phon+natom+11,ikpt,isppol,mkmem,mpi_enreg,nband_k,nsppol,&
 842 &           debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf,enl=chi_ij,ffnl1=ffnl1,ffnl1_test=ffnl1_test)
 843            call dotprod_g(enlout1(1),enlout1(2),gs_hamkq%istwf_k,npw_k*nspinor,2,cgj,s_cwave,&
 844 &           mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 845            sum_psi0H2psi1a   = sum_psi0H2psi1a   + half*dtset%wtk(ikpt)*occ_k(jband)*enlout1(1) ! be careful : factor 0.5
 846            sum_psi0H2psi1a_i = sum_psi0H2psi1a_i + half*dtset%wtk(ikpt)*occ_k(jband)*enlout1(2) ! be careful : factor 0.5
 847 
 848 !          Read ddk file
 849            if(idir_elfd==i2dir) then
 850              call wfk_read_bks(ddk_f(2), jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
 851            else
 852              call wfk_read_bks(ddk_f(4), jband, ikpt, isppol, xmpio_single, cg_bks=iddk)
 853            end if
 854            s_cwave = iddk
 855            iddk(1,:) = -s_cwave(2,:)
 856            iddk(2,:) =  s_cwave(1,:)
 857            call rf2_apply_hamiltonian(cg_jband,cprj_jband,cgj,cprj_empty,s_cwave,dummy_array2,eig0_k,eig1_k_i2pert,&
 858 &           jband,gs_hamkq,iddk,idir_getgh2c,ipert_phon+natom+11,ikpt,isppol,mkmem,mpi_enreg,nband_k,nsppol,&
 859 &           debug_mode,dtset%prtvol,rf_hamkq_i2pert,size_cprj,size_wf,enl=chi_ij,ffnl1=ffnl1,ffnl1_test=ffnl1_test)
 860            call dotprod_g(enlout2(1),enlout2(2),gs_hamkq%istwf_k,npw_k*nspinor,2,cwavef2,s_cwave,&
 861 &           mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 862            sum_psi0H2psi1b   = sum_psi0H2psi1b   + half*dtset%wtk(ikpt)*occ_k(jband)*enlout2(1) ! be careful : factor 0.5
 863            sum_psi0H2psi1b_i = sum_psi0H2psi1b_i + half*dtset%wtk(ikpt)*occ_k(jband)*enlout2(2) ! be careful : factor 0.5
 864 
 865          end if ! end if compute_rho21
 866 
 867        end if
 868      end do   ! end loop over jband
 869 
 870 ! **************************************************************************************************
 871 !    END OF BAND LOOP
 872 ! **************************************************************************************************
 873 
 874      ABI_DEALLOCATE(cgi)
 875      ABI_DEALLOCATE(cgj)
 876      ABI_DEALLOCATE(iddk)
 877      ABI_DEALLOCATE(h_cwave)
 878      ABI_DEALLOCATE(s_cwave)
 879 
 880      ABI_DEALLOCATE(cwave_right)
 881      ABI_DEALLOCATE(cwave_left)
 882      ABI_DEALLOCATE(eig1_k_i2pert)
 883      ABI_DEALLOCATE(eig1_k_stored)
 884 
 885      bandtot = bandtot + nband_k
 886      icg0 = icg0 + npw_k*nspinor*nband_k
 887      ikg = ikg + npw_k
 888      ikg1 = ikg1 + npw1_k
 889 
 890      ABI_DEALLOCATE(cwavef1)
 891      if (compute_rho21) then
 892        ABI_DEALLOCATE(cwavef2)
 893      end if
 894      ABI_DEALLOCATE(cwavef3)
 895      ABI_DEALLOCATE(dkinpw)
 896      ABI_DEALLOCATE(kg_k)
 897      ABI_DEALLOCATE(kg1_k)
 898      ABI_DEALLOCATE(kinpw1)
 899      ABI_DEALLOCATE(dudk)
 900      ABI_DEALLOCATE(dudkde)
 901      ABI_DEALLOCATE(ylm_k)
 902      ABI_DEALLOCATE(ylm1_k)
 903      ABI_DEALLOCATE(ylmgr1_k)
 904      ABI_DEALLOCATE(ffnl1)
 905      if (compute_rho21.and.debug_mode/=0) then
 906        ABI_DEALLOCATE(ffnl1_test)
 907      end if
 908      ABI_DEALLOCATE(kpg_k)
 909      ABI_DEALLOCATE(kpg1_k)
 910      ABI_DEALLOCATE(cg_jband)
 911      ABI_DEALLOCATE(occ_k)
 912      ABI_DEALLOCATE(ph3d)
 913      if (has_cprj_jband) call pawcprj_free(cprj_jband)
 914      ABI_DATATYPE_DEALLOCATE(cprj_jband)
 915 
 916    end do   ! end loop over k-points
 917 
 918  end do   ! end loop over spins
 919 
 920  call destroy_rf_hamiltonian(rf_hamkq_i2pert)
 921 
 922 ! **************************************************************************************************
 923 !    GATHER BAND-BY-BAND AND XC CONTRIBUTIONS
 924 ! **************************************************************************************************
 925 
 926  if (xmpi_paral == 1) then
 927 
 928 !  Real parts
 929    buffer(1)  = sum_psi1H1psi1
 930    buffer(2)  = sum_lambda1psi1psi1
 931    buffer(3)  = sum_lambda1psi0S1psi1
 932    buffer(4)  = sum_psi0H2psi1a
 933    buffer(5)  = sum_psi0H2psi1b
 934 
 935 !  Imaginary parts
 936    buffer(6)  = sum_psi1H1psi1_i
 937    buffer(7)  = sum_lambda1psi1psi1_i
 938    buffer(8)  = sum_lambda1psi0S1psi1_i
 939    buffer(9)  = sum_psi0H2psi1a_i
 940    buffer(10) = sum_psi0H2psi1b_i
 941 
 942    call xmpi_sum(buffer,spaceComm,ierr)
 943 
 944 !  Real parts
 945    sum_psi1H1psi1          = buffer(1)
 946    sum_lambda1psi1psi1     = buffer(2)
 947    sum_lambda1psi0S1psi1   = buffer(3)
 948    sum_psi0H2psi1a         = buffer(4)
 949    sum_psi0H2psi1b         = buffer(5)
 950 
 951 !  Imaginary parts
 952    sum_psi1H1psi1_i        = buffer(6)
 953    sum_lambda1psi1psi1_i   = buffer(7)
 954    sum_lambda1psi0S1psi1_i = buffer(8)
 955    sum_psi0H2psi1a_i       = buffer(9)
 956    sum_psi0H2psi1b_i       = buffer(10)
 957 
 958 !  Accumulate PAW occupancies
 959    if (compute_rho21) then
 960      call pawrhoij_mpisum_unpacked(pawrhoij21_unsym,spaceComm)
 961    end if
 962 
 963  end if
 964 
 965 ! **************************************************************************************************
 966 !      Compute E_xc^(3) (NOTE : E_H^(3) = 0)
 967 ! **************************************************************************************************
 968 ! This term is equal to the term "exc3" in pead_nl_loop.F90. However, we treat the nonlinear xc core correction in
 969 ! a slightly different way, in order to keep the symmetry between pert1,pert2 and pert3. It helps for debugging.
 970 
 971  ABI_ALLOCATE(xc_tmp,(cplex*nfftf,nspden))
 972  ABI_ALLOCATE(rho1r1_tot,(cplex*nfftf,nspden))
 973  if (nspden==1)then
 974 
 975    if (cplex==1) then
 976      do ifft=1,nfftf
 977        rho1r1_tot(ifft,1) = rho1r1(ifft,1) + xccc3d1(ifft)
 978        rho2ur = rho2r1(ifft,1)+xccc3d2(ifft)
 979        rho3ur = rho3r1(ifft,1)+xccc3d3(ifft)
 980        xc_tmp(ifft,1)= k3xc(ifft,1)*rho2ur*rho3ur
 981      end do
 982    else
 983      do ifft=1,nfftf
 984        ifft_re = 2*ifft-1 ! Real part
 985        ifft_im = 2*ifft   ! Imaginary part
 986        rho1r1_tot(ifft_re,1) = rho1r1(ifft_re,1) + xccc3d1(ifft_re)
 987        rho1r1_tot(ifft_im,1) = rho1r1(ifft_im,1) + xccc3d1(ifft_im)
 988        rho2ur = rho2r1(ifft_re,1)+xccc3d2(ifft_re)
 989        rho2ui = rho2r1(ifft_im,1)+xccc3d2(ifft_im)
 990        rho3ur = rho3r1(ifft_re,1)+xccc3d3(ifft_re)
 991        rho3ui = rho3r1(ifft_im,1)+xccc3d3(ifft_im)
 992        xc_tmp(ifft_re,1)= k3xc(ifft,1)*(rho2ur*rho3ur-rho2ui*rho3ui)
 993        xc_tmp(ifft_im,1)= k3xc(ifft,1)*(rho2ur*rho3ui+rho2ui*rho3ur)
 994      end do
 995    end if
 996 
 997  else if (nspden==2) then
 998 
 999 !  Remember : rhor(...,1) = total density  ( n_up(r) + n_down(r) )
1000 !             rhor(...,2) =    up density  ( n_up(r) )
1001 !       But :  pot(...,1) =   up potential ( v_up(r) )
1002 !              pot(...,2) = down potential ( v_down(r) )
1003 
1004    if (cplex==1) then
1005      do ifft=1,nfftf
1006        rho1r1_tot(ifft,1) = rho1r1(ifft,1)    + xccc3d1(ifft)      ! 1 tot
1007        rho1r1_tot(ifft,2) = rho1r1(ifft,2)    + xccc3d1(ifft)*half ! 1 up
1008        rho2ur = rho2r1(ifft,2)                + xccc3d2(ifft)*half ! 2 up
1009        rho2dr = rho2r1(ifft,1)-rho2r1(ifft,2) + xccc3d2(ifft)*half ! 2 down
1010        rho3ur = rho3r1(ifft,2)                + xccc3d3(ifft)*half ! 3 up
1011        rho3dr = rho3r1(ifft,1)-rho3r1(ifft,2) + xccc3d3(ifft)*half ! 3 down
1012 !                      uuu                          uud
1013        xc_tmp(ifft,1)= k3xc(ifft,1)*rho2ur*rho3ur + k3xc(ifft,2)*rho2ur*rho3dr + &
1014 !                      udu                          udd
1015 &       k3xc(ifft,2)*rho2dr*rho3ur + k3xc(ifft,3)*rho2dr*rho3dr
1016 !                      duu                          dud
1017        xc_tmp(ifft,2)= k3xc(ifft,2)*rho2ur*rho3ur + k3xc(ifft,3)*rho2ur*rho3dr + &
1018 !                      ddu                          ddd
1019 &       k3xc(ifft,3)*rho2dr*rho3ur + k3xc(ifft,4)*rho2dr*rho3dr
1020      end do
1021 
1022    else ! cplex = 2
1023 
1024      do ifft=1,nfftf
1025        ifft_re = 2*ifft-1 ! Real part
1026        ifft_im = 2*ifft   ! Imaginary part
1027        rho1r1_tot(ifft_re,1) = rho1r1(ifft_re,1)    + xccc3d1(ifft_re)      ! 1 tot re
1028        rho1r1_tot(ifft_im,1) = rho1r1(ifft_im,1)    + xccc3d1(ifft_im)      ! 1 tot im
1029        rho1r1_tot(ifft_re,2) = rho1r1(ifft_re,2)    + xccc3d1(ifft_re)*half ! 1 up re
1030        rho1r1_tot(ifft_im,2) = rho1r1(ifft_im,2)    + xccc3d1(ifft_im)*half ! 1 up im
1031        rho2ur = rho2r1(ifft_re,2)                   + xccc3d2(ifft_re)*half ! 2 up re
1032        rho2ur = rho2r1(ifft_im,2)                   + xccc3d2(ifft_im)*half ! 2 up im
1033        rho2dr = rho2r1(ifft_re,1)-rho2r1(ifft_re,2) + xccc3d2(ifft_re)*half ! 2 down re
1034        rho2dr = rho2r1(ifft_im,1)-rho2r1(ifft_im,2) + xccc3d2(ifft_im)*half ! 2 down im
1035        rho3ur = rho3r1(ifft_re,2)                   + xccc3d3(ifft_re)*half ! 3 up re
1036        rho3ur = rho3r1(ifft_im,2)                   + xccc3d3(ifft_im)*half ! 3 up im
1037        rho3dr = rho3r1(ifft_re,1)-rho3r1(ifft_re,2) + xccc3d3(ifft_re)*half ! 3 down re
1038        rho3dr = rho3r1(ifft_im,1)-rho3r1(ifft_im,2) + xccc3d3(ifft_im)*half ! 3 down im
1039 !      Real part:
1040 !                         uuu                                          uud
1041        xc_tmp(ifft_re,1)= k3xc(ifft,1)*(rho2ur*rho3ur-rho2ui*rho3ui) + k3xc(ifft,2)*(rho2ur*rho3dr-rho2ui*rho3di) + &
1042 !                         udu                                          udd
1043 &       k3xc(ifft,2)*(rho2dr*rho3ur-rho2di*rho3ui) + k3xc(ifft,3)*(rho2dr*rho3dr-rho2di*rho3di)
1044 !                         duu                                          dud
1045        xc_tmp(ifft_re,2)= k3xc(ifft,2)*(rho2ur*rho3ur-rho2ui*rho3ui) + k3xc(ifft,3)*(rho2ur*rho3dr-rho2ui*rho3di) + &
1046 !                         ddu                                          ddd
1047 &       k3xc(ifft,3)*(rho2dr*rho3ur-rho2di*rho3ui) + k3xc(ifft,4)*(rho2dr*rho3dr-rho2di*rho3di)
1048 !      Imaginary part:
1049 !                         uuu                                          uud
1050        xc_tmp(ifft_im,1)= k3xc(ifft,1)*(rho2ur*rho3ui+rho2ui*rho3ur) + k3xc(ifft,2)*(rho2ur*rho3di+rho2ui*rho3dr) + &
1051 !                         udu                                          udd
1052 &       k3xc(ifft,2)*(rho2dr*rho3ui+rho2di*rho3ur) + k3xc(ifft,3)*(rho2dr*rho3di+rho2di*rho3dr)
1053 !                         duu                                          dud
1054        xc_tmp(ifft_im,2)= k3xc(ifft,2)*(rho2ur*rho3ui+rho2ui*rho3ur) + k3xc(ifft,3)*(rho2ur*rho3di+rho2ui*rho3dr) + &
1055 !                         ddu                                          ddd
1056 &       k3xc(ifft,3)*(rho2dr*rho3ui+rho2di*rho3ur) + k3xc(ifft,4)*(rho2dr*rho3di+rho2di*rho3dr)
1057      end do
1058 
1059    end if
1060 
1061  else
1062    MSG_BUG('DFPTNL_PERT is implemented only for nspden=1 or 2')
1063  end if
1064 
1065  call dotprod_vn(cplex,rho1r1_tot,exc3(1),exc3(2),nfftf,nfftotf,nspden,2,xc_tmp,ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1066  ABI_DEALLOCATE(xc_tmp)
1067  ABI_DEALLOCATE(rho1r1_tot)
1068 
1069  exc3_paw = zero
1070  if (usepaw==1) then
1071 
1072    call paw_dfptnl_energy(exc3_paw,dtset%ixc,natom,natom,psps%ntypat,paw_an0,pawang,dtset%pawprtvol,pawrad,&
1073 &   pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,pawtab,dtset%pawxcdev,mpi_enreg%my_atmtab,mpi_enreg%comm_atom)
1074 
1075  end if
1076 
1077 ! **************************************************************************************************
1078 !      Compute E_Hxc^(2:1)
1079 ! **************************************************************************************************
1080 
1081  eHxc21_paw = zero
1082  eHxc21_nhat = zero
1083  if (compute_rho21) then
1084 
1085    if (pawfgr%nfft/=nfftf) then
1086      write(msg,'(2(a,i10))') 'pawfgr%nfft/=nfftf : pawfgr%nfft=',pawfgr%nfft,' nfftf = ',nfftf
1087      MSG_ERROR(msg)
1088    end if
1089 
1090    call pawnhatfr(0,idir_phon,ipert_phon,natom,dtset%natom,nspden,psps%ntypat,&
1091 &   pawang,pawfgrtab,pawrhoij11,pawtab,rprimd)
1092 
1093    ABI_ALLOCATE(nhat21,(cplex*nfftf,nspden))
1094    call pawmkrho(0,arg,cplex,gs_hamkq%gprimd,idir_phon,indsy1,ipert_phon,mpi_enreg,&
1095 &   natom,natom,nspden,nsym1,psps%ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
1096 &   dtset%pawprtvol,pawrhoij21,pawrhoij21_unsym,pawtab,dtset%qptn,dummy_array2,dummy_array2,&
1097 &   dummy_array2,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,&
1098 &   pawang_sym=pawang1,pawnhat=nhat21,pawrhoij0=pawrhoij0)
1099 
1100 !   if (paral_atom) then
1101 !     call pawrhoij_free(pawrhoij21_unsym)
1102 !     ABI_DATATYPE_DEALLOCATE(pawrhoij21_unsym)
1103 !   end if
1104    nzlmopt = 0
1105    call pawdfptenergy(eHxc21_paw,i2pert,ipert_phon,dtset%ixc,natom,dtset%natom,dtset%ntypat,&
1106 &   nzlmopt,nzlmopt,paw_an0,paw_an1_i2pert,paw_ij1_i2pert,pawang,dtset%pawprtvol,&
1107 &   pawrad,pawrhoij1_i2pert,pawrhoij21,pawtab,dtset%pawxcdev,dtset%xclevel)
1108 !&   mpi_atmtab=my_atmtab,comm_atom=my_comm_atom
1109 
1110    ABI_ALLOCATE(v_i2pert,(cplex*nfftf,nspden))
1111    v_i2pert(:,1) = vhartr1_i2pert(:)
1112    if(nspden>1) then
1113      v_i2pert(:,2) = vhartr1_i2pert(:)
1114    end if
1115 
1116    typat_ipert_phon = dtset%typat(ipert_phon)
1117    if (debug_mode/=0) then
1118      write(msg,'(2(a,i6))') ' DFPTNL_PERT : pawtab(',typat_ipert_phon,')%usexcnhat = ',pawtab(ipert_phon)%usexcnhat
1119      call wrtout(std_out,msg,'COLL')
1120    end if
1121    if (pawtab(typat_ipert_phon)%usexcnhat>0) then
1122      v_i2pert(:,:) = v_i2pert(:,:) + vxc1_i2pert(:,:)
1123    end if
1124 
1125    call dotprod_vn(cplex,nhat21,eHxc21_nhat(1),eHxc21_nhat(2),nfftf,nfftotf,nspden,2,v_i2pert,&
1126 &   ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft)
1127 
1128    ABI_DEALLOCATE(v_i2pert)
1129    ABI_DEALLOCATE(nhat21)
1130 
1131  end if
1132 
1133 ! **************************************************************************************************
1134 !    ALL TERMS HAVE BEEN COMPUTED
1135 ! **************************************************************************************************
1136 
1137 !Real part
1138  e3tot = sum_psi1H1psi1 + sum_lambda1psi1psi1 + sum_lambda1psi0S1psi1 + sum_psi0H2psi1a + sum_psi0H2psi1b
1139  e3tot = e3tot + half * (eHxc21_paw(1)+eHxc21_nhat(1)) + sixth * (exc3(1) + exc3_paw(1))
1140 
1141 !Imaginary part
1142  sumi = sum_psi1H1psi1_i + sum_lambda1psi1psi1_i + sum_lambda1psi0S1psi1_i + sum_psi0H2psi1a_i + sum_psi0H2psi1b_i
1143  sumi = sumi   + half * (eHxc21_paw(2)+eHxc21_nhat(2)) + sixth * (exc3(2) + exc3_paw(2))
1144 
1145 !Real parts
1146  d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)   = e3tot
1147  d3etot_1(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi1H1psi1
1148  d3etot_2(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi1psi1
1149  d3etot_3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi0S1psi1
1150  d3etot_4(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1a
1151  d3etot_5(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1b
1152  d3etot_6(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_paw(1)
1153  d3etot_7(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_nhat(1)
1154  d3etot_8(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3(1)
1155  d3etot_9(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3_paw(1)
1156 !Imaginary parts
1157  d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)   = sumi
1158  d3etot_1(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi1H1psi1_i
1159  d3etot_2(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi1psi1_i
1160  d3etot_3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_lambda1psi0S1psi1_i
1161  d3etot_4(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1a_i
1162  d3etot_5(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sum_psi0H2psi1b_i
1163  d3etot_6(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_paw(2)
1164  d3etot_7(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = half * eHxc21_nhat(2)
1165  d3etot_8(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3(2)
1166  d3etot_9(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sixth * exc3_paw(2)
1167 
1168 !Before printing, set small contributions to zero
1169 !Real parts
1170  if (abs(sum_psi1H1psi1)       <tol8) sum_psi1H1psi1        = zero
1171  if (abs(sum_lambda1psi1psi1)  <tol8) sum_lambda1psi1psi1   = zero
1172  if (abs(sum_lambda1psi0S1psi1)<tol8) sum_lambda1psi0S1psi1 = zero
1173  if (abs(sum_psi0H2psi1a)      <tol8) sum_psi0H2psi1a       = zero
1174  if (abs(sum_psi0H2psi1b)      <tol8) sum_psi0H2psi1b       = zero
1175  if (abs(eHxc21_paw(1))        <tol8) eHxc21_paw(1)         = zero
1176  if (abs(eHxc21_nhat(1))       <tol8) eHxc21_nhat(1)        = zero
1177  if (abs(exc3(1))              <tol8) exc3(1)               = zero
1178  if (abs(exc3_paw(1))          <tol8) exc3_paw(1)           = zero
1179  if (abs(e3tot)                <tol8) e3tot                 = zero
1180 
1181 !Imaginary parts
1182  if (abs(sum_psi1H1psi1_i)       <tol8) sum_psi1H1psi1_i        = zero
1183  if (abs(sum_lambda1psi1psi1_i)  <tol8) sum_lambda1psi1psi1_i   = zero
1184  if (abs(sum_lambda1psi0S1psi1_i)<tol8) sum_lambda1psi0S1psi1_i = zero
1185  if (abs(sum_psi0H2psi1a_i)      <tol8) sum_psi0H2psi1a_i       = zero
1186  if (abs(sum_psi0H2psi1b_i)      <tol8) sum_psi0H2psi1b_i       = zero
1187  if (abs(eHxc21_paw(2))          <tol8) eHxc21_paw(2)           = zero
1188  if (abs(eHxc21_nhat(2))         <tol8) eHxc21_nhat(2)          = zero
1189  if (abs(exc3(2))                <tol8) exc3(2)                 = zero
1190  if (abs(exc3_paw(2))            <tol8) exc3_paw(2)             = zero
1191  if (abs(sumi)                   <tol8) sumi                    = zero
1192 
1193  write(msg,'(2a,3(a,i2,a,i1))') ch10,'NONLINEAR : ',&
1194  ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir
1195  call wrtout(std_out,msg,'COLL')
1196  call wrtout(ab_out,msg,'COLL')
1197  if (dtset%nonlinear_info>0) then
1198    write(msg,'(10(a,2(a,f18.8)),a)') &
1199     ch10,'        sum_psi1H1psi1 = ',sum_psi1H1psi1,        ',',sum_psi1H1psi1_i,&
1200     ch10,'   sum_lambda1psi1psi1 = ',sum_lambda1psi1psi1,   ',',sum_lambda1psi1psi1_i,&
1201     ch10,' sum_lambda1psi0S1psi1 = ',sum_lambda1psi0S1psi1, ',',sum_lambda1psi0S1psi1_i,&
1202     ch10,'       sum_psi0H2psi1a = ',sum_psi0H2psi1a,       ',',sum_psi0H2psi1a_i,&
1203     ch10,'       sum_psi0H2psi1b = ',sum_psi0H2psi1b,       ',',sum_psi0H2psi1b_i,&
1204     ch10,'          eHxc21_paw/2 = ',half*eHxc21_paw(1),    ',',half*eHxc21_paw(2),&
1205     ch10,'         eHxc21_nhat/2 = ',half*eHxc21_nhat(1),   ',',half*eHxc21_nhat(2),&
1206     ch10,'                exc3/6 = ',sixth*exc3(1),         ',',sixth*exc3(2),&
1207     ch10,'            exc3_paw/6 = ',sixth*exc3_paw(1),     ',',sixth*exc3_paw(2),&
1208     ch10,' >>>>>>>>>>>>>>> e3tot = ',e3tot,                 ',',sumi,ch10
1209    call wrtout(std_out,msg,'COLL')
1210    call wrtout(ab_out,msg,'COLL')
1211  end if
1212 
1213  if (compute_rho21) then
1214    call pawcprj_free(cwaveprj0)
1215    call pawcprj_free(cwaveprj1)
1216    call pawrhoij_free(pawrhoij21)
1217  end if
1218  ABI_DEALLOCATE(chi_ij)
1219  ABI_DEALLOCATE(dummy_array)
1220  ABI_DEALLOCATE(dummy_array2)
1221  ABI_DEALLOCATE(phkxred)
1222  ABI_DATATYPE_DEALLOCATE(cwaveprj0)
1223  ABI_DATATYPE_DEALLOCATE(cwaveprj1)
1224  ABI_DATATYPE_DEALLOCATE(pawrhoij21)
1225  ABI_DEALLOCATE(ylm)
1226  ABI_DEALLOCATE(ylm1)
1227  ABI_DEALLOCATE(ylmgr)
1228  ABI_DEALLOCATE(ylmgr1)
1229  ABI_DEALLOCATE(vlocal)
1230  ABI_DEALLOCATE(vlocal1_i2pert)
1231  ABI_DEALLOCATE(wfraug)
1232 
1233  call status(0,dtfil%filstat,iexit,level,'exit          ')
1234 
1235  DBG_EXIT("COLL")
1236 
1237 end subroutine dfptnl_pert

ABINIT/m_dfptnl_pert [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfptnl_pert

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

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