TABLE OF CONTENTS


ABINIT/d2frnl [ Functions ]

[ Top ] [ Functions ]

NAME

 d2frnl

FUNCTION

 Compute the frozen-wavefunction non-local contribution for response functions
 (strain and/or phonon)

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of WF
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
  dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere
  has_allddk= True if all ddk file are present on disk
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis
   primitive translations
  mgfftf=maximum size of 1D FFTs for the fine FFT grid (PAW)
  mpi_enreg=information about MPI parallelization
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  my_natom=number of atoms treated by current processor
  natom=number of atoms in unit cell
  nfftf= -PAW ONLY- number of FFT grid points for the fine grid
         (nfftf=nfft for norm-conserving potential runs)
  ngfft(18)=contain all needed information about 3D FFT,
     see ~abinit/doc/variables/vargs.htm#ngfft
  ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
              (ngs_rbzfftf=ngfft for norm-conserving potential runs)
  npwarr(nkpt)=number of planewaves at each k point, and boundary
  ntypat=integer specification of atom type (1, 2, ...)
  occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2) at each k point
  rfphon=1   if non local contribution of dynamical matrix have to be computed
  rfstrs!=0  if non local contribution of elastic tensor have to be computed
  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawbec= flag for the computation of Born Effective Charge within PAW ; set to 1 if yes
  pawpiezo= flag for the computation of piezoelectric tensor  within PAW ; set to 1 if yes
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=phase information related to structure factor
  ph1df(2,3*(2*mgfftf+1)*natom)=phase information related to structure factor on the fine FFT grid (PAW)
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless)
  vtrial(nfftf,nspden)=total potential (Hartree+XC+loc)
  vxc(nfftf,nspden)=XC potential
  xred(3,natom)=reduced coordinates of atoms (dimensionless)
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylmgr(mpw*mkmem,9,mpsang*mpsang*useylm)= gradients of real spherical harmonics for each G and k point

OUTPUT

  becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only)
                            (3,natom) = derivative wr to the displ. of one atom in one direction
                            (3)       = derivative wr to electric field in one direction
  piezofrnl(3,6*pawpiezo)=NL frozen contribution to piezoelectric tensor (PAW only)
  dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
         non-symmetrized non-local contribution to the dynamical matrix
         If NCPP, it depends on one atom
         If PAW,  it depends on two atoms
  eltfrnl(6+3*natom,6)=non-symmetrized non-local contribution to the
                    elastic tensor

SIDE EFFECTS

  ===== if psps%usepaw==1
  pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
                          pawfgrtab(:)%gylmgr2 are deallocated here
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
    (gradients of rhoij for each atom with respect to atomic positions are computed here)

PARENTS

      respfn

CHILDREN

      appdig,check_degeneracies,destroy_hamiltonian,dotprod_g
      init_hamiltonian,load_k_hamiltonian,load_spin_hamiltonian,metric,mkffnl
      mkkin,mkkpg,nonlop,paw_ij_free,paw_ij_init,paw_ij_nullify
      paw_ij_reset_flags,pawaccrhoij,pawcprj_alloc,pawcprj_free,pawdij2e1kb
      pawdijfr,pawfgrtab_free,pawfgrtab_init,pawgrnl,pawrhoij_free
      pawrhoij_gather,pawrhoij_nullify,pawtab_get_lsize,strconv,symrhoij
      timab,wfk_close,wfk_open_read,wfk_read_bks,wrtout,xmpi_sum

SOURCE

 126 subroutine d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,dyfr_nondiag,efmasdeg,efmasval,eigen,eltfrnl,&
 127 &          gsqcut,has_allddk,indsym,kg,mgfftf,mpi_enreg,mpsang,my_natom,natom,nfftf,ngfft,ngfftf,npwarr,&
 128 &          occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,&
 129 &          rprimd,rfphon,rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr)
 130 
 131  use defs_basis
 132  use defs_datatypes
 133  use defs_abitypes
 134  use m_abicore
 135  use m_xmpi
 136  use m_errors
 137  use m_cgtools
 138  use m_nctk
 139  use m_hamiltonian
 140  use m_efmas_defs
 141  use m_wfk
 142 
 143  use m_time,     only : timab
 144  use m_geometry, only : metric, strconv
 145  use m_efmas,    only : check_degeneracies
 146  use m_io_tools, only : file_exists
 147  use m_hdr,      only : hdr_skip
 148  use m_pawang,   only : pawang_type
 149  use m_pawrad,   only : pawrad_type
 150  use m_pawtab,   only : pawtab_type,pawtab_get_lsize
 151  use m_pawfgrtab,only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
 152  use m_paw_ij,   only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags
 153  use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free, pawrhoij_gather, pawrhoij_nullify, symrhoij
 154  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_copy, pawcprj_free
 155  use m_pawdij,   only : pawdijfr
 156  use m_paw_dfpt, only : pawgrnl
 157  use m_kg,       only : mkkin, mkkpg
 158  use m_mkffnl,   only : mkffnl
 159  use m_mpinfo,   only : proc_distrb_cycle
 160  use m_nonlop,   only : nonlop
 161  use m_paw_occupancies, only : pawaccrhoij
 162 
 163 !This section has been created automatically by the script Abilint (TD).
 164 !Do not modify the following lines by hand.
 165 #undef ABI_FUNC
 166 #define ABI_FUNC 'd2frnl'
 167 !End of the abilint section
 168 
 169  implicit none
 170 
 171 !Arguments ------------------------------------
 172 !scalars
 173  integer,intent(in) :: dyfr_cplex,dyfr_nondiag,mgfftf,mpsang,my_natom,natom
 174  integer,intent(in) :: nfftf,pawbec,pawpiezo,rfphon,rfstrs
 175  real(dp),intent(in) :: gsqcut
 176  type(MPI_type),intent(in) :: mpi_enreg
 177  type(datafiles_type),intent(in) :: dtfil
 178  type(dataset_type),intent(in) :: dtset
 179  type(pawang_type),intent(in) :: pawang
 180  type(pseudopotential_type),intent(in) :: psps
 181 !arrays
 182  integer,intent(in) :: indsym(4,dtset%nsym,natom),kg(3,dtset%mpw*dtset%mkmem)
 183  integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt)
 184  integer,intent(in) :: symrec(3,3,dtset%nsym)
 185  real(dp),intent(in) :: cg(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol)
 186  real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
 187  real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 188  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*natom)
 189  real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*natom),rprimd(3,3)
 190  real(dp),intent(in) :: vxc(nfftf,dtset%nspden),xred(3,natom)
 191  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm)
 192  real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,9,mpsang*mpsang*psps%useylm)
 193  real(dp),intent(in),target :: vtrial(nfftf,dtset%nspden)
 194  real(dp),intent(out) :: becfrnl(3,natom,3*pawbec),piezofrnl(6,3*pawpiezo)
 195  real(dp),intent(out) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)
 196  real(dp),intent(out) :: eltfrnl(6+3*natom,6)
 197  logical,intent(inout):: has_allddk
 198  type(efmasdeg_type),allocatable,intent(out):: efmasdeg(:)
 199  type(efmasval_type),allocatable,intent(out):: efmasval(:,:)
 200  type(paw_ij_type),intent(in) :: paw_ij(my_natom)
 201  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
 202  type(pawrad_type),intent(in) :: pawrad(psps%ntypat)
 203  type(pawrhoij_type),intent(inout),target :: pawrhoij(my_natom*psps%usepaw)
 204  type(pawtab_type),intent(in) :: pawtab(psps%ntypat)
 205 
 206 !Local variables-------------------------------
 207 !scalars
 208  integer,parameter :: formeig1=1,usecprj=0
 209  integer :: bandmin,bandmax,bdtot_index,bufdim
 210  integer :: choice_bec2,choice_bec54,choice_efmas,choice_phon,choice_strs,choice_piez3,choice_piez55
 211  integer :: cplex,cplx,cpopt,cpopt_bec,ddkcase,deg_dim
 212  integer :: dimffnl,dimffnl_str,dimnhat,ia,iatom,iashift,iband,jband,ibg,icg,icplx,ideg,ider,idir
 213  integer :: ider_str,idir_ffnl,idir_str,ielt,ieltx,ierr,ii,ikg,ikpt,ilm,ipw
 214  integer :: ispinor,isppol,istwf_k,isub,itypat,jj,jsub,klmn,master,me,mu
 215  integer :: my_comm_atom,n1,n2,n3,nband_k,ncpgr,nfftot,ngrhoij,nkpg,nnlout_bec1,nnlout_bec2,nnlout_efmas
 216  integer :: nnlout_piez1,nnlout_piez2,nnlout_phon,nnlout_strs,npw_,npw_k,nsp,nsploop,nu
 217  integer :: optgr,optgr2,option,option_rhoij,optstr,optstr2,paw_opt,paw_opt_1,paw_opt_3,paw_opt_efmas
 218  integer :: shift_rhoij,signs,signs_field,spaceworld,sz2,sz3,tim_nonlop
 219  real(dp) :: arg,eig_k,enl,enlk,occ_k,ucvol,wtk_k
 220  logical :: has_ddk_file,need_becfr,need_efmas,need_piezofr,paral_atom,t_test,usetimerev
 221  character(len=500) :: msg
 222  type(gs_hamiltonian_type) :: gs_ham
 223 !arrays
 224  integer :: ik_ddk(3),ddkfil(3)
 225  integer,allocatable :: dimlmn(:),kg_k(:,:),l_size_atm(:)
 226  integer,pointer :: my_atmtab(:)
 227  real(dp) :: dotprod(2),dummy(0),gmet(3,3),gprimd(3,3),grhoij(3),kpoint(3),nonlop_dum(1,1)
 228  real(dp) :: rmet(3,3),tsec(2)
 229  real(dp),allocatable :: becfrnl_tmp(:,:,:),becfrnlk(:,:,:),becij(:,:,:,:,:),cg_left(:,:)
 230  real(dp),allocatable :: cwavef(:,:),ddk(:,:),ddkinpw(:,:,:),dyfrnlk(:,:)
 231  real(dp),allocatable :: elt_work(:,:),eltfrnlk(:,:),enlout_bec1(:),enlout_bec2(:),enlout_efmas(:)
 232  real(dp),allocatable :: enlout_piez1(:),enlout_piez2(:),enlout_phon(:),enlout_strs(:)
 233  real(dp),allocatable :: gh2c(:,:),gs2c(:,:)
 234  real(dp),allocatable :: kpg_k(:,:),mpibuf(:),nhat_dum(:,:),piezofrnlk(:,:),ph3d(:,:,:)
 235  real(dp),allocatable :: svectout(:,:),ylm_k(:,:),ylmgr_k(:,:,:)
 236  real(dp),allocatable,target :: ffnl(:,:,:,:),ffnl_str(:,:,:,:,:)
 237  character(len=fnlen) :: fiwfddk(3)
 238  type(paw_ij_type),allocatable :: paw_ij_tmp(:)
 239  type(pawcprj_type),allocatable,target :: cwaveprj(:,:)
 240  type(pawfgrtab_type),allocatable :: pawfgrtab_tmp(:)
 241  type(pawrhoij_type),pointer :: pawrhoij_tot(:)
 242  type(wfk_t) :: ddkfiles(3)
 243 
 244 ! *************************************************************************
 245 
 246  DBG_ENTER("COLL")
 247 
 248  call timab(159,1,tsec)
 249 
 250  write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivatives, this can take some time...',ch10
 251  call wrtout(std_out,msg,'COLL')
 252 
 253 !Set up parallelism
 254  spaceworld=mpi_enreg%comm_cell
 255  me=mpi_enreg%me_kpt
 256  master=0
 257  paral_atom=(my_natom/=natom)
 258  my_comm_atom=mpi_enreg%comm_atom
 259  my_atmtab=>mpi_enreg%my_atmtab
 260 
 261 !Compute gmet, gprimd and ucvol from rprimd
 262  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
 263 
 264 !If needed, check for ddk files (used for effective charges)
 265  if (pawbec==1.or.pawpiezo==1) then
 266    ddkfil(:)=0
 267    do ii=1,3
 268      ddkcase=ii+natom*3
 269      call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(ii))
 270      t_test = file_exists(fiwfddk(ii))
 271      ! Trick needed to run Abinit test suite in netcdf mode.
 272      if (.not. t_test .and. file_exists(nctk_ncify(fiwfddk(ii)))) then
 273        t_test = .True.; fiwfddk(ii) = nctk_ncify(fiwfddk(ii))
 274        write(msg,"(3a)")"- File: ",trim(fiwfddk(ii))," does not exist but found netcdf file with similar name."
 275        call wrtout(std_out,msg,'COLL')
 276      end if
 277      if (t_test) ddkfil(ii)=20+ii ! Note the use of unit numbers 21, 22 and 23
 278    end do
 279    has_ddk_file=(any(ddkfil(:)>0))
 280    has_allddk  =(all(ddkfil(:)>0))
 281  else
 282    has_ddk_file=.FALSE.
 283    has_allddk  =.FALSE.
 284  end if
 285 
 286  if(pawbec==1.or.pawpiezo==1.and.has_ddk_file) then
 287    if(.not.has_allddk) then
 288      write(msg,'(5a)')ch10,&
 289 &     ' WARNING: All ddk perturbations are needed to compute',ch10,&
 290 &     ' the frozen part of Born effective charges and/or piezoelectric tensor.',ch10
 291      call wrtout(std_out,msg,'COLL')
 292    else
 293      write(msg,'(5a)')ch10,&
 294 &     ' All ddk perturbations are available.',ch10,&
 295 &     ' The frozen part of Born effective charges and/or piezoelectric tensor will be computed',ch10
 296      call wrtout(std_out,msg,'COLL')
 297    end if
 298  end if
 299 
 300  need_becfr=(pawbec==1.and.has_ddk_file)
 301  need_piezofr=(pawpiezo==1.and.has_ddk_file)
 302 
 303 !Initialization of frozen non local array
 304  if(rfphon==1) then
 305    dyfrnl(:,:,:,:,:)=zero
 306    ABI_ALLOCATE(dyfrnlk,(6,natom))
 307  end if
 308  if(rfstrs/=0)then
 309    eltfrnl(:,:)=zero;enl=zero
 310    ABI_ALLOCATE(eltfrnlk,(6+3*natom,6))
 311  end if
 312  if (need_becfr) then
 313    becfrnl(:,:,:)=zero
 314    ABI_ALLOCATE(becfrnlk,(3,natom,3))
 315  end if
 316  if (need_piezofr) then
 317    piezofrnl(:,:)=zero
 318    ABI_ALLOCATE(piezofrnlk,(6,3))
 319  end if
 320  need_efmas=dtset%efmas>0
 321  if(need_efmas.and.(rfphon==1.or.rfstrs/=0.or.need_becfr.or.need_piezofr)) then
 322    write(msg,'(5a)')ch10,&
 323 &   ' ERROR: Efmas calculation is incompatible with phonons, elastic tensor, Born effective charges,',ch10,&
 324 &   ' and piezoelectric tensor calculations. Please revise your input.',ch10
 325    MSG_ERROR(msg)
 326  end if
 327 
 328 !Common initialization
 329  bdtot_index=0;ibg=0;icg=0
 330  nsploop=dtset%nsppol;if (dtset%nspden==4) nsploop=4
 331  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
 332  nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
 333 
 334 !Common data for "nonlop" routine
 335  tim_nonlop=6
 336  signs=1 ; signs_field = 2 ; eig_k=zero ; idir=0
 337 ! ffnl are in cartesian coordinates for EFMAS (idir==4),
 338 ! in contrast to reduced coordinates for the other responses (idir==0).
 339  idir_ffnl=0 ; if(need_efmas) idir_ffnl=4
 340  choice_phon=0;choice_strs=0
 341  if(rfphon==1)then
 342    shift_rhoij=0
 343    choice_phon=4
 344    nnlout_phon=max(1,6*natom)
 345    ABI_ALLOCATE(enlout_phon,(nnlout_phon))
 346  end if
 347  if(rfstrs/=0)then
 348    shift_rhoij=6
 349    choice_strs=6
 350    nnlout_strs=6*(3*natom+6)
 351    ABI_ALLOCATE(enlout_strs,(nnlout_strs))
 352  end if
 353  if (psps%usepaw==0) then
 354    paw_opt=0 ; cpopt=-1
 355  else
 356    paw_opt=2 ; cpopt=1+2*usecprj
 357  end if
 358  if(need_piezofr)then
 359    choice_piez3  =  3
 360    choice_piez55 = 55
 361    nnlout_piez1  =  6
 362    nnlout_piez2  = 36
 363    paw_opt_1     = 1
 364    paw_opt_3     = 3
 365    ABI_ALLOCATE(enlout_piez1,(nnlout_piez1))
 366    ABI_ALLOCATE(enlout_piez2,(nnlout_piez2))
 367  end if
 368  if (need_becfr) then
 369    choice_bec2=2 ; choice_bec54=54
 370    nnlout_bec1=max(1,3*natom) ; nnlout_bec2=max(1,18*natom);
 371    paw_opt_1=1 ; paw_opt_3=3 ; cpopt_bec=-1
 372    ABI_ALLOCATE(enlout_bec1,(nnlout_bec1))
 373    ABI_ALLOCATE(enlout_bec2,(nnlout_bec2))
 374  else
 375    choice_bec2=0 ; choice_bec54=0
 376    nnlout_bec1=0
 377  end if
 378  if(need_efmas) then
 379    ABI_MALLOC(enlout_efmas,(0))
 380    ABI_DATATYPE_ALLOCATE(efmasdeg,(dtset%nkpt))
 381    ABI_DATATYPE_ALLOCATE(efmasval,(dtset%mband,dtset%nkpt))
 382  end if
 383 
 384 !Initialize Hamiltonian (k-independent terms)
 385  call init_hamiltonian(gs_ham,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
 386 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
 387 & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
 388 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom)
 389 
 390 !===== PAW specific section
 391  if (psps%usepaw==1) then
 392 
 393 !  Define several sizes & flags
 394    ncpgr=0;ngrhoij=0
 395    if(rfphon==1)then
 396      ncpgr=3;ngrhoij=3
 397    end if
 398    if(rfphon==1.or.need_becfr)then
 399      ncpgr=6;ngrhoij=6
 400    end if
 401    if(rfstrs/=0.and.rfphon==1)then
 402      ncpgr=9;ngrhoij=9
 403    end if
 404    if(rfstrs/=0.or.need_piezofr)then
 405      ncpgr=9;ngrhoij=9
 406    end if
 407 
 408 !  If PAW and Born Eff. Charges, one has to compute some additional data:
 409 !  For each atom and for electric field direction k:
 410 !  becij(k)=<Phi_i|r_k-R_k|Phi_j>-<tPhi_i|r_k-R_k|tPhi_j> + sij.R_k
 411    if (need_becfr.or.need_piezofr) then
 412      ABI_ALLOCATE(becij,(gs_ham%dimekb1,gs_ham%dimekb2,dtset%nspinor**2,1,3))
 413      becij=zero
 414      ABI_DATATYPE_ALLOCATE(paw_ij_tmp,(my_natom))
 415      ABI_DATATYPE_ALLOCATE(pawfgrtab_tmp,(my_natom))
 416      call paw_ij_nullify(paw_ij_tmp)
 417      cplex=1;nsp=1 ! Force nsppol/nspden to 1 because Dij^(1) due to electric field is spin-independent
 418      call paw_ij_init(paw_ij_tmp,cplex,dtset%nspinor,nsp,nsp,dtset%pawspnorb,natom,psps%ntypat,&
 419 &     dtset%typat,pawtab,has_dijfr=1,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab )
 420      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,mpi_atmtab=my_atmtab)
 421      call pawfgrtab_init(pawfgrtab_tmp,1,l_size_atm,dtset%nspden,dtset%typat,&
 422 &     mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
 423      ABI_DEALLOCATE(l_size_atm)
 424      do ii=1,3 ! Loop over direction of electric field
 425        call paw_ij_reset_flags(paw_ij_tmp,all=.True.)
 426        call pawdijfr(cplex,gprimd,ii,natom+2,my_natom,natom,nfftf,ngfftf,nsp,nsp,psps%ntypat,&
 427 &       0,paw_ij_tmp,pawang,pawfgrtab_tmp,pawrad,pawtab,&
 428 &       (/zero,zero,zero/),rprimd,ucvol,vtrial,vtrial,vxc,xred,&
 429 &       comm_atom=my_comm_atom, mpi_atmtab=my_atmtab ) ! vtrial not used here
 430        do isppol=1,dtset%nspinor**2
 431          call pawdij2e1kb(paw_ij_tmp(:),nsp,my_comm_atom,e1kbfr=becij(:,:,:,:,ii),mpi_atmtab=my_atmtab)
 432        end do
 433      end do
 434      call paw_ij_free(paw_ij_tmp)
 435      call pawfgrtab_free(pawfgrtab_tmp)
 436      ABI_DATATYPE_DEALLOCATE(paw_ij_tmp)
 437      ABI_DATATYPE_DEALLOCATE(pawfgrtab_tmp)
 438    end if
 439 
 440 !  PAW occupancies: need to communicate when paral atom is activated
 441    if (paral_atom) then
 442      ABI_DATATYPE_ALLOCATE(pawrhoij_tot,(natom))
 443      call pawrhoij_nullify(pawrhoij_tot)
 444      call pawrhoij_gather(pawrhoij,pawrhoij_tot,-1,my_comm_atom)
 445    else
 446      pawrhoij_tot => pawrhoij
 447    end if
 448 
 449 !  Projected WF (cprj) and PAW occupancies (& gradients)
 450    ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,dtset%nspinor))
 451    call pawcprj_alloc(cwaveprj,ncpgr,gs_ham%dimcprj)
 452    do iatom=1,natom
 453      sz2=pawrhoij_tot(iatom)%cplex*pawrhoij_tot(iatom)%lmn2_size
 454      sz3=pawrhoij_tot(iatom)%nspden
 455      ABI_ALLOCATE(pawrhoij_tot(iatom)%grhoij,(ngrhoij,sz2,sz3))
 456      pawrhoij_tot(iatom)%ngrhoij=ngrhoij
 457      pawrhoij_tot(iatom)%grhoij=zero
 458    end do
 459    usetimerev=(dtset%kptopt>0.and.dtset%kptopt<3)
 460 
 461  else
 462    ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0))
 463  end if !PAW
 464 
 465 !If needed, manage ddk files
 466 !Open ddk WF file(s)
 467  if (need_becfr.or.need_piezofr) then
 468    do ii=1,3 ! Loop over elect. field directions
 469      if (ddkfil(ii)/=0) then
 470        write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(ii))
 471        call wrtout(std_out,msg,'COLL')
 472        call wrtout(ab_out,msg,'COLL')
 473        call wfk_open_read(ddkfiles(ii),fiwfddk(ii),formeig1,dtset%iomode,ddkfil(ii),spaceworld)
 474      end if
 475    end do
 476  end if
 477 
 478 !LOOP OVER SPINS
 479  do isppol=1,dtset%nsppol
 480 
 481 !  Continue to initialize the Hamiltonian (PAW DIJ coefficients)
 482    call load_spin_hamiltonian(gs_ham,isppol,with_nonlocal=.true.)
 483 
 484 !  Rewind (k+G) data if needed
 485    ikg=0
 486 
 487 !  Loop over k points
 488    do ikpt=1,dtset%nkpt
 489      nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 490      istwf_k=dtset%istwfk(ikpt)
 491      npw_k=npwarr(ikpt)
 492      wtk_k=dtset%wtk(ikpt)
 493      kpoint(:)=dtset%kptns(:,ikpt)
 494 
 495 !    Skip this k-point if not the proper processor
 496      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
 497        bdtot_index=bdtot_index+nband_k
 498        cycle
 499      end if
 500 
 501 !    If needed, manage ddk files
 502      if (need_becfr.or.need_piezofr) then
 503        do ii=1,3 ! Loop over elect. field directions
 504          if (ddkfil(ii)/=0)then
 505 !        Number of k points to skip in the full set of k pointsp
 506            ik_ddk(ii) = wfk_findk(ddkfiles(ii), kpoint)
 507            ABI_CHECK(ik_ddk(ii) /= -1, "Cannot find k-point in DDK")
 508            npw_ = ddkfiles(ii)%hdr%npwarr(ik_ddk(ii))
 509            if (npw_/=npw_k) then
 510              write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')&
 511 &             'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',ii,ch10,&
 512 &             'the number of plane waves in the ddk file is equal to', npw_,ch10,&
 513 &             'while it should be ',npw_k
 514              MSG_ERROR(msg)
 515            end if
 516 
 517          end if
 518        end do
 519      end if
 520 
 521      ABI_ALLOCATE(cwavef,(2,npw_k*dtset%nspinor))
 522      if (need_becfr.or.need_piezofr) then
 523        ABI_ALLOCATE(svectout,(2,npw_k*dtset%nspinor))
 524      end if
 525      if (need_efmas) then
 526        ABI_MALLOC(cg_left,(2,npw_k*dtset%nspinor))
 527        ABI_MALLOC(gh2c,(2,npw_k*dtset%nspinor))
 528        ABI_MALLOC(gs2c,(2,npw_k*dtset%nspinor))
 529      end if
 530 
 531      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
 532      if(rfstrs/=0.or.need_becfr.or.need_piezofr.or.need_efmas)then
 533        ABI_ALLOCATE(ylmgr_k,(npw_k,9,mpsang*mpsang*psps%useylm))
 534      else
 535        ABI_ALLOCATE(ylmgr_k,(0,0,0))
 536      end if
 537 
 538      ABI_ALLOCATE(kg_k,(3,npw_k))
 539      kg_k(:,:) = 0
 540 !$OMP PARALLEL DO
 541      do ipw=1,npw_k
 542        kg_k(1,ipw)=kg(1,ipw+ikg)
 543        kg_k(2,ipw)=kg(2,ipw+ikg)
 544        kg_k(3,ipw)=kg(3,ipw+ikg)
 545      end do
 546      if (psps%useylm==1) then
 547 !SOMP PARALLEL DO COLLAPSE(2)
 548        do ilm=1,mpsang*mpsang
 549          do ipw=1,npw_k
 550            ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
 551          end do
 552        end do
 553        if(rfstrs/=0.or.need_becfr.or.need_piezofr.or.need_efmas)then
 554 !SOMP PARALLEL DO COLLAPSE(3)
 555          do ilm=1,mpsang*mpsang
 556            do ii=1,9
 557              do ipw=1,npw_k
 558                ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm)
 559              end do
 560            end do
 561          end do
 562        end if
 563      end if
 564 
 565      cplex=2;if (istwf_k>1) cplex=1
 566 
 567 !    Compute (k+G) vectors (only if useylm=1)
 568      nkpg=0
 569      if (rfstrs/=0.or.need_efmas.or.pawpiezo==1) nkpg=3*dtset%nloalg(3)
 570      ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 571      if (nkpg>0) then
 572        call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
 573      end if
 574 
 575      !EFMAS: Compute second order derivatives w/r to k for all direction for this k-point.
 576      if (need_efmas) then
 577        ABI_ALLOCATE(ddkinpw,(npw_k,3,3))
 578        do mu=1,3
 579          do nu=1,3
 580 !           call d2kpg(ddkinpw(:,mu,nu),dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,mu,nu,kg_k,kpoint,npw_k)
 581            call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw(:,mu,nu),kpoint,npw_k,mu,nu)
 582          end do
 583        end do
 584      end if
 585 
 586 !    Compute nonlocal form factors ffnl at all (k+G):
 587      ider=0;dimffnl=1;
 588      if(need_becfr) then
 589        ider=1;dimffnl=4
 590      end if
 591      if(rfstrs/=0.or.need_piezofr.or.need_efmas)then
 592        ider=2;dimffnl=3+7*psps%useylm
 593      end if
 594      ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
 595      call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
 596 &     gmet,gprimd,ider,idir_ffnl,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
 597 &     psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,&
 598 &     psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
 599 
 600 !    For piezoelectric tensor need additional ffnl derivatives
 601      if(need_piezofr)then
 602        ider_str=1 ; dimffnl_str=2
 603        ABI_ALLOCATE(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,psps%ntypat,6))
 604        do mu=1,6 !loop over strain
 605          idir_str=-mu
 606          call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str(:,:,:,:,mu),&
 607 &         psps%ffspl,gmet,gprimd,ider_str,idir_str,psps%indlmn,kg_k,kpg_k,&
 608 &         kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,&
 609 &         psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
 610        end do
 611      end if
 612 
 613 !    Load k-dependent part in the Hamiltonian datastructure
 614      ABI_ALLOCATE(ph3d,(2,npw_k,gs_ham%matblk))
 615      call load_k_hamiltonian(gs_ham,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,&
 616 &     kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.)
 617 
 618 !    Initialize contributions from current k point
 619      if(rfphon==1) dyfrnlk(:,:)=zero
 620      if(rfstrs/=0)then
 621        enlk=zero;eltfrnlk(:,:)=zero
 622      end if
 623      if (need_becfr) becfrnlk(:,:,:)=zero
 624      if (need_piezofr) piezofrnlk(:,:)=zero
 625      if(need_efmas) then
 626        call check_degeneracies(efmasdeg(ikpt),dtset%efmas_bands(:,ikpt),nband_k,eigen(bdtot_index+1:bdtot_index+nband_k), &
 627 &       dtset%efmas_deg_tol)
 628        do ideg=1,efmasdeg(ikpt)%ndegs
 629          if( efmasdeg(ikpt)%deg_range(1) <= ideg .and. ideg <= efmasdeg(ikpt)%deg_range(2) ) then
 630            deg_dim=efmasdeg(ikpt)%degs_bounds(2,ideg) - efmasdeg(ikpt)%degs_bounds(1,ideg) + 1
 631            ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(3,3,deg_dim,deg_dim))
 632            ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(3,3,deg_dim,deg_dim))
 633            efmasval(ideg,ikpt)%ch2c=zero
 634            efmasval(ideg,ikpt)%eig2_diag=zero
 635          else
 636            ABI_MALLOC(efmasval(ideg,ikpt)%ch2c,(0,0,0,0))
 637            ABI_MALLOC(efmasval(ideg,ikpt)%eig2_diag,(0,0,0,0))
 638          end if
 639        end do
 640      end if
 641 
 642 !    Loop over bands
 643      do iband=1,nband_k
 644 
 645        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) then
 646          cycle
 647        end if
 648 
 649        occ_k=occ(iband+bdtot_index)
 650        cwavef(:,1:npw_k*dtset%nspinor) = cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg)
 651 
 652 !      Compute non-local contributions from n,k
 653        if (psps%usepaw==1) eig_k=eigen(iband+bdtot_index)
 654 
 655 !      === Dynamical matrix
 656        if(rfphon==1) then
 657          call nonlop(choice_phon,cpopt,cwaveprj,enlout_phon,gs_ham,idir,(/eig_k/),mpi_enreg,1,&
 658 &         nnlout_phon,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
 659 !        Accumulate non-local contributions from n,k
 660          dyfrnlk(:,:)=dyfrnlk(:,:)+occ_k*reshape(enlout_phon(:),(/6,natom/))
 661        end if
 662 
 663 !      === Elastic tensor
 664        if(rfstrs/=0) then
 665          call nonlop(choice_strs,cpopt,cwaveprj,enlout_strs,gs_ham,idir,(/eig_k/),mpi_enreg,1,&
 666 &         nnlout_strs,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
 667 !        Accumulate non-local contribut ions from n,k
 668          eltfrnlk(:,:)=eltfrnlk(:,:)+occ_k*reshape(enlout_strs(:),(/3*natom+6,6/))
 669        end if !endo if strs
 670 
 671 !      PAW: accumulate gradients of rhoij
 672        !EFMAS: Bug with efmas currently; to be looked into...
 673        if (psps%usepaw==1.and.(.not.need_efmas)) then
 674          call pawaccrhoij(gs_ham%atindx,cplex,cwaveprj,cwaveprj,0,isppol,natom,&
 675 &         natom,dtset%nspinor,occ_k,3,pawrhoij_tot,usetimerev,wtk_k)
 676        end if
 677 
 678 !      PAW: Compute frozen contribution to piezo electric tensor
 679        if (need_piezofr) then
 680          do ii=1,3 ! Loop over elect. field directions
 681            call nonlop(choice_piez3,cpopt,cwaveprj,enlout_piez1,gs_ham,0,(/zero/),mpi_enreg,1,&
 682 &           nnlout_piez1,paw_opt_1,signs,nonlop_dum,tim_nonlop,cwavef,cwavef,enl=becij(:,:,:,:,ii))
 683            piezofrnlk(:,ii)=piezofrnlk(:,ii)+occ_k*enlout_piez1(:)
 684          end do !end do ii
 685        end if
 686 
 687 !      PAW: Compute frozen contribution to Born Effective Charges
 688        if (need_becfr) then
 689          do ii=1,3 ! Loop over elect. field directions
 690            call nonlop(choice_bec2,cpopt,cwaveprj,enlout_bec1,gs_ham,0,(/zero/),mpi_enreg,1,&
 691 &           nnlout_bec1,paw_opt_1,signs,nonlop_dum,tim_nonlop,cwavef,cwavef,enl=becij(:,:,:,:,ii))
 692            becfrnlk(:,:,ii)=becfrnlk(:,:,ii)+occ_k*reshape(enlout_bec1(:),(/3,natom/))
 693          end do !end do ii
 694        end if
 695 
 696        if (need_becfr.or.need_piezofr) then
 697          do ii=1,3 ! Loop over elect. field directions
 698 !          Not able to compute if ipert=(Elect. field) and no ddk WF file
 699            if (ddkfil(ii)==0) cycle
 700 !            Read ddk wave function
 701            ABI_ALLOCATE(ddk,(2,npw_k*dtset%nspinor))
 702            if (ddkfil(ii)/=0) then
 703              call wfk_read_bks(ddkfiles(ii), iband, ik_ddk(ii), isppol, xmpio_single, cg_bks=ddk)
 704 !            Multiply ddk by +i
 705              do jj=1,npw_k*dtset%nspinor
 706                arg=ddk(1,jj)
 707                ddk(1,jj)=-ddk(2,jj);ddk(2,jj)=arg
 708              end do
 709            else
 710              ddk=zero
 711            end if
 712 
 713            if(need_becfr)then
 714              do iatom=1,natom !Loop over atom
 715                ia=gs_ham%atindx(iatom)
 716                do mu=1,3 !loop over atom direction
 717                  call nonlop(choice_bec2,cpopt_bec,cwaveprj,enlout_bec1,gs_ham,mu,(/zero/),&
 718 &                 mpi_enreg,1,nnlout_bec1,paw_opt_3,signs_field,svectout,tim_nonlop,&
 719 &                 cwavef,cwavef,iatom_only=iatom)
 720                  call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,svectout,ddk,&
 721 &                 mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 722                  becfrnlk(mu,ia,ii)=becfrnlk(mu,ia,ii)+occ_k*dotprod(1)
 723                end do
 724              end do
 725            end if
 726 
 727            if(need_piezofr)then
 728              do mu=1,6 !loop over strain
 729                call load_k_hamiltonian(gs_ham,ffnl_k=ffnl_str(:,:,:,:,mu))
 730                call nonlop(choice_piez3,cpopt,cwaveprj,enlout_piez1,gs_ham,mu,(/zero/),mpi_enreg,1,&
 731 &               nnlout_piez1,paw_opt_3,signs_field,svectout,tim_nonlop,cwavef,svectout)
 732                call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,svectout,ddk,&
 733 &               mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
 734                piezofrnlk(mu,ii)=piezofrnlk(mu,ii)+occ_k*dotprod(1)
 735              end do
 736              call load_k_hamiltonian(gs_ham,ffnl_k=ffnl)
 737            end if
 738 
 739            ABI_DEALLOCATE(ddk)
 740          end do ! End loop ddk file
 741        end if
 742 
 743        if(need_piezofr)then
 744          enlout_piez2 = zero
 745          call nonlop(choice_piez55,cpopt,cwaveprj,enlout_piez2,gs_ham,0,(/zero/),mpi_enreg,1,&
 746 &         nnlout_piez2,paw_opt_3,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
 747 !         Multiply enlout by +i
 748          iashift = 1
 749          do mu=1,6     ! strain
 750            do nu=1,3   ! k
 751              piezofrnlk(mu,nu)=piezofrnlk(mu,nu)-occ_k*(enlout_piez2(iashift+1)) ! Real part
 752 !            piezofrnlk(mu,nu)=piezofrnlk(mu,nu)+occ_k*(enlout_piez2(iashift  ))! Imaginary part
 753              iashift = iashift + 2
 754            end do
 755          end do
 756        end if
 757 
 758        if(need_becfr)then
 759          call nonlop(choice_bec54,cpopt,cwaveprj,enlout_bec2,gs_ham,0,(/zero/),mpi_enreg,1,&
 760 &         nnlout_bec2,paw_opt_3,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
 761 !        Multiply enlout by +i
 762          iashift = 1
 763          do iatom=1,natom ! atm
 764            do mu=1,3     ! atm pos.
 765              do nu=1,3   ! k
 766                becfrnlk(mu,iatom,nu)=becfrnlk(mu,iatom,nu)-occ_k*(enlout_bec2(iashift+1)) ! Real part
 767 !               becfrnlk(mu,iatom,nu)=becfrnlk(mu,iatom,nu)+occ_k*(enlout_bec2(iashift  ))! Imaginary part
 768                iashift = iashift + 2
 769              end do
 770            end do
 771          end do
 772        end if
 773 
 774        if(need_efmas) then
 775          bandmin=efmasdeg(ikpt)%degs_bounds(1, efmasdeg(ikpt)%deg_range(1) )
 776          bandmax=efmasdeg(ikpt)%degs_bounds(2, efmasdeg(ikpt)%deg_range(2) )
 777          if ( iband>=bandmin .and. iband<=bandmax ) then
 778            choice_efmas=8; signs=2
 779            cpopt=-1  !To prevent re-use of stored dgxdt, which are not for all direction required for EFMAS.
 780            paw_opt_efmas=0; if(psps%usepaw/=0) paw_opt_efmas=4 !To get both gh2c and gs2c
 781            nnlout_efmas=0; tim_nonlop=0 ! No tim_nonlop for efmas, currently.
 782            do mu=1,3
 783              do nu=1,3
 784                idir=3*(mu-1)+nu !xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9, (xyz,xyz)=(mu,nu)
 785                gh2c=zero; gs2c=zero
 786                call nonlop(choice_efmas,cpopt,cwaveprj,enlout_efmas,gs_ham,idir,(/eig_k/),mpi_enreg,&
 787                1,nnlout_efmas,paw_opt_efmas,signs,gs2c,tim_nonlop,cwavef,gh2c)
 788                do ispinor=1,dtset%nspinor
 789                  ii = 1+(ispinor-1)*npw_k
 790                  do icplx=1,2
 791                    gh2c(icplx,ii:ispinor*npw_k) = gh2c(icplx,ii:ispinor*npw_k) +  &
 792 &                   ddkinpw(1:npw_k,mu,nu)*cwavef(icplx,ii:ispinor*npw_k)
 793                  end do
 794                end do
 795                gh2c = gh2c - eig_k*gs2c
 796                ideg = efmasdeg(ikpt)%ideg(iband)
 797                do jband=efmasdeg(ikpt)%degs_bounds(1,ideg),efmasdeg(ikpt)%degs_bounds(2,ideg)
 798                  cg_left(:,1:npw_k*dtset%nspinor) = cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg)
 799                  dotprod=0
 800                  call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,cg_left,gh2c,mpi_enreg%me_g0,&
 801 &                 mpi_enreg%comm_spinorfft)
 802                  isub = iband-efmasdeg(ikpt)%degs_bounds(1,ideg)+1
 803                  jsub = jband-efmasdeg(ikpt)%degs_bounds(1,ideg)+1
 804                  efmasval(ideg,ikpt)%ch2c(mu,nu,jsub,isub)=cmplx(dotprod(1),dotprod(2),kind=dpc)
 805                end do
 806              end do
 807            end do
 808          end if
 809        end if
 810 
 811      end do ! End of loop on bands
 812 
 813      if(rfphon==1) then
 814        do iatom=1,natom
 815          ia=iatom;if (dyfr_nondiag==0) ia=1
 816          dyfrnl(1,1,1,iatom,ia)=dyfrnl(1,1,1,iatom,ia)+wtk_k*dyfrnlk(1,iatom)
 817          dyfrnl(1,2,2,iatom,ia)=dyfrnl(1,2,2,iatom,ia)+wtk_k*dyfrnlk(2,iatom)
 818          dyfrnl(1,3,3,iatom,ia)=dyfrnl(1,3,3,iatom,ia)+wtk_k*dyfrnlk(3,iatom)
 819          dyfrnl(1,2,3,iatom,ia)=dyfrnl(1,2,3,iatom,ia)+wtk_k*dyfrnlk(4,iatom)
 820          dyfrnl(1,1,3,iatom,ia)=dyfrnl(1,1,3,iatom,ia)+wtk_k*dyfrnlk(5,iatom)
 821          dyfrnl(1,1,2,iatom,ia)=dyfrnl(1,1,2,iatom,ia)+wtk_k*dyfrnlk(6,iatom)
 822        end do
 823      end if ! end if rfphon
 824      if(rfstrs/=0) then
 825        eltfrnl(:,:)=eltfrnl(:,:)+dtset%wtk(ikpt)*eltfrnlk(:,:)
 826      end if
 827      if(need_becfr) then
 828        becfrnl(:,:,:)=becfrnl(:,:,:)+dtset%wtk(ikpt)*becfrnlk(:,:,:)
 829      end if
 830      if(need_piezofr) then
 831        piezofrnl(:,:)=piezofrnl(:,:)+dtset%wtk(ikpt)*piezofrnlk(:,:)
 832      end if
 833 !    Increment indexes
 834      bdtot_index=bdtot_index+nband_k
 835      if (dtset%mkmem/=0) then
 836        ibg=ibg+nband_k*dtset%nspinor
 837        icg=icg+npw_k*dtset%nspinor*nband_k
 838        ikg=ikg+npw_k
 839      end if
 840 
 841      ABI_DEALLOCATE(ffnl)
 842      ABI_DEALLOCATE(kpg_k)
 843      ABI_DEALLOCATE(ph3d)
 844      ABI_DEALLOCATE(ylm_k)
 845      ABI_DEALLOCATE(ylmgr_k)
 846      ABI_DEALLOCATE(cwavef)
 847      ABI_DEALLOCATE(kg_k)
 848      if (need_becfr.or.need_piezofr) then
 849        ABI_DEALLOCATE(svectout)
 850      end if
 851      if (need_piezofr) then
 852        ABI_DEALLOCATE(ffnl_str)
 853      end if
 854      if (need_efmas) then
 855        ABI_DEALLOCATE(ddkinpw)
 856        ABI_DEALLOCATE(cg_left)
 857        ABI_DEALLOCATE(gh2c)
 858        ABI_DEALLOCATE(gs2c)
 859      end if
 860 
 861    end do ! End loops on isppol and ikpt
 862  end do
 863  if(rfphon==1) then
 864    ABI_DEALLOCATE(dyfrnlk)
 865    ABI_DEALLOCATE(enlout_phon)
 866  end if
 867  if(rfstrs/=0) then
 868    ABI_DEALLOCATE(eltfrnlk)
 869    ABI_DEALLOCATE(enlout_strs)
 870  end if
 871  if (need_becfr)  then
 872    ABI_DEALLOCATE(becfrnlk)
 873    ABI_DEALLOCATE(enlout_bec1)
 874    ABI_DEALLOCATE(enlout_bec2)
 875  end if
 876  if(need_piezofr)then
 877    ABI_DEALLOCATE(enlout_piez1)
 878    ABI_DEALLOCATE(enlout_piez2)
 879    ABI_DEALLOCATE(piezofrnlk)
 880  end if
 881  if(need_efmas) then
 882    ABI_DEALLOCATE(enlout_efmas)
 883  end if
 884  if (psps%usepaw==1) then
 885    if (need_becfr.or.need_piezofr)  then
 886      ABI_DEALLOCATE(becij)
 887    end if
 888    call pawcprj_free(cwaveprj)
 889  end if
 890  ABI_DATATYPE_DEALLOCATE(cwaveprj)
 891 
 892 !Fill in lower triangle of matrixes
 893  if (rfphon==1) then
 894    do iatom=1,natom
 895      ia=iatom;if (dyfr_nondiag==0) ia=1
 896      dyfrnl(1,3,2,iatom,ia)=dyfrnl(1,2,3,iatom,ia)
 897      dyfrnl(1,3,1,iatom,ia)=dyfrnl(1,1,3,iatom,ia)
 898      dyfrnl(1,2,1,iatom,ia)=dyfrnl(1,1,2,iatom,ia)
 899    end do
 900  end if
 901  if(rfstrs/=0)then
 902    do jj=2,6
 903      do ii=1,jj-1
 904        eltfrnl(jj,ii)=eltfrnl(ii,jj)
 905      end do
 906    end do
 907  end if
 908 
 909 !Parallel case: accumulate (n,k) contributions
 910  if (xmpi_paral==1) then
 911    call timab(48,1,tsec)
 912 !  Accumulate dyfrnl
 913    if(rfphon==1)then
 914      call xmpi_sum(dyfrnl,spaceworld,ierr)
 915    end if
 916 !  Accumulate eltfrnl.
 917    if(rfstrs/=0)then
 918      call xmpi_sum(eltfrnl,spaceworld,ierr)
 919    end if
 920 !  Accumulate becfrnl
 921    if (need_becfr) then
 922      call xmpi_sum(becfrnl,spaceworld,ierr)
 923    end if
 924 !  Accumulate piezofrnl
 925    if (need_piezofr) then
 926      call xmpi_sum(piezofrnl,spaceworld,ierr)
 927    end if
 928 
 929 !  PAW: accumulate gradients of rhoij
 930    if (psps%usepaw==1) then
 931      ABI_ALLOCATE(dimlmn,(natom))
 932      dimlmn(1:natom)=pawrhoij_tot(1:natom)%cplex*pawrhoij_tot(1:natom)%lmn2_size
 933      bufdim=ncpgr*sum(dimlmn)*nsploop
 934      ABI_ALLOCATE(mpibuf,(bufdim))
 935      ii=0;mpibuf=zero
 936      do iatom=1,natom
 937        do isppol=1,nsploop
 938          do mu=1,ncpgr
 939            mpibuf(ii+1:ii+dimlmn(iatom))=pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)
 940            ii=ii+dimlmn(iatom)
 941          end do
 942        end do
 943      end do
 944      call xmpi_sum(mpibuf,spaceworld,ierr)
 945      ii=0
 946      do iatom=1,natom
 947        do isppol=1,nsploop
 948          do mu=1,ncpgr
 949            pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)=mpibuf(ii+1:ii+dimlmn(iatom))
 950            ii=ii+dimlmn(iatom)
 951          end do
 952        end do
 953      end do
 954      ABI_DEALLOCATE(mpibuf)
 955      ABI_DEALLOCATE(dimlmn)
 956    end if
 957    call timab(48,2,tsec)
 958  end if
 959 
 960 !====== PAW: Additional steps
 961  if (psps%usepaw==1) then
 962 
 963 !  Symmetrize rhoij gradients and transfer to cartesian (reciprocal space) coord.
 964 !  This symetrization is necessary in the antiferromagnetic case...
 965    if (rfphon==1.and.rfstrs==0) then
 966      option_rhoij=2;option=0
 967      call symrhoij(pawrhoij_tot,pawrhoij_tot,option_rhoij,gprimd,indsym,0,natom,dtset%nsym,&
 968 &     psps%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,&
 969 &     comm_atom=my_comm_atom, mpi_atmtab=my_atmtab)
 970    else if (rfphon==1.and.rfstrs==1) then
 971      option_rhoij=23;option=0
 972      call symrhoij(pawrhoij_tot,pawrhoij_tot,option_rhoij,gprimd,indsym,0,natom,dtset%nsym,&
 973 &     psps%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,&
 974 &     comm_atom=my_comm_atom, mpi_atmtab=my_atmtab)
 975    end if
 976 
 977 !  Translate coordinates
 978    do iatom=1,natom
 979      cplx=pawrhoij_tot(iatom)%cplex
 980      do isppol=1,nsploop
 981        do klmn=1,pawrhoij_tot(iatom)%lmn2_size
 982          do ii=1,cplx
 983            if(rfphon==1.or.rfstrs/=0)then
 984              grhoij(1:3)=pawrhoij_tot(iatom)%grhoij(shift_rhoij+1:shift_rhoij+3,cplx*(klmn-1)+ii,isppol)
 985              do mu=1,3
 986                pawrhoij_tot(iatom)%grhoij(shift_rhoij+mu,cplx*(klmn-1)+ii,isppol)=gprimd(mu,1)*grhoij(1)&
 987 &               +gprimd(mu,2)*grhoij(2)+gprimd(mu,3)*grhoij(3)
 988              end do
 989            end if
 990            if(rfstrs/=0)then
 991              call strconv(pawrhoij_tot(iatom)%grhoij(1:6,cplx*(klmn-1)+ii,isppol),gprimd,&
 992 &             pawrhoij_tot(iatom)%grhoij(1:6,cplx*(klmn-1)+ii,isppol))
 993            end if
 994          end do
 995        end do
 996      end do
 997    end do
 998 
 999 !  In case of elastic tensor computation, add diagonal contribution:
1000 !     -delta_{alphabeta} rhoi_{ij} to drhoij/d_eps
1001    if(rfstrs/=0)then
1002      do iatom=1,natom
1003        cplx=pawrhoij_tot(iatom)%cplex
1004        do isppol=1,nsploop
1005          do nu=1,pawrhoij_tot(iatom)%nrhoijsel
1006            klmn=pawrhoij_tot(iatom)%rhoijselect(nu)
1007            do ii=1,cplx
1008              pawrhoij_tot(iatom)%grhoij(1:3,cplx*(klmn-1)+ii,isppol)= &
1009 &             pawrhoij_tot(iatom)%grhoij(1:3,cplx*(klmn-1)+ii,isppol)&
1010 &             -pawrhoij_tot(iatom)%rhoijp(cplx*(nu-1)+ii,isppol)
1011            end do
1012          end do
1013        end do
1014      end do
1015    end if
1016 
1017 !  Add gradients due to Dij derivatives to dynamical matrix/stress tensor
1018    dimnhat=0;optgr=0;optgr2=0;optstr=0;optstr2=0
1019    if (rfphon==1) optgr2=1
1020    if (rfstrs/=0) optstr2=1
1021    ABI_ALLOCATE(nhat_dum,(1,0))
1022    call pawgrnl(gs_ham%atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,dummy,gsqcut,mgfftf,my_natom,natom,&
1023 &   gs_ham%nattyp,nfftf,ngfftf,nhat_dum,dummy,dtset%nspden,dtset%nsym,psps%ntypat,optgr,optgr2,optstr,optstr2,&
1024 &   pawang,pawfgrtab,pawrhoij_tot,pawtab,ph1df,psps,dtset%qptn,rprimd,symrec,dtset%typat,ucvol,vtrial,vxc,xred,&
1025 &   mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
1026    ABI_DEALLOCATE(nhat_dum)
1027  end if !PAW
1028 
1029 !The indexing array atindx is used to reestablish the correct order of atoms
1030  if (rfstrs/=0)then
1031    ABI_ALLOCATE(elt_work,(6+3*natom,6))
1032    elt_work(1:6,1:6)=eltfrnl(1:6,1:6)
1033    do ia=1,natom
1034      ielt=7+3*(ia-1)
1035      ieltx=7+3*(gs_ham%atindx(ia)-1)
1036      elt_work(ielt:ielt+2,1:6)=eltfrnl(ieltx:ieltx+2,1:6)
1037    end do
1038    eltfrnl(:,:)=elt_work(:,:)
1039    ABI_DEALLOCATE(elt_work)
1040  end if
1041 
1042 !Born Effective Charges and PAW:
1043 !1-Re-order atoms -- 2-Add diagonal contribution from rhoij
1044 !3-Multiply by -1 because that the effective charges
1045  !  are minus the second derivatives of the energy
1046  if (need_becfr) then
1047    ABI_ALLOCATE(becfrnl_tmp,(3,natom,3))
1048    becfrnl_tmp=-becfrnl
1049    do ia=1,natom         ! Atom (sorted by type)
1050      iatom=gs_ham%atindx1(ia)   ! Atom (not sorted)
1051      itypat=dtset%typat(iatom)
1052      do ii=1,3           ! Direction of electric field
1053        do jj=1,3         ! Direction of atom
1054          becfrnl(jj,iatom,ii)=becfrnl_tmp(jj,ia,ii)
1055        end do
1056      end do
1057    end do
1058    ABI_DEALLOCATE(becfrnl_tmp)
1059  end if
1060 
1061 !Piezoelectric Tensor
1062 !-Multiply by -1 because that the piezoelectric tensor
1063 !  are minus the second derivatives of the energy
1064  if (need_piezofr) then
1065    piezofrnl=-piezofrnl
1066  end if
1067 
1068 !Close the ddk files
1069  do ii=1,3
1070    call wfk_close(ddkfiles(ii))
1071  end do
1072 
1073 !Release now useless memory
1074  if (psps%usepaw==1) then
1075    do iatom=1,natom
1076      ABI_DEALLOCATE(pawrhoij_tot(iatom)%grhoij)
1077      pawrhoij_tot(iatom)%ngrhoij=0
1078    end do
1079    if (paral_atom) then
1080      call pawrhoij_free(pawrhoij_tot)
1081      ABI_DATATYPE_DEALLOCATE(pawrhoij_tot)
1082    end if
1083  end if
1084  call destroy_hamiltonian(gs_ham)
1085  call timab(159,2,tsec)
1086 
1087  write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivative done',ch10
1088  call wrtout(std_out,msg,'COLL')
1089 
1090  DBG_EXIT("COLL")
1091 
1092 end subroutine d2frnl

ABINIT/m_d2frnl [ Modules ]

[ Top ] [ Modules ]

NAME

  m_d2frnl

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GM, AR, MB, 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 .

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_d2frnl
28 
29  implicit none
30 
31  private