TABLE OF CONTENTS


ABINIT/berryphase_new [ Functions ]

[ Top ] [ Functions ]

NAME

 berryphase_new

FUNCTION

 This routine computes the Berry Phase polarization
  and the finite difference expression of the ddk.
  See for example Na Sai et al., PRB 66, 104108 (2002) [[cite:Sai2002]]

INPUTS

 atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
 cg(2,mcg)=planewave coefficients of wavefunctions
 cprj(natom,mcprj*usecrpj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
 dtfil <type(datafiles_type)>=variables related to files
 dtset <type(dataset_type)>=all input variables in this dataset
 psps <type(pseudopotential_type)>=variables related to pseudopotentials
 gprimd(3,3)=reciprocal space dimensional primitive translations
 hdr <type(hdr_type)>=the header of wf, den and pot files
 indlmn(6,lmnmax,ntypat)
   array giving l,m,n,lm,ln,spin for i=ln  (if useylm=0)
                                  or i=lmn (if useylm=1)
 kg(3,mpw*mkmem)=reduced planewave coordinates
 lmnmax  If useylm=0, max number of (l,m,n) comp. over all type of psps (lnproj)
         If useylm=1, max number of (l,n)   comp. over all type of psps (lmnproj)
 mband=maximum number of bands
 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
 mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
 mkmem=number of k points treated by this node
 mpi_enreg=information about MPI parallelization
 mpw=maximum dimensioned size of npw
 my_natom=number of atoms treated by current processor
 natom=number of atoms in cell
 nkpt=number of k points
 npwarr(nkpt)=number of planewaves in basis at this k point
 nsppol=1 for unpolarized, 2 for spin-polarized
 ntypat=number of types of atoms in unit cell
 nkpt=number of k-points
 calc_pol_ddk = 1: compute Berryphase polarization
                2: compute finite difference expression of the ddk
                3: compute polarization & ddk
 pawrhoij(natom*usepaw) <type(pawrhoij_type)> atomic occupancies
 pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
 pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
 pwind_alloc = first dimension of pwind
 pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
 rprimd(3,3)=dimensional real space primitive translations (bohr)
 typat(natom)=type integer for each atom in cell
 ucvol=unit cell volume in bohr**3.
 unit_out= unit for output of the results (usually the .out file of ABINIT)
   The option unit_out = 0 is allowed. In this case, no information is written
   to the output file but only to the log file.
 usecprj=1 if cprj datastructure has been allocated
 usepaw= 1: use paw framework. 0:do not use paw.
 xred(3,natom)=reduced atomic coordinates
 zion(ntypat)=valence charge of each type of atom

OUTPUT

 ptot(3) = total polarization including correction for jumps
 red_ptot(3) = total polarization including correction for jumps reduced units
 pel(3) = reduced coordinates of the electronic polarization (a. u.)
 pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates (already contained in pel)
 pion(3)= reduced coordinates of the ionic polarization (a. u.)

SIDE EFFECTS

 Input/Output
 dtefield <type(efield_type)> = variables related to Berry phase
       and electric field calculations (see initberry.f).
       In case berryopt = 4, the overlap matrices computed
       in this routine are stored in dtefield%smat in order
       to be used in the electric field calculation.

TODO

  - Use the analytical relation between the overlap matrices
    S(k,k+dk) and S(k+dk,k) to avoid to recompute them
    when ifor = 2.

NOTES

 - pel and pion do not take into account the factor 1/ucvol
 - In case of a ddk calculation, the eigenvalues are not computed.
 - The ddk computed by this routine should not be used to
   compute the electronic dielectric tensor.

PARENTS

      elpolariz,update_e_field_vars

CHILDREN

      appdig,dsdr_k_paw,outwf,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawcprj_get,pawcprj_getdim,pawcprj_mpi_allgather,pawcprj_mpi_recv
      pawcprj_mpi_send,pawcprj_put,pawcprj_symkn,pawpolev,polcart,rhophi
      smatrix,smatrix_k_paw,wrtout,xmpi_allgather,xmpi_bcast,xmpi_recv
      xmpi_send,xmpi_sum,xred2xcart

SOURCE

 171 subroutine berryphase_new(atindx1,cg,cprj,dtefield,dtfil,dtset,psps,&
 172 &  gprimd,hdr,indlmn,kg,lmnmax,mband,mcg,mcprj,&
 173 &  mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,nsppol,ntypat,&
 174 &  nkpt,calc_pol_ddk,pawrhoij,pawtab,pel,pelev,pion,ptot,red_ptot,pwind,&  !!REC
 175 &  pwind_alloc,pwnsfac,&
 176 &  rprimd,typat,ucvol,unit_out,usecprj,usepaw,xred,zion)
 177 
 178 
 179 !This section has been created automatically by the script Abilint (TD).
 180 !Do not modify the following lines by hand.
 181 #undef ABI_FUNC
 182 #define ABI_FUNC 'berryphase_new'
 183 !End of the abilint section
 184 
 185  implicit none
 186 
 187 !Arguments ------------------------------------
 188  integer, intent(in) :: lmnmax,mband,mcg,mcprj,mkmem,mpw,my_natom,natom,nkpt
 189  integer, intent(in) :: nsppol,ntypat,calc_pol_ddk
 190  integer, intent(in) :: pwind_alloc,unit_out,usecprj,usepaw
 191  real(dp), intent(in) :: ucvol
 192  type(MPI_type), intent(in) :: mpi_enreg
 193  type(datafiles_type), intent(in) :: dtfil
 194  type(dataset_type), intent(in) :: dtset
 195  type(pseudopotential_type),intent(in) :: psps
 196  type(efield_type), intent(inout) :: dtefield
 197  type(hdr_type), intent(inout) :: hdr
 198 !arrays
 199  integer, intent(in) :: atindx1(natom),indlmn(6,lmnmax,ntypat),kg(3,mpw*mkmem)
 200  integer, intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
 201  integer, intent(in) :: typat(natom)
 202  real(dp), intent(in) :: cg(2,mcg),gprimd(3,3)
 203  real(dp), intent(in) :: pwnsfac(2,pwind_alloc)
 204  real(dp), intent(in) :: rprimd(3,3),zion(ntypat)
 205  real(dp), intent(inout) :: xred(3,natom)
 206  real(dp), intent(out) :: pel(3),pelev(3),pion(3)
 207  real(dp), intent(out) :: ptot(3),red_ptot(3) !!REC
 208  type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw)
 209  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
 210  type(pawcprj_type),intent(in) ::  cprj(natom,mcprj*usecprj)
 211 
 212 !Local variables -------------------------
 213  integer :: count,count1,dest,fdir,unt
 214  integer :: iatom,iband,icg,icg1,idir,idum,ikpt1i_sp
 215  integer :: ierr,ifor,ikg,ikpt,ikpt1,ikpt_loc!,ikpt2,ikpt2i,npw_k2, itrs
 216  integer :: icp1, icp2,icpgr_offset,iproc
 217 ! integer :: ii ! appears commented out below in a debug section
 218  integer :: inibz,ikpt1i
 219  integer :: isppol,istr,itypat,jband,jkpt,jkstr,jsppol
 220  integer :: det_inv_smat, det_smat, inv_smat
 221  integer :: maxbd,mcg1_k
 222  integer :: minbd,my_nspinor,nband_k,ncpgr,nfor,npw_k1,ntotcp,n2dim,nproc,pertcase
 223  integer :: response,shiftbd,source,spaceComm,tag
 224  integer :: jj,jstr,kk,ineigh_str
 225  integer :: istep,jstep,kpt_mark(dtefield%fnkpt),nkstr,nstr,iunmark,berrystep
 226  integer :: jkpt2, jkpt2i
 227  real(dp) :: det_mod,dkinv,dphase,dtm_real,dtm_imag,fac,gmod,phase0
 228  real(dp) :: pol,polbtot,polion,politot,poltot,rho
 229  logical :: calc_epaw3_force,calc_epaw3_stress,efield_flag
 230  integer :: polflag, ddkflag
 231 
 232 !!REC start
 233  integer :: jump
 234  real(dp),save ::pol0(3)
 235  logical, save :: first=.true.
 236  logical :: lexist
 237 !!REC end
 238  real(dp) :: dphase_new,dphase_init
 239  character(len=fnlen) :: fiwf1o
 240  character(len=500) :: message
 241  type(wvl_wf_type) :: wfs
 242  type(wvl_internal_type) :: wvl
 243  integer,allocatable :: dimlmn(:),ikpt1_recv(:), sflag_k(:)!,pwind_k(:)
 244  integer,allocatable :: ikpt3(:), ikpt3i(:), sflag_k_mult(:,:),nattyp_dum(:),npw_k3(:)
 245  integer,allocatable :: idxkstr_mult(:,:), pwind_k_mult(:,:),itrs_mult(:)
 246  real(dp) :: det_average(2),dk(3),dtm_k(2),gpard(3),pel_cart(3),pion_cart(3)
 247  real(dp) :: polb(nsppol),ptot_cart(3),rel_string(2),xcart(3,natom)
 248  real(dp) :: delta_str(2),dist_,dstr(2)
 249  real(dp),allocatable :: buffer(:,:),buffer1(:),buffer2(:)
 250  real(dp),allocatable :: cg1(:,:),cg1_k(:,:),cgq(:,:)
 251  real(dp),allocatable :: det_string(:,:),dsdr(:,:,:,:,:),dsdr_sum(:,:,:),dsds_sum(:,:,:)
 252  real(dp),allocatable :: eig_dum(:),epawf3_str(:,:,:),epaws3_str(:,:,:)
 253  real(dp),allocatable :: occ_dum(:),polberry(:),resid(:),pwnsfac_k(:,:)
 254  real(dp),allocatable :: smat_inv(:,:,:),smat_k(:,:,:),smat_k_paw(:,:,:)
 255  real(dp),allocatable :: str_flag(:)
 256 ! real(dp),allocatable :: dist_str(:,:),det_string_test(:,:)
 257  real(dp),allocatable :: dtm_mult(:,:,:), coef(:,:), polb_mult(:,:)
 258  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:),cprj_buf(:,:),cprj_gat(:,:)
 259  type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:)
 260 
 261 ! integer :: bband,bbs,bra_start,bra_end,ilmn,ipw,ispinor,jlmn,kband,kbs,ket_start,ket_end,klmn,npw_k
 262 ! integer :: nspinor,spnipw,spnshft
 263 ! real(dp) :: err_ovlp,mag_ovlp,max_err_ovlp, ovlp_r, ovlp_i, paw_r, paw_i
 264 ! real(dp) :: tot_r, tot_i
 265 ! real(dp),allocatable :: bra(:,:),ket(:,:)
 266 ! complex(dpc) :: cpb,cpk,cterm
 267 
 268 !BEGIN TF_CHANGES
 269  integer :: me
 270 !END TF_CHANGES
 271 
 272 !no_abirules
 273 
 274 ! ***********************************************************************
 275 
 276 !DEBUG
 277 !write(std_out,*)' berryphase_new : enter'
 278 !do ii=1,pwind_alloc
 279 !write(std_out,*)ii,pwnsfac(:,ii)
 280 !end do
 281 !stop
 282 !ENDDEBUG
 283 
 284 !Init MPI
 285  spaceComm=mpi_enreg%comm_cell
 286  nproc=xmpi_comm_size(spaceComm)
 287  me=mpi_enreg%me_kpt
 288 
 289  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 290 
 291  polflag = 1
 292  ddkflag = 1
 293  if (calc_pol_ddk == 1) then
 294    ddkflag = 0
 295  else if (calc_pol_ddk == 2) then
 296    polflag = 0
 297  end if
 298 
 299 !allocate(pwind_k(mpw))
 300  ABI_ALLOCATE(pwnsfac_k,(4,mpw))
 301  ABI_ALLOCATE(sflag_k,(dtefield%mband_occ))
 302 !pwind_k(:) = 0
 303  pwnsfac_k(1,:) = 1.0_dp ! bra real
 304  pwnsfac_k(2,:) = 0.0_dp ! bra imag
 305  pwnsfac_k(3,:) = 1.0_dp ! ket real
 306  pwnsfac_k(4,:) = 0.0_dp ! ket imag
 307 
 308  if (maxval(dtset%istwfk(:)) /= 1) then
 309    write(message, '(a,a,a)' )&
 310 &   'This routine does not work yet with istwfk /= 1.',ch10,&
 311 &   'This should have been tested previously ...'
 312    MSG_BUG(message)
 313  end if
 314 
 315  if (usepaw == 1 .and. usecprj /= 1) then
 316    message = ' PAW calculation but cprj datastructure has not been allocated !'
 317    MSG_BUG(message)
 318  end if
 319 
 320 ! useful flags for various efield possibilities
 321  efield_flag = ( dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. &
 322 & dtset%berryopt ==14 .or. dtset%berryopt ==16 .or. dtset%berryopt ==17 )
 323  calc_epaw3_force = ( efield_flag .and. dtset%optforces /= 0 .and. usepaw == 1 )
 324  calc_epaw3_stress = ( efield_flag .and. dtset%optstress /= 0  .and. usepaw == 1 )
 325 
 326  mcg1_k = mpw*mband
 327  shiftbd = 1
 328  if (ddkflag==1) then
 329    ABI_ALLOCATE(cg1,(2,mcg))
 330    ABI_ALLOCATE(eig_dum,(2*mband*mband*nkpt*nsppol))
 331    ABI_ALLOCATE(occ_dum,(mband*nkpt*nsppol))
 332    eig_dum(:) = zero
 333    occ_dum(:) = dtefield%sdeg
 334  end if
 335 
 336 !initialize variable tied to multiple step computation
 337  berrystep=dtset%berrystep
 338  ABI_ALLOCATE(ikpt3,(berrystep))
 339  ABI_ALLOCATE(ikpt3i,(berrystep))
 340  ABI_ALLOCATE(sflag_k_mult,(dtefield%mband_occ,berrystep))
 341  ABI_ALLOCATE(npw_k3,(berrystep))
 342  ABI_ALLOCATE(pwind_k_mult,(mpw,berrystep))
 343  ABI_ALLOCATE(itrs_mult,(berrystep))
 344  ABI_ALLOCATE(coef,(berrystep,berrystep))
 345  ABI_ALLOCATE(polb_mult,(nsppol,berrystep))
 346 !coefficient for berryphase computation
 347  coef(:,:) = 0.0_dp
 348  do jstep = 1, berrystep
 349    coef(jstep,1) = 1.d0/real(jstep*jstep,dp)
 350    if(jstep/=1)coef(jstep,1)=coef(jstep,1)/real(1-jstep*jstep,dp)
 351  end do
 352  do istep = 2, berrystep
 353    do jstep = 1, berrystep
 354      coef(jstep, istep) = real(istep*istep,dp)*coef(jstep,istep-1)
 355      if(jstep /= istep)coef(jstep, istep)=coef(jstep,istep)/real(istep*istep-jstep*jstep,dp)
 356    end do
 357  end do
 358 !the berryphase using the strings of steps dk, 2*dk, ..., istep*dk is :
 359 !coef(1,istep)*berryphase(dk) + coef(2,istep)*berryphase(2*dk) + ... + coef(istep,istep)*berryphase(istep*dk)
 360 !DEBUG
 361 !write(std_out,*)'coef, sum coef'
 362 !do istep=1,step
 363 !write(std_out,*)coef(:,istep), sum(coef(1:istep,istep))
 364 !end do
 365 !ENDDEBUG
 366 
 367 !allocate(dtm(2,dtefield%fnkpt*nsppol))
 368  ABI_ALLOCATE(dtm_mult,(2,dtefield%fnkpt*nsppol,berrystep))
 369  ABI_ALLOCATE(cg1_k,(2,mcg1_k))
 370 
 371  if (usepaw == 1) then ! cprj allocation
 372    ncpgr = cprj(1,1)%ncpgr
 373    if ( calc_epaw3_force ) then
 374      ABI_ALLOCATE(dsdr_sum,(natom,3,dtefield%fnkpt*nsppol))
 375      ABI_ALLOCATE(epawf3_str,(natom,3,3))
 376    end if
 377    if ( calc_epaw3_stress ) then
 378      ABI_ALLOCATE(dsds_sum,(natom,6,dtefield%fnkpt*nsppol))
 379      ABI_ALLOCATE(epaws3_str,(natom,3,6))
 380    end if
 381    ABI_ALLOCATE(dimlmn,(natom))
 382    call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R')
 383    ABI_DATATYPE_ALLOCATE(cprj_k,(natom,dtefield%nspinor*mband))
 384    ABI_DATATYPE_ALLOCATE(cprj_kb,(natom,dtefield%nspinor*mband))
 385    ABI_DATATYPE_ALLOCATE(cprj_gat,(natom,nproc*dtefield%nspinor*mband))
 386    call pawcprj_alloc(cprj_k,ncpgr,dimlmn)
 387    call pawcprj_alloc(cprj_kb,ncpgr,dimlmn)
 388    call pawcprj_alloc(cprj_gat,ncpgr,dimlmn)
 389    if (dtset%kptopt /= 3) then
 390      ABI_DATATYPE_ALLOCATE(cprj_ikn,(natom,dtefield%nspinor*mband))
 391      ABI_DATATYPE_ALLOCATE(cprj_fkn,(natom,dtefield%nspinor*mband))
 392      call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn)
 393      call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn)
 394    end if
 395 
 396    n2dim = dtefield%nspinor*mband
 397    ntotcp = n2dim*SUM(dimlmn(:))
 398    if (nproc>1) then
 399      ABI_DATATYPE_ALLOCATE(cprj_buf,(natom,dtefield%nspinor*mband))
 400      call pawcprj_alloc(cprj_buf,ncpgr,dimlmn)
 401    end if
 402 
 403    if ( efield_flag ) then
 404      write(message,'(2a,i5,2a)')ch10,&
 405 &     ' nkpt = ',nkpt,ch10,' copy cprj to dtefield%cprj '
 406      call wrtout(std_out,message,'COLL')
 407 
 408      do isppol = 1, nsppol
 409 
 410        ikpt_loc = 0
 411        ikpt1 = 0
 412        do while (ikpt_loc < mkmem)
 413 
 414          if (ikpt_loc < mkmem) ikpt1 = ikpt1 + 1
 415          if ((ikpt1 > nkpt).and.(ikpt_loc < mkmem)) exit
 416          nband_k = dtset%nband(ikpt1)
 417 
 418          if ( (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt1,1,nband_k,isppol,me)) .and. &
 419 &         (ikpt_loc <= mkmem) ) cycle
 420 
 421          ikpt_loc = ikpt_loc + 1
 422 
 423          ABI_ALLOCATE(ikpt1_recv,(nproc))
 424          call xmpi_allgather(ikpt1,ikpt1_recv,spaceComm,ierr)
 425          call pawcprj_get(atindx1,cprj_k,cprj,natom,1,(ikpt_loc-1)*nband_k*my_nspinor,ikpt1,0,isppol,mband,&
 426 &         mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,0,&
 427 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 428          call pawcprj_mpi_allgather(cprj_k,cprj_gat,natom,n2dim,dimlmn,ncpgr,nproc,spaceComm,ierr,rank_ordered=.true.)
 429          do iproc = 1, nproc
 430            icp2=nband_k*(iproc-1)*my_nspinor
 431            call pawcprj_get(atindx1,cprj_k,cprj_gat,natom,1,icp2,ikpt1,0,isppol,mband,&
 432 &           nproc,natom,nband_k,nband_k,my_nspinor,1,0,&
 433 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 434            icp1 = nband_k*(ikpt1_recv(iproc)-1)*my_nspinor
 435            call pawcprj_put(atindx1,cprj_k,dtefield%cprj,natom,1,icp1,ikpt1,0,isppol,&
 436 &           mband,dtefield%fnkpt,natom,nband_k,nband_k,dimlmn,my_nspinor,nsppol,0,&
 437 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 438          end do
 439          ABI_DEALLOCATE(ikpt1_recv)
 440 
 441        end do ! close loop over k-points
 442 
 443      end do ! end loop over nsppol
 444 
 445    end if ! end check on efield
 446  end if
 447 
 448 !!=======================================
 449 !! code to test orthonormality of cg_k
 450 !!=======================================
 451 !
 452 !ikpt = 3
 453 !npw_k = npwarr(ikpt)
 454 !isppol = 1
 455 !nband_k = dtefield%mband_occ
 456 !ABI_ALLOCATE(bra,(2,npw_k*my_nspinor))
 457 !ABI_ALLOCATE(ket,(2,npw_k*my_nspinor))
 458 !max_err_ovlp=0.0
 459 !call pawcprj_get(atindx1,cprj_k,cprj,natom,1,dtefield%cprjindex(ikpt,isppol),ikpt,0,isppol,mband,&
 460 !&         mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,0)
 461 !do bband = 1, nband_k
 462 !bra_start = dtefield%cgindex(ikpt,nsppol)+1+(bband-1)*npw_k*my_nspinor
 463 !bra_end = bra_start + npw_k*my_nspinor - 1
 464 !bra(1:2,1:npw_k*my_nspinor) = cg(1:2,bra_start:bra_end)
 465 !do kband = 1, nband_k
 466 !ket_start = dtefield%cgindex(ikpt,nsppol)+1+(kband-1)*npw_k*my_nspinor
 467 !ket_end = ket_start + npw_k*my_nspinor - 1
 468 !ket(1:2,1:npw_k*my_nspinor) = cg(1:2,ket_start:ket_end)
 469 !
 470 !tot_r = 0.0; tot_i = 0.0
 471 !do ispinor = 1, my_nspinor
 472 !ovlp_r = 0.0; ovlp_i = 0.0
 473 !spnshft = (ispinor-1)*npw_k
 474 !do ipw = 1, npw_k
 475 !spnipw = ipw + spnshft
 476 !ovlp_r = ovlp_r + bra(1,spnipw)*ket(1,spnipw)+bra(2,spnipw)*ket(2,spnipw)
 477 !ovlp_i = ovlp_i - bra(2,spnipw)*ket(1,spnipw)+bra(1,spnipw)*ket(2,spnipw)
 478 !end do ! end loop over ipw
 479 !paw_r = 0.0; paw_i = 0.0
 480 !do iatom = 1, natom
 481 !itypat = typat(iatom)
 482 !do ilmn = 1, dtefield%lmn_size(itypat)
 483 !do jlmn = 1, dtefield%lmn_size(itypat)
 484 !klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
 485 !bbs = my_nspinor*(bband-1)+ispinor
 486 !kbs = my_nspinor*(kband-1)+ispinor
 487 !cpb=cmplx(cprj_k(iatom,bbs)%cp(1,ilmn),cprj_k(iatom,bbs)%cp(2,ilmn))
 488 !cpk=cmplx(cprj_k(iatom,kbs)%cp(1,jlmn),cprj_k(iatom,kbs)%cp(2,jlmn))
 489 !cterm = conjg(cpb)*pawtab(itypat)%sij(klmn)*cpk
 490 !paw_r = paw_r + real(cterm)
 491 !paw_i = paw_i + aimag(cterm)
 492 !end do ! end loop over jlmn
 493 !end do ! end loop over ilmn
 494 !end do ! end loop over iatom
 495 !tot_r = tot_r + ovlp_r + paw_r
 496 !tot_i = tot_i + ovlp_i + paw_i
 497 !end do ! end loop over ispinor
 498 !
 499 !!     write(std_out,'(a,2i4,2es16.8)')' JWZ Debug: berryphase_new bband kband ovlp : ',&
 500 !!&           bband,kband,tot_r,tot_i
 501 !mag_ovlp =  tot_r*tot_r + tot_i*tot_i
 502 !if(bband==kband) then
 503 !err_ovlp=abs(mag_ovlp-1.0)
 504 !else
 505 !err_ovlp=abs(mag_ovlp)
 506 !end if
 507 !max_err_ovlp=MAX(max_err_ovlp,err_ovlp)
 508 !end do ! end loop over kband
 509 !end do ! end loop over bband
 510 !write(std_out,'(a,i4,es16.8)')' JWZ Debug: berrphase_new ikpt ovlp err : ',&
 511 !&           ikpt,max_err_ovlp
 512 !ABI_DEALLOCATE(bra)
 513 !ABI_DEALLOCATE(ket)
 514 !
 515 !!=========================================
 516 !! end code to test orthonormality of cg_k
 517 !!=========================================
 518 
 519  pel(:) = zero ; pelev(:)=zero ; pion(:) = zero ; ptot(:)=zero ; red_ptot(:)=zero
 520 
 521  minbd = 1   ;  maxbd = dtefield%mband_occ
 522 
 523  if(calc_epaw3_force) dtefield%epawf3(:,:,:) = zero
 524  if(calc_epaw3_stress) dtefield%epaws3(:,:,:) = zero
 525 
 526  do idir = 1, 3
 527 
 528 !  dtm(:,:) = zero
 529    dtm_mult(:,:,:) = zero
 530    if (calc_epaw3_force) dsdr_sum(:,:,:) = zero
 531    if (calc_epaw3_stress) dsds_sum(:,:,:) = zero
 532 
 533    if (dtset%rfdir(idir) /= 1) cycle
 534 
 535    if (abs(dtefield%efield_dot(idir)) < tol12) dtefield%sflag(:,:,:,idir) = 0
 536 
 537 ! calculate vector steps in k space
 538    dk(:) = dtefield%dkvecs(:,idir)
 539    gpard(:) = dk(1)*gprimd(:,1) + dk(2)*gprimd(:,2) + dk(3)*gprimd(:,3)
 540    gmod = sqrt(dot_product(gpard,gpard))
 541 
 542    write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
 543 &   ' Computing the polarization (Berry phase) for reciprocal vector:',ch10,&
 544 &   dk(:),' (in reduced coordinates)',ch10,&
 545 &   gpard(1:3),' (in cartesian coordinates - atomic units)'
 546    call wrtout(std_out,message,'COLL')
 547    if (unit_out /= 0) then
 548      call wrtout(unit_out,message,'COLL')
 549    end if
 550 
 551    write(message,'(a,i5,a,a,i5)')&
 552 &   ' Number of strings: ',dtefield%nstr(idir),ch10,&
 553 &   ' Number of k points in string:', dtefield%nkstr(idir)
 554    call wrtout(std_out,message,'COLL')
 555    if (unit_out /= 0) then
 556      call wrtout(unit_out,message,'COLL')
 557    end if
 558 
 559 !  Check whether the polarization or the ddk must be computed
 560 
 561 !  nfor = 1 : to compute P, I only need the WF at k + dk
 562 !  nfor = 2 : to compute the ddk, I need the WF at k + dk and k - dk
 563 !  dkinv    : +-1/2dk
 564 
 565 
 566 !  default for polarization
 567    nfor = 1
 568    if (ddkflag == 1) then
 569      nfor = 2
 570    end if
 571 
 572    if (ddkflag == 1) then
 573 
 574      cg1(:,:) = zero
 575      dkinv = one/(two*dk(idir))
 576 
 577      write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
 578 &     ' Computing the ddk (Berry phase) for reciprocal vector:',ch10,&
 579 &     dk(:),' (in reduced coordinates)',ch10,&
 580 &     gpard(1:3),' (in cartesian coordinates - atomic units)'
 581      call wrtout(std_out,message,'COLL')
 582      if (unit_out /= 0) then
 583        call wrtout(unit_out,message,'COLL')
 584      end if
 585 
 586    end if
 587 
 588 ! From smatrix routine: det_inv_smat = type of calculation
 589 !        1 : compute inverse of the overlap matrix
 590 !       10 : compute determinant of the overlap matrix
 591 !       11 : compute determinant and inverse of the overlap matrix
 592    inv_smat = 0
 593    det_smat = 0
 594 
 595 !  for ddk need inverse matrix
 596    if (ddkflag == 1) then
 597      inv_smat = 1
 598    end if
 599 
 600 !   if polarization is requested need smat determinant as well
 601    if (polflag == 1) then
 602      det_smat = 1
 603    end if
 604 
 605 ! electric fields with PAW also needs S_inverse for forces and stresses, even just for polarization
 606    if (calc_epaw3_force .or. calc_epaw3_stress) then
 607      inv_smat = 1
 608    end if
 609 
 610    det_inv_smat = 10*det_smat  + inv_smat
 611 
 612 !--------------------------------------------------------------------
 613 !  for each dk we require, calculate the smatrix, derivatives etc...
 614 !--------------------------------------------------------------------
 615    do ifor = 1, nfor
 616 
 617      if (ifor == 2) then
 618        dk(:) = -1_dp*dk(:)
 619 !      only the inverse of the overlap matrix is required on second pass, speeds things up a bit
 620        det_inv_smat = 1
 621        dkinv = -1_dp*dkinv
 622      end if
 623 
 624 
 625 !    Compute the determinant and/or the inverse of the overlap matrix
 626 !    for each pair of k-points < u_nk | u_nk+dk >
 627 
 628      icg = 0 ; icg1 = 0
 629      ABI_ALLOCATE(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ))
 630      ABI_ALLOCATE(smat_inv,(2,dtefield%mband_occ,dtefield%mband_occ))
 631      ABI_ALLOCATE(smat_k_paw,(2,usepaw*dtefield%mband_occ,usepaw*dtefield%mband_occ))
 632      if (calc_epaw3_force .or. calc_epaw3_stress) then ! dsdr needed for forces and stresses in electric field with PAW
 633        ABI_ALLOCATE(dsdr,(2,natom,ncpgr,usepaw*dtefield%mband_occ,usepaw*dtefield%mband_occ))
 634        dsdr = zero
 635      end if
 636 
 637 
 638 !    Loop on the values of ikpt_loc and ikpt1 :
 639 !    ikpt1 is incremented one by one, and number the k points in the FBZ
 640 !    ikpt1i refer to the k point numbering in the IBZ
 641 !    ikpt_loc differs from ikpt1 only in the parallel case, and gives
 642 !    the index of the k point in the FBZ, in the set treated by the present processor
 643 !    NOTE : in order to allow synchronisation, ikpt_loc contain information about
 644 !    ikpt AND ISPPOL !
 645 !    It means that the following loop is equivalent to a double loop :
 646 !    do isppol = 1, nsppol
 647 !    do ikpt1 =  1, dtefield%fmkmem
 648 !
 649      do ikpt_loc = 1, dtefield%fmkmem_max*nsppol
 650 
 651        ikpt1=mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc,1)
 652        isppol=mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc,2)
 653 
 654 !      if this k and spin are for me do it
 655        if (ikpt1 > 0 .and. isppol > 0) then
 656 
 657          ikpt1i = dtefield%indkk_f2ibz(ikpt1,1)
 658          nband_k = dtset%nband(ikpt1i + (isppol-1)*dtset%nkpt)
 659 
 660 !        DEBUG
 661 !        Please keep this debugging feature
 662 !        write(std_out,'(a,5i4)' )' berryphase_new : ikpt_loc,ikpt1,isppol,idir,ifor=',&
 663 !        &                                  ikpt_loc,ikpt1,isppol,idir,ifor
 664 !        ENDDEBUG
 665 
 666          inibz=0
 667          if (dtset%kptns(1,ikpt1i) == dtefield%fkptns(1,ikpt1) .and. &
 668 &         dtset%kptns(2,ikpt1i) == dtefield%fkptns(2,ikpt1) .and. &
 669 &         dtset%kptns(3,ikpt1i) == dtefield%fkptns(3,ikpt1)) inibz=1
 670 
 671          ikg = dtefield%fkgindex(ikpt1)
 672 !        ikpt2 = dtefield%ikpt_dk(ikpt1,ifor,idir)
 673 !        ikpt2i = dtefield%indkk_f2ibz(ikpt2,1)
 674 
 675 !        ikpt3(istep) : index of kpt1 + istep*dk in the FBZ
 676 !        ikpt3i(istep) : index of kpt1 + istep*dk in the IBZ
 677          ikpt3(1) = dtefield%ikpt_dk(ikpt1,ifor,idir)
 678          ikpt3i(1) = dtefield%indkk_f2ibz(ikpt3(1),1)
 679          do istep = 1, berrystep-1
 680            ikpt3(istep+1) = dtefield%ikpt_dk(ikpt3(istep),ifor,idir)
 681            ikpt3i(istep+1) = dtefield%indkk_f2ibz(ikpt3(istep+1),1)
 682          end do
 683 
 684 !        itrs = 0
 685 !        if (dtefield%indkk_f2ibz(ikpt1,6) == 1 ) itrs = itrs + 1
 686 !        if (dtefield%indkk_f2ibz(ikpt2,6) == 1 ) itrs = itrs + 10
 687 
 688          itrs_mult(:)=0
 689          if (dtefield%indkk_f2ibz(ikpt1,6) == 1 ) itrs_mult(:) = itrs_mult(:) + 1
 690          do istep=1,berrystep
 691            if (dtefield%indkk_f2ibz(ikpt3(istep),6) == 1 ) itrs_mult(istep) = itrs_mult(istep) + 10
 692          end do
 693 
 694          npw_k1 = npwarr(ikpt1i)
 695 !        npw_k2 = npwarr(ikpt2i)
 696 
 697          do istep = 1, berrystep
 698            npw_k3(istep)=npwarr(ikpt3i(istep))
 699          end do
 700 
 701 !        ji: the loop is over the FBZ, but sflag and smat only apply to the IBZ
 702          if ( efield_flag .and. inibz == 1) then  !!HONG
 703            ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt
 704            smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt1i_sp,ifor,idir)
 705          else
 706            smat_k(:,:,:) = zero
 707          end if
 708 
 709 !        pwind_k(1:npw_k1) = pwind(ikg+1:ikg+npw_k1,ifor,idir)
 710          pwnsfac_k(1,1:npw_k1) = pwnsfac(1,ikg+1:ikg+npw_k1)
 711          pwnsfac_k(2,1:npw_k1) = pwnsfac(2,ikg+1:ikg+npw_k1)
 712 
 713 !        the array needed to compute the overlap matrix between k and k+istep*dk (with multiple steps)
 714 !        the 0-case (no corresponding pw in k and k+dk) could be handled better (k+2*dk could have a corresponding pw ?)
 715          pwind_k_mult(1:npw_k1,1)=pwind(ikg+1:ikg+npw_k1,ifor,idir)
 716          do istep = 1, berrystep-1
 717            do jj=1, npw_k1
 718              if(pwind_k_mult(jj,istep)/=0)then
 719                pwind_k_mult(jj,istep+1) = pwind(dtefield%fkgindex(ikpt3(istep))+pwind_k_mult(jj,istep),ifor,idir)
 720              else
 721                pwind_k_mult(jj,istep+1) = 0
 722              end if
 723            end do
 724          end do
 725 
 726 !        DEBUG
 727 !        write(std_out,*)' berryphase_new : dtset%berryopt,inibz,ikpt1i,isppol,dtset%nkpt,ifor,idir', &
 728 !        &          dtset%berryopt,inibz,ikpt1i,isppol,dtset%nkpt,ifor,idir
 729 !        write(std_out,'(a,4i4)' )' berryphase_new : sflag_k(:)=',sflag_k(:)
 730 !        ENDDEBUG
 731 
 732          if ( efield_flag .and. inibz == 1) then  !!HONG
 733            ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt
 734            sflag_k(:) = dtefield%sflag(:,ikpt1i_sp,ifor,idir)
 735          else
 736            sflag_k(:) = 0
 737          end if
 738 
 739          if (usepaw == 1) then
 740            icp1=dtefield%cprjindex(ikpt1i,isppol)
 741            call pawcprj_get(atindx1,cprj_k,cprj,natom,1,icp1,ikpt1i,0,isppol,&
 742 &           mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,&
 743 &           my_nspinor,nsppol,0,mpicomm=mpi_enreg%comm_kpt,&
 744 &           proc_distrb=mpi_enreg%proc_distrb)
 745 
 746            if ( ikpt1i /= ikpt1 ) then
 747              call pawcprj_copy(cprj_k,cprj_ikn)
 748              call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,indlmn,&
 749 &             dtefield%indkk_f2ibz(ikpt1,2),dtefield%indkk_f2ibz(ikpt1,6),&
 750 &             dtefield%fkptns(:,dtefield%i2fbz(ikpt1i)),&
 751 &             dtefield%lmax,dtefield%lmnmax,mband,natom,dtefield%mband_occ,my_nspinor,&
 752 &             dtefield%nsym,ntypat,typat,dtefield%zarot)
 753              call pawcprj_copy(cprj_fkn,cprj_k)
 754            end if
 755 
 756          end if ! end if usepaw
 757 
 758 !        DEBUG
 759 !        write(std_out,'(a,4i4)' )' berryphase_new : sflag_k(:)=',sflag_k(:)
 760 !        ENDDEBUG
 761 
 762 !        DEBUG
 763 !        write(std_out,'(a,7i4)')'me, idir,ifor, ikpt_loc, ikpt1, isppol = ',&
 764 !        & me,idir,ifor,ikpt_loc,ikpt1,isppol
 765 !        write(std_out,'(a,10i3)')'pwind_k(1:10) = ',pwind_k(1:10)
 766 !        ENDDEBUG
 767 
 768          do istep=1,berrystep
 769            sflag_k_mult(:,istep) = sflag_k(:)
 770          end do
 771 
 772        end if ! end check that ikpt1 > 0 and isppol > 0
 773 
 774 !      --------------------------------------------------------------------------------
 775 !      Communication
 776 !      --------------------------------------------------------------------------------
 777 
 778        do istep=1,berrystep
 779 
 780 !        if(ikpt_loc <= nsppol*dtefield%fmkmem) then
 781          if (ikpt1 > 0 .and. isppol > 0) then ! I currently have a true kpt to use
 782 
 783            count = npw_k3(istep)*my_nspinor*nband_k
 784            ABI_ALLOCATE(cgq,(2,count))
 785            cgq = zero
 786            source = me
 787            if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt3i(istep),1,nband_k,isppol,me)) then
 788 !            I need the datas from someone else
 789              source = mpi_enreg%proc_distrb(ikpt3i(istep),1,isppol)
 790            end if
 791          else
 792            source = -1 ! I do not have a kpt to use
 793          end if
 794 
 795          do dest = 0, nproc-1
 796 
 797            if ((dest==me) .and. (ikpt1>0) .and. (isppol>0)) then
 798 !            I am destination and I have something to do
 799 !            if (mpi_enreg%paral_compil_kpt == 1) write(std_out,*) &
 800 !            &               'coucou 2, mpi_enreg%proc_distrb(ikpt3i(istep),1:nband_k,isppol) : ', &
 801 !            &               mpi_enreg%proc_distrb(ikpt3i(istep),1:nband_k,isppol)
 802 !            write(std_out,*)'ikpt3i(istep) ', ikpt3i(istep)
 803 !            write(std_out,*)'nband_k ',nband_k
 804 !            write(std_out,*)'isppol ', isppol
 805 !            write(std_out,*)'mpi_enreg%proc_distrb',mpi_enreg%proc_distrb
 806 
 807              if (source == me) then
 808 !              I am destination and source
 809 !              DEBUG
 810 !              write(std_out,*)'copying ... '
 811 !              write(std_out,*)'me: ',me, 'ikpt3i(istep) ', ikpt3i(istep), 'isppol ', isppol
 812 !              ENDDEBUG
 813 
 814 !              pwnsfac
 815                idum = dtefield%fkgindex(ikpt3(istep))
 816                pwnsfac_k(3,1:npw_k3(istep)) = pwnsfac(1,idum+1:idum+npw_k3(istep))
 817                pwnsfac_k(4,1:npw_k3(istep)) = pwnsfac(2,idum+1:idum+npw_k3(istep))
 818 
 819 !              cgq (and cprj)
 820                icg1 = dtefield%cgindex(ikpt3i(istep),isppol)
 821 
 822                if (usepaw == 1) then
 823                  icp2=dtefield%cprjindex(ikpt3i(istep),isppol)
 824                  call pawcprj_get(atindx1,cprj_kb,cprj,natom,1,icp2,ikpt3i(istep),0,isppol,&
 825 &                 mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,my_nspinor,&
 826 &                 nsppol,0,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
 827                end if
 828 
 829                cgq(:,1:count)  = cg(:,icg1+1:icg1+count)
 830 !              if (usepaw == 1) then
 831 !                call pawcprj_copy(cprj_buf,cprj_kb)
 832 !              end if
 833 
 834 !              if ((source /= me)) then
 835              else
 836 !              I am the destination but not the source -> receive
 837 !              DEBUG
 838 !              write(std_out,'(a)')'receiving ...'
 839 !              write(std_out,'(a,i4,a,i4,a,i4,a,i4)')'me: ',me, 'source ', source,'ikpt3i(istep) ', ikpt3i(istep), 'isppol ', isppol
 840 !              ENDDEBUG
 841 
 842 !              receive pwnsfac
 843                ABI_ALLOCATE(buffer,(2,npw_k3(istep)))
 844                tag = ikpt3(istep) + (isppol - 1)*dtefield%fnkpt
 845                call xmpi_recv(buffer,source,tag,spaceComm,ierr)
 846                pwnsfac_k(3,1:npw_k3(istep)) = buffer(1,1:npw_k3(istep))
 847                pwnsfac_k(4,1:npw_k3(istep)) = buffer(2,1:npw_k3(istep))
 848                ABI_DEALLOCATE(buffer)
 849 
 850 !              receive cgq (and cprj)
 851                tag = ikpt3i(istep) + (isppol - 1)*nkpt
 852                call xmpi_recv(cgq,source,tag,spaceComm,ierr)
 853 
 854                if (usepaw == 1) then
 855                  call pawcprj_mpi_recv(natom,n2dim,dimlmn,ncpgr,cprj_kb,source,spaceComm,ierr)
 856                end if
 857 
 858              end if
 859 
 860            else if (dest /= me) then
 861 
 862 !            jkpt is the kpt which is being treated by dest
 863 !            jsppol is his isppol
 864              jkpt = mpi_enreg%kpt_loc2fbz_sp(dest, ikpt_loc,1)
 865              jsppol = mpi_enreg%kpt_loc2fbz_sp(dest, ikpt_loc,2)
 866 
 867              if (jkpt > 0 .and. jsppol > 0) then ! dest is treating a true kpt
 868 
 869                jkpt2 = dtefield%ikpt_dk(jkpt,ifor,idir)
 870                jkpt2i = dtefield%indkk_f2ibz(jkpt2,1)
 871 
 872 !              check if I am his source
 873                if((mpi_enreg%proc_distrb(jkpt2i,1,jsppol) == me))  then
 874 !                I know something about jkpt3i and I must send it
 875 !                DEBUG
 876 !                write(std_out,'(a)')'sending ...'
 877 !                write(std_out,'(a,i4,a,i4,a,i4,a,i4)')'dest: ',dest,' me: ',me,&
 878 !                &                          ' jkpt2i ',jkpt2i,' jsppol: ',jsppol
 879 !                ENDDEBUG
 880 
 881 !                pwnsfac
 882                  tag = jkpt2 + (jsppol - 1)*dtefield%fnkpt
 883                  count1 = npwarr(jkpt2i)
 884                  ABI_ALLOCATE(buffer,(2,count1))
 885                  idum = dtefield%fkgindex(jkpt2)
 886                  buffer(1,1:count1)  = pwnsfac(1,idum+1:idum+count1)
 887                  buffer(2,1:count1)  = pwnsfac(2,idum+1:idum+count1)
 888                  call xmpi_send(buffer,dest,tag,spaceComm,ierr)
 889                  ABI_DEALLOCATE(buffer)
 890 
 891 !                cgq (and cprj)
 892                  icg1 = dtefield%cgindex(jkpt2i,jsppol)
 893 
 894                  if (usepaw == 1) then
 895                    icp2=dtefield%cprjindex(jkpt2i,jsppol)
 896                    call pawcprj_get(atindx1,cprj_buf,cprj,natom,1,icp2,jkpt2i,0,jsppol,&
 897 &                   mband,mkmem,natom,dtefield%mband_occ,dtefield%mband_occ,&
 898 &                   my_nspinor,nsppol,0,mpicomm=mpi_enreg%comm_kpt,&
 899 &                   proc_distrb=mpi_enreg%proc_distrb)
 900                  end if
 901 
 902                  tag = jkpt2i + (jsppol - 1)*nkpt
 903                  count1 = npwarr(jkpt2i)*my_nspinor*nband_k
 904                  ABI_ALLOCATE(buffer,(2,count1))
 905                  buffer(:,1:count1)  = cg(:,icg1+1:icg1+count1)
 906                  call xmpi_send(buffer,dest,tag,spaceComm,ierr)
 907                  ABI_DEALLOCATE(buffer)
 908 
 909                  if (usepaw == 1 ) then
 910                    call pawcprj_mpi_send(natom,n2dim,dimlmn,ncpgr,cprj_buf,dest,spaceComm,ierr)
 911                  end if
 912 
 913                end if ! end check that I am his source
 914 
 915              end if ! end check that jkpt > 0 and jsppol > 0
 916 
 917            end if ! end if statements on dest == me or dest /= me
 918 
 919          end do  ! end loop over dest = 0, nproc - 1
 920 
 921          if (ikpt1 > 0 .and. isppol > 0) then ! if I am treating a kpt, compute the smatrix
 922 
 923            if (usepaw == 1) then
 924              if (ikpt3(istep) /= ikpt3i(istep)) then ! cprj_kb refers to ikpt3i(istep), must compute ikpt3(istep) value
 925                call pawcprj_copy(cprj_kb,cprj_ikn)
 926 
 927                call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,indlmn,&
 928 &               dtefield%indkk_f2ibz(ikpt3(istep),2),dtefield%indkk_f2ibz(ikpt3(istep),6),&
 929 &               dtefield%fkptns(:,dtefield%i2fbz(ikpt3i(istep))),&
 930 &               dtefield%lmax,dtefield%lmnmax,mband,natom,&
 931 &               dtefield%mband_occ,my_nspinor,dtefield%nsym,ntypat,typat,&
 932 &               dtefield%zarot)
 933                call pawcprj_copy(cprj_fkn,cprj_kb)
 934              end if
 935              call smatrix_k_paw(cprj_k,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,typat)
 936 !            write(std_out,'(a,5i4)')' JWZ berryphase_new : ikpt_loc,ikpt1,ikpt1i,ikpt2,ikpt2i ',ikpt_loc,ikpt1,ikpt1i,ikpt3(istep),ikpt3i(istep)
 937 !            call smatrix_k0_paw(atindx1,cprj_k,cprj_k,dtefield,ikpt1i,idir,ifor,&
 938 !            &                                  mband,mpi_enreg,natom,ntypat,pawtab,smat_k_paw,typat)
 939              if (calc_epaw3_force .or. calc_epaw3_stress) then
 940                call dsdr_k_paw(cprj_k,cprj_kb,dsdr,dtefield,idir,ifor,mband,natom,ncpgr,typat)
 941              end if
 942            end if
 943 
 944            icg1 = 0
 945            icg = dtefield%cgindex(ikpt1i,isppol)
 946 !          DEBUG
 947 !          if(istep<=2)then
 948 !          if(ikpt1==1)then
 949 !          write(std_out,'(a,2i4,3e15.4)')'istep ikpt3, kpt, cgq', istep, ikpt3(istep), dtefield%fkptns(:,ikpt3(istep))
 950 !          write(std_out,*) cgq
 951 !          write(std_out,*)
 952 !          end if
 953 !          end if
 954 !          ENDDEBUG
 955            call smatrix(cg,cgq,cg1_k,ddkflag,dtm_k,icg,icg1,itrs_mult(istep),det_inv_smat,maxbd,&
 956 &           mcg,count,mcg1_k,minbd,&
 957 &           mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),&
 958 &           npw_k1,npw_k3(istep),my_nspinor,pwind_k_mult(:,istep),pwnsfac_k,sflag_k_mult(:,istep),&
 959 &           shiftbd,smat_inv,smat_k,smat_k_paw,usepaw)
 960 
 961 ! in finite electric field case with paw must save additional F3 term in forces
 962            if(calc_epaw3_force) then
 963 ! when ncpgr = 3, gradients are wrt to atom displacements
 964 ! but when ncpgr = 9, first 6 gradients are wrt strains, last three are displacements
 965              icpgr_offset = 0
 966              if (ncpgr == 9) icpgr_offset = 6
 967              do iatom = 1, natom
 968                do fdir = 1, 3
 969                  dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = zero
 970                  do iband = 1, dtefield%nband_occ(isppol)
 971                    do jband = 1, dtefield%nband_occ(isppol)
 972 ! collect Im{Trace{S^{-1}.dS/dR}} for this k point
 973                      dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = &
 974 &                     dsdr_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) + &
 975 &                     smat_inv(2,iband,jband)*dsdr(1,iatom,icpgr_offset+fdir,jband,iband) + &
 976 &                     smat_inv(1,iband,jband)*dsdr(2,iatom,icpgr_offset+fdir,jband,iband)
 977                    end do ! end sum over jband
 978                  end do ! end sum over iband
 979                end do ! end sum over fdir
 980              end do ! end sum over iatom
 981            end if ! end check on calc_epaw3_force
 982 
 983 ! in finite electric field case with paw must save additional F3 term in stress
 984 ! note that when strains are present they are always saved before forces
 985 ! therefore no need for icpgr_offset in this case
 986            if(calc_epaw3_stress) then
 987              do iatom = 1, natom
 988                do fdir = 1, 6
 989                  dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = zero
 990                  do iband = 1, dtefield%nband_occ(isppol)
 991                    do jband = 1, dtefield%nband_occ(isppol)
 992 ! collect Im{Trace{S^{-1}.dS/de}} for this k point
 993                      dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) = &
 994 &                     dsds_sum(iatom,fdir,ikpt1+(isppol-1)*dtefield%fnkpt) + &
 995 &                     smat_inv(2,iband,jband)*dsdr(1,iatom,fdir,jband,iband) + &
 996 &                     smat_inv(1,iband,jband)*dsdr(2,iatom,fdir,jband,iband)
 997                    end do ! end sum over jband
 998                  end do ! end sum over iband
 999                end do ! end sum over fdir
1000              end do ! end sum over iatom
1001            end if ! end check on calc_epaw3_stress
1002 
1003            if ((det_inv_smat == 10).or.(det_inv_smat == 11)) then
1004 
1005              if (sqrt(dtm_k(1)*dtm_k(1) + dtm_k(2)*dtm_k(2)) < tol12) then
1006                write(message,'(a,i5,a,a,a)')&
1007 &               '  For k-point #',ikpt1,',',ch10,&
1008 &               '  the determinant of the overlap matrix is found to be 0.'
1009                MSG_BUG(message)
1010              end if
1011 
1012              dtm_mult(1,ikpt1+(isppol-1)*dtefield%fnkpt,istep) = dtm_k(1)
1013              dtm_mult(2,ikpt1+(isppol-1)*dtefield%fnkpt,istep) = dtm_k(2)
1014 
1015            end if
1016 
1017            if ( efield_flag .and. inibz == 1 .and. istep == 1)  then  !!HONG
1018              ikpt1i_sp=ikpt1i+(isppol-1)*dtset%nkpt
1019              dtefield%smat(:,:,:,ikpt1i_sp,ifor,idir) = &
1020 &             smat_k(:,:,:)
1021              dtefield%sflag(:,ikpt1i_sp,ifor,idir) = &
1022 &             sflag_k_mult(:,1)
1023            end if
1024 
1025 ! for IBZ k-points and first step, add
1026            if ((ddkflag==1 .and.((det_inv_smat == 1).or.(det_inv_smat == 11))) .and. inibz == 1 .and. istep == 1) then
1027              cg1(:,icg + 1: icg + npw_k1*my_nspinor*nband_k) = &
1028              cg1(:,icg + 1:icg + npw_k1*my_nspinor*nband_k) + &
1029              dkinv*cg1_k(:,1:npw_k1*my_nspinor*nband_k)
1030            end if
1031 
1032            ABI_DEALLOCATE(cgq)
1033 
1034          end if ! end if ikpt1 > 0 and isppol > 0
1035 
1036        end do ! end loop over istep
1037 
1038 !      if (ikpt_loc <= dtefield%fmkmem) sflag_k(:) = sflag_k_mult(:,1)
1039        if (ikpt1 > 0) sflag_k(:) = sflag_k_mult(:,1)
1040 
1041      end do ! close loop over ikpt_loc (k-points, isppol)
1042 
1043      ABI_DEALLOCATE(smat_inv)
1044      ABI_DEALLOCATE(smat_k)
1045      ABI_DEALLOCATE(smat_k_paw)
1046      if (calc_epaw3_force .or. calc_epaw3_stress) then
1047        ABI_DEALLOCATE(dsdr)
1048      end if
1049 
1050    end do   ! close loop over ifor
1051 
1052 !  MPI communicate stuff between everyone
1053    if (nproc>1) then
1054      count = 2*dtefield%fnkpt*nsppol*berrystep
1055      ABI_ALLOCATE(buffer1,(count))
1056      ABI_ALLOCATE(buffer2,(count))
1057      buffer1(:) = reshape(dtm_mult(:,:,:),(/count/))
1058      call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr)
1059      dtm_mult(:,:,:) = reshape(buffer2(:),(/2,dtefield%fnkpt*nsppol,berrystep/))
1060      ABI_DEALLOCATE(buffer1)
1061      ABI_DEALLOCATE(buffer2)
1062      if (calc_epaw3_force) then
1063        count = natom*3*dtefield%fnkpt*nsppol
1064        ABI_ALLOCATE(buffer1,(count))
1065        ABI_ALLOCATE(buffer2,(count))
1066        buffer1(:) = reshape(dsdr_sum(:,:,:),(/count/))
1067        call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr)
1068        dsdr_sum(:,:,:) = reshape(buffer2(:),(/natom,3,dtefield%fnkpt*nsppol/))
1069        ABI_DEALLOCATE(buffer1)
1070        ABI_DEALLOCATE(buffer2)
1071      end if
1072      if (calc_epaw3_stress) then
1073        count = natom*6*dtefield%fnkpt*nsppol
1074        ABI_ALLOCATE(buffer1,(count))
1075        ABI_ALLOCATE(buffer2,(count))
1076        buffer1(:) = reshape(dsds_sum(:,:,:),(/count/))
1077        call xmpi_sum(buffer1,buffer2,count,spaceComm,ierr)
1078        dsds_sum(:,:,:) = reshape(buffer2(:),(/natom,6,dtefield%fnkpt*nsppol/))
1079        ABI_DEALLOCATE(buffer1)
1080        ABI_DEALLOCATE(buffer2)
1081      end if
1082    end if ! if parallel
1083 
1084 !  DEBUG
1085 !  write(std_out,*)
1086 !  write(std_out,*)'istep = 1, nsppol =',nsppol
1087 !  istep=1
1088 !  isppol=1
1089 !  do jkpt = 1, dtefield%fnkpt
1090 !  write(std_out,'(a,i4,3e15.4,2e15.4)')'jkpt, kpt, dtm_mult(:,kpt,1)', jkpt, dtefield%fkptns(:,jkpt),  dtm_mult(:,jkpt+(isppol-1)*dtefield%fnkpt,istep)
1091 !  end do
1092 !  write(std_out,*)
1093 !  write(std_out,*) "istep = 2"
1094 !  if(berrystep>=2)then
1095 !  istep=2
1096 !  isppol=1
1097 !  do jkpt = 1, dtefield%fnkpt
1098 !  write(std_out,'(a,i4,3e15.4,2e15.4)')'jkpt, kpt, dtm_mult(:,kpt,2)', jkpt, dtefield%fkptns(:,jkpt),  dtm_mult(:,jkpt+(isppol-1)*dtefield%fnkpt,istep)
1099 !  end do
1100 !  end if
1101 !  ENDDEBUG
1102 
1103 !  ===========================================================================
1104 !  in DDK case everything has been calculated above from finite difference
1105 !  Now write the ddk WF to a file
1106 !  ===========================================================================
1107 
1108    if (ddkflag == 1) then
1109 
1110      pertcase = idir + 3*natom
1111      response = 1
1112      call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o)
1113      ABI_ALLOCATE(resid,(mband*nkpt*nsppol))
1114      resid(:) = zero
1115 
1116      call outwf(cg1,dtset,psps,eig_dum,fiwf1o,hdr,kg,dtset%kptns,&
1117 &     mband,mcg,mkmem,mpi_enreg,mpw,natom,dtset%nband,&
1118 &     nkpt,npwarr,nsppol,&
1119 &     occ_dum,resid,response,dtfil%unwff2,wfs,wvl)
1120 
1121      ABI_DEALLOCATE(resid)
1122    end if  ! ddkflag == 1
1123 ! end of ddk part for this idir
1124 
1125 
1126 !  ===========================================================================
1127 !  Compute the Berry phase polarization
1128 !  ===========================================================================
1129 
1130    if (polflag == 1) then
1131 
1132 !    Compute the electronic Berry phase
1133 
1134      polb_mult(:,:)=zero
1135      do istep = 1,berrystep
1136 
1137        if(berrystep==1) then
1138          write(message,'(a,a)')ch10,&
1139 &         ' Compute the electronic contribution to polarization'
1140          call wrtout(std_out,message,'COLL')
1141        else
1142          write(message,'(a,a,i4,a)')ch10,&
1143 &         ' Compute the electronic contribution to polarization for a step of istep=',&
1144 &         istep,'*dk'
1145          call wrtout(std_out,message,'COLL')
1146        end if
1147 
1148        if(istep /= 1) then
1149 !        construct the strings for a step of istep*dk
1150 !        string length
1151          istr=1
1152          nkstr=1
1153          ikpt1=1
1154          do ikpt=1,dtefield%fnkpt
1155            do jstep = 1,istep
1156              ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
1157            end do
1158            if (ikpt1 == 1) exit
1159            nkstr = nkstr + 1
1160          end do
1161 !        Check that the string length is a divisor of nkpt
1162          if(mod(dtefield%fnkpt,nkstr) /= 0) then
1163            write(message,'(a,i5,a,i5,a,i7)')&
1164 &           '  For istep = ', istep,&
1165 &           '  The string length = ',nkstr,&
1166 &           ', is not a divisor of fnkpt =',dtefield%fnkpt
1167            MSG_BUG(message)
1168          end if
1169          nstr = dtefield%fnkpt/nkstr
1170 
1171          write(message,'(a,i1,a,i2,a,i3,a,i6)')&
1172 &         '  berryphase_new: for direction ',idir, ' and istep ', istep, ', nkstr = ',nkstr,&
1173 &         ', nstr = ',nstr
1174          call wrtout(std_out,message,'COLL')
1175          call wrtout(ab_out,message,'COLL')
1176 
1177          ABI_ALLOCATE(idxkstr_mult,(nkstr,nstr))
1178          iunmark = 1
1179          kpt_mark(:)=0
1180          do istr=1,nstr
1181            do while(kpt_mark(iunmark) /= 0)
1182              iunmark = iunmark + 1
1183            end do
1184            idxkstr_mult(1,istr) = iunmark
1185            kpt_mark(iunmark)=1
1186 
1187            ikpt1 = idxkstr_mult(1,istr)
1188            do jkstr=2, nkstr
1189              do jstep = 1, istep
1190                ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
1191              end do
1192              idxkstr_mult(jkstr,istr) = ikpt1
1193              kpt_mark(ikpt1) = 1
1194            end do
1195          end do
1196        else
1197          nstr = dtefield%nstr(idir)
1198          nkstr = dtefield%nkstr(idir)
1199          ABI_ALLOCATE(idxkstr_mult,(nkstr,nstr))
1200          idxkstr_mult(:,:) = dtefield%idxkstr(1:nkstr,1:nstr,idir)
1201        end if
1202 !      DEBUG
1203 !      do istr=1,nstr
1204 !      write(std_out,*)'string ', idxkstr_mult(:,istr)
1205 !      end do
1206 !      ENDBEBUG
1207 
1208        ABI_ALLOCATE(det_string,(2,nstr))
1209        ABI_ALLOCATE(polberry,(nstr))
1210        write(message,'(a,10x,a,10x,a)')ch10,&
1211 &       'istr','polberry(istr)'
1212        call wrtout(std_out,message,'COLL')
1213 
1214        polbtot = zero
1215        do isppol = 1, nsppol
1216 
1217          det_string(1,:) = one ; det_string(2,:) = zero
1218          dtm_k(:) = one
1219          det_average(:) = zero
1220 
1221 
1222          do istr = 1, nstr
1223 
1224            if(calc_epaw3_force) epawf3_str(:,:,:) = zero
1225            if(calc_epaw3_stress) epaws3_str(:,:,:) = zero
1226 
1227            do jkstr = 1, nkstr
1228 
1229              ikpt=idxkstr_mult(jkstr,istr)
1230 
1231              dtm_real=dtm_mult(1,ikpt+(isppol-1)*dtefield%fnkpt,istep)
1232              dtm_imag=dtm_mult(2,ikpt+(isppol-1)*dtefield%fnkpt,istep)
1233 
1234              dtm_k(1) = det_string(1,istr)*dtm_real - &
1235 &             det_string(2,istr)*dtm_imag
1236              dtm_k(2) = det_string(1,istr)*dtm_imag + &
1237 &             det_string(2,istr)*dtm_real
1238              det_string(1:2,istr) = dtm_k(1:2)
1239 !            DEBUG
1240 !            write(std_out,'(a,i4,3e15.4,2e15.4)')'ikpt, kpt, dtm', ikpt, dtefield%fkptns(:,ikpt),  dtm_k
1241 !            ENDDEBUG
1242 
1243              if(calc_epaw3_force) then
1244                do iatom = 1, natom
1245                  do fdir = 1, 3
1246                    epawf3_str(iatom,idir,fdir) = epawf3_str(iatom,idir,fdir) + &
1247 &                   dsdr_sum(iatom,fdir,ikpt+(isppol-1)*dtefield%fnkpt)
1248                  end do ! end loop over fdir
1249                end do ! end loop over natom
1250              end if ! end check on calc_epaw3_force
1251              if(calc_epaw3_stress) then
1252                do iatom = 1, natom
1253                  do fdir = 1, 6
1254                    epaws3_str(iatom,idir,fdir) = epaws3_str(iatom,idir,fdir) + &
1255 &                   dsds_sum(iatom,fdir,ikpt+(isppol-1)*dtefield%fnkpt)
1256                  end do ! end loop over fdir
1257                end do ! end loop over natom
1258              end if ! end check on calc_epaw3_stress
1259 
1260            end do
1261 
1262            if(calc_epaw3_force) then
1263              do iatom = 1, natom
1264                do fdir = 1, 3
1265                  dtefield%epawf3(iatom,idir,fdir) = dtefield%epawf3(iatom,idir,fdir) + &
1266 &                 epawf3_str(iatom,idir,fdir)
1267                end do ! end loop over fdir
1268              end do ! end loop over natom
1269            end if ! end check on calc_epaw3_force
1270            if(calc_epaw3_stress) then
1271              do iatom = 1, natom
1272                do fdir = 1, 6
1273                  dtefield%epaws3(iatom,idir,fdir) = dtefield%epaws3(iatom,idir,fdir) + &
1274 &                 epaws3_str(iatom,idir,fdir)
1275                end do ! end loop over fdir
1276              end do ! end loop over natom
1277            end if ! end check on calc_epaw3_stress
1278 
1279            det_average(:) = det_average(:) + &
1280 &           det_string(:,istr)/dble(nstr)
1281 
1282          end do
1283 
1284 
1285 !        correction to obtain a smouth logarithm of the determinant
1286          ABI_ALLOCATE(str_flag,(nstr))
1287 !        DEBUG
1288 !        since we don't have any case of non-nul Chern number,
1289 !        we must change the det_string value "by brute force" if we want debug this
1290 !        allocate(det_string_test(2,dtefield%nstr(idir)))
1291 !        det_string_test(:,:)=det_string(:,:)
1292 !        kk=0
1293 !        det_string(1,1)=cos(2._dp*Pi*real(kk,dp)/four)
1294 !        det_string(2,1)=sin(2._dp*Pi*real(kk,dp)/four)
1295 !        jj=dtefield%str_neigh(1,1,idir)
1296 !        ll=dtefield%str_neigh(2,1,idir)
1297 !        do while (jj/=1)
1298 !        kk=kk+1
1299 !        det_string(1,jj)=cos(2._dp*Pi*real(kk,dp)/four)
1300 !        det_string(2,jj)=sin(2._dp*Pi*real(kk,dp)/four)
1301 !        det_string(1,ll)=cos(-2._dp*Pi*real(kk,dp)/four)
1302 !        det_string(2,ll)=sin(-2._dp*Pi*real(kk,dp)/four)
1303 !        jj=dtefield%str_neigh(1,jj,idir)
1304 !        ll=dtefield%str_neigh(2,ll,idir)
1305 !        enddo
1306 !        ENDDEBUG
1307          if (istep==1) then
1308            do ineigh_str = 1,2
1309              str_flag(:)=0
1310              delta_str(:) = &
1311 &             dtefield%coord_str(:,dtefield%str_neigh(ineigh_str,1,idir),idir) - dtefield%coord_str(:,1,idir)
1312              dstr(:)= delta_str(:) - nint(delta_str(:)) - real(dtefield%strg_neigh(ineigh_str,1,:,idir),dp)
1313              dist_=0._dp
1314              do kk = 1,2
1315                do jj = 1,2
1316                  dist_ = dist_ + dstr(kk)*dtefield%gmet_str(kk,jj,idir)*dstr(jj)
1317                end do
1318              end do
1319              dist_=sqrt(dist_)
1320              do istr = 1,dtefield%nstr(idir)
1321                if(str_flag(istr)==0)then
1322 !                write(std_out,*)'new string'
1323                  str_flag(istr)=1
1324                  call rhophi(det_string(:,istr),dphase,rho)
1325 !                write(std_out,'(i4,e15.4,e15.4,e15.4)')istr, det_string(:,istr),dphase
1326                  dphase_init=dphase
1327                  jstr = dtefield%str_neigh(ineigh_str,istr,idir)
1328                  do while (istr/=jstr)
1329                    str_flag(jstr)=1
1330                    call rhophi(det_string(:,jstr),dphase_new,rho)
1331                    jj=nint((dphase_new-dphase)/(2._dp*Pi))
1332 !                  DEBUG
1333 !                  write(std_out,'(i4,e15.4,e15.4,e15.4,e15.4,i4)')jstr, det_string(:,jstr),dphase_new,dphase_new-dphase,jj
1334 !                  ENDDEBUG
1335                    dphase_new=dphase_new-two*Pi*real(jj,dp)
1336                    if(jj/=0)then
1337                      write(message,'(6a)') ch10,&
1338 &                     ' berryphase_new : WARNING -',ch10,&
1339 &                     '  the berry phase has some huge variation in the space of strings of k-points',ch10,&
1340 &                     '  ABINIT is trying to correct the berry phase, but it is highly experimental'
1341                      call wrtout(std_out,message,'PERS')
1342                    end if
1343 !                  if(jj/=0)write(std_out,'(i4,e15.4,e15.4,e15.4,e15.4)')jstr, det_string(:,jstr),dphase_new,dphase_new-dphase
1344                    dphase=dphase_new
1345                    jstr=dtefield%str_neigh(ineigh_str,jstr,idir)
1346                  end do
1347 !                write(std_out,*)dphase_init, dphase, (dphase-dphase_init)/(2._dp*Pi),nint((dphase-dphase_init)/(2._dp*Pi))
1348                end if
1349              end do
1350            end do
1351          end if
1352          ABI_DEALLOCATE(str_flag)
1353 !        DEBUG
1354 !        deallocate(dist_str)
1355 !        det_string(:,:)=det_string_test(:,:)
1356 !        deallocate(det_string_test)
1357 !        ENDDEBUG
1358 
1359 !        First berry phase that corresponds to det_average
1360 !        phase0 = atan2(det_average(2),det_average(1))
1361          call rhophi(det_average,phase0,rho)
1362          det_mod = det_average(1)**2+det_average(2)**2
1363 
1364 !        Then berry phase that corresponds to each string relative to the average
1365          do istr = 1, nstr
1366 
1367            rel_string(1) = (det_string(1,istr)*det_average(1) + &
1368            det_string(2,istr)*det_average(2))/det_mod
1369            rel_string(2) = (det_string(2,istr)*det_average(1) - &
1370            det_string(1,istr)*det_average(2))/det_mod
1371 !          dphase = atan2(rel_string(2),rel_string(1))
1372            call rhophi(rel_string,dphase,rho)
1373            polberry(istr) = dtefield%sdeg*(phase0 + dphase)/two_pi
1374            polb_mult(isppol,istep) = polb_mult(isppol,istep) + polberry(istr)/(istep*dtefield%nstr(idir))
1375            polb(isppol) = zero
1376            do jstep=1, istep
1377              polb(isppol)=polb(isppol)+coef(jstep,istep)*polb_mult(isppol,jstep)
1378            end do
1379 
1380            write(message,'(10x,i6,7x,e16.9)')istr,polberry(istr)
1381            call wrtout(std_out,message,'COLL')
1382 
1383          end do
1384 
1385          if(berrystep>1)then
1386            write(message,'(9x,a,7x,e16.9,1x,a,i4,a,i4,a)')&
1387 &           'total',polb_mult(isppol,istep),'(isppol=',isppol,', istep=',istep,')'!,ch10
1388            call wrtout(std_out,message,'COLL')
1389 
1390            write(message,'(3x,a,7x,e16.9,1x,a,i4,a,i4,a,a)')&
1391 &           '+correction',polb(isppol),'(isppol=',isppol,', istep=1..',istep,')',ch10
1392            call wrtout(std_out,message,'COLL')
1393 
1394          else
1395 
1396            write(message,'(9x,a,7x,e16.9,1x,a,i4,a)')&
1397 &           'total',polb_mult(isppol,istep),'(isppol=',isppol,')'!,ch10
1398            call wrtout(std_out,message,'COLL')
1399          end if
1400 
1401          polbtot = polbtot + polb(isppol)
1402 
1403        end do    ! isppol
1404 
1405 !      Fold into interval [-1,1]
1406        polbtot = polbtot - 2_dp*nint(polbtot/2_dp)
1407 
1408        ABI_DEALLOCATE(det_string)
1409        ABI_DEALLOCATE(polberry)
1410 
1411 !      ==========================================================================
1412 
1413 !      Compute the ionic Berry phase
1414 
1415        call xred2xcart(natom,rprimd,xcart,xred)
1416        politot = zero
1417        write(message,'(a)')' Compute the ionic contributions'
1418        call wrtout(std_out,message,'COLL')
1419 
1420        write(message,'(a,2x,a,2x,a,15x,a)')ch10,&
1421 &       'itom', 'itypat', 'polion'
1422        call wrtout(std_out,message,'COLL')
1423 
1424        do iatom = 1, natom
1425          itypat = typat(iatom)
1426 
1427 !        The ionic phase can be computed much easier
1428          polion = zion(itypat)*xred(idir,iatom)
1429 
1430 !        Fold into interval (-1,1)
1431          polion = polion - 2_dp*nint(polion/2_dp)
1432          politot = politot + polion
1433          write(message,'(2x,i2,5x,i2,10x,e16.9)') iatom,itypat,polion
1434          call wrtout(std_out,message,'COLL')
1435        end do
1436 
1437 !      Fold into interval [-1,1] again
1438        politot = politot - 2_dp*nint(politot/2_dp)
1439        pion(idir) = politot
1440 
1441        write(message,'(9x,a,7x,es19.9)') 'total',politot
1442        call wrtout(std_out,message,'COLL')
1443 
1444 
1445 !      ==========================================================================
1446 
1447 !      Compute the total polarization
1448 
1449        poltot = politot + polbtot
1450 
1451        if (berrystep==1)then
1452          write(message,'(a,a)')ch10,&
1453 &         ' Summary of the results'
1454          call wrtout(std_out,message,'COLL')
1455          if (unit_out /= 0) then
1456            call wrtout(unit_out,message,'COLL')
1457          end if
1458        else
1459          write(message,'(a,a,i4)')ch10,&
1460 &         ' Summary of the results for istep =',istep
1461          call wrtout(std_out,message,'COLL')
1462          if (unit_out /= 0) then
1463            call wrtout(unit_out,message,'COLL')
1464          end if
1465        end if
1466 
1467        write(message,'(a,es19.9)')&
1468 &       ' Electronic Berry phase ' ,polbtot
1469        call wrtout(std_out,message,'COLL')
1470        if (unit_out /= 0) then
1471          call wrtout(unit_out,message,'COLL')
1472        end if
1473 
1474        write(message,'(a,es19.9)') &
1475 &       '            Ionic phase ', politot
1476        call wrtout(std_out,message,'COLL')
1477        if (unit_out /= 0) then
1478          call wrtout(unit_out,message,'COLL')
1479        end if
1480 
1481        write(message,'(a,es19.9)') &
1482 &       '            Total phase ', poltot
1483        call wrtout(std_out,message,'COLL')
1484        if (unit_out /= 0) then
1485          call wrtout(unit_out,message,'COLL')
1486        end if
1487 
1488 !      REC start
1489        if(abs(dtset%polcen(idir))>tol8)then
1490          poltot = poltot-dtset%polcen(idir)
1491          write(message,'(a,f15.10)') &
1492 &         '    Translating Polarization by P0 for centrosymmetric cell: ',&
1493 &         dtset%polcen(idir)
1494          call wrtout(std_out,message,'COLL')
1495          if (unit_out /= 0) then
1496            call wrtout(unit_out,message,'COLL')
1497          end if
1498        end if
1499 !      REC end
1500 
1501        poltot = poltot - 2.0_dp*nint(poltot/2._dp)
1502        write(message,'(a,es19.9)') &
1503 &       '    Remapping in [-1,1] ', poltot
1504        call wrtout(std_out,message,'COLL')
1505        if (unit_out /= 0) then
1506          call wrtout(unit_out,message,'COLL')
1507        end if
1508 
1509 !      ! REC and HONG
1510 !      =====================================================================================
1511 !      Polarization branch control  (start)
1512 !      -------------------------------------------------------------------------------------
1513 !      berrysav == 0,  for non fixed D/d calculation, polarizaion is in [-1,1],done above
1514 !      for fixed D/d calculation, choose polarization to minimize internal
1515 !      energy, or minimize |red_efiled|. (red_dfield=red_efiled+red_ptot)
1516 !      (d=e+p, as (26) of Stengel, Suppl.) [[cite:Stengel2009]]
1517 !      This is default value.
1518 !
1519 !      berrysav == 1,  keep the polarization on the same branch, which saved in file POLSAVE
1520 !      ======================================================================================
1521 
1522 !      for fixed D/d calculation, choose polarization to minimize internal energy, or to minimize reduced electric field |red_efield|
1523        if(dtset%berrysav ==0 .and. (dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. &
1524 &       dtset%berryopt == 16 .or. dtset%berryopt == 17))  then
1525 
1526          jump=-nint(dtset%red_dfield(idir) - poltot)   ! red_efield = red_dfield - poltot
1527 
1528          if(jump /= 0)then
1529            write(message,'(a,i1,a,es19.9,a,i2)') &
1530 &           ' P(',idir,') Shifted polarization branch to minimize red_efield &
1531 &           k from ',poltot, ' by ',jump
1532            call wrtout(std_out,message,'COLL')
1533            if (unit_out /= 0) then
1534              call wrtout(unit_out,message,'COLL')
1535            end if
1536            poltot=poltot-jump
1537          end if
1538          pol0(idir)=poltot
1539 
1540        end if
1541 
1542 
1543 !      keep the polarization on the same branch.
1544        if (dtset%berrysav == 1) then
1545 
1546 !        use saved polarization to keep on same branch
1547          inquire(file='POLSAVE',exist=lexist)
1548          if(lexist)then
1549            if(idir==1)then
1550              if(mpi_enreg%me==0)then
1551                if (open_file('POLSAVE',message,newunit=unt,status='OLD') /= 0) then
1552                  MSG_ERROR(message)
1553                end if
1554                read(unt,*)pol0
1555                write(message,'(a,3f20.12)')'Reading old polarization:',pol0
1556                call wrtout(std_out,message,'COLL')
1557                if (unit_out /= 0) then
1558                  call wrtout(unit_out,message,'COLL')
1559                end if
1560                close(unt)
1561              end if
1562              call xmpi_bcast(pol0,0,spaceComm,ierr)
1563            end if
1564          else
1565            pol0(idir)=poltot
1566          end if
1567          jump=nint(poltot-pol0(idir))
1568          if(jump /= 0)then
1569            write(message,'(a,i1,a,es19.9,a,i2)') &
1570 &           ' P(',idir,') jumped to new branch. Shifting bac&
1571 &           k from ',poltot, ' by ',jump
1572            call wrtout(std_out,message,'COLL')
1573            if (unit_out /= 0) then
1574              call wrtout(unit_out,message,'COLL')
1575            end if
1576            poltot=poltot-jump
1577          end if
1578 
1579          pol0(idir)=poltot
1580 
1581        end if
1582 
1583 !      =====================================================================================
1584 !      Polarization branch control  (end)
1585 !      =====================================================================================
1586 
1587 
1588 !      Transform the phase into a polarization
1589        fac = 1._dp/(gmod*dtefield%nkstr(idir))
1590 !      !REC         fac = fac/ucvol
1591 !      !REC         pol = fac*poltot
1592        red_ptot(idir)=poltot !!REC
1593        pol = fac*red_ptot(idir)/ucvol !!REC
1594        ptot(idir)=red_ptot(idir)/ucvol !!REC
1595        write(message,'(a,a,es19.9,a,a,a,es19.9,a,a)')ch10,&
1596 &       '           Polarization ', pol,' (a.u. of charge)/bohr^2',ch10,&
1597 &       '           Polarization ', pol*(e_Cb)/(Bohr_Ang*1d-10)**2,&
1598 &       ' C/m^2',ch10
1599        call wrtout(std_out,message,'COLL')
1600        if (unit_out /= 0) then
1601          call wrtout(unit_out,message,'COLL')
1602        end if
1603 
1604 
1605        ABI_DEALLOCATE(idxkstr_mult)
1606 
1607      end do !istep
1608 
1609      pel(idir) = polbtot
1610 
1611    end if   ! if calculate polarization polflag==1
1612 
1613  end do    ! Close loop over idir
1614 
1615 !!REC start
1616  if (dtset%berrysav == 1) then
1617    if(mpi_enreg%me==0)then
1618      if (open_file('POLSAVE',message,newunit=unt,status='UNKNOWN') /= 0) then
1619        MSG_ERROR(message)
1620      end if
1621      write(unt,'(3F20.12)') pol0
1622      close(unt)
1623    end if
1624    first=.false.
1625  end if
1626 !!REC end
1627 
1628 !-------------------------------------------------
1629 !   Compute polarization in cartesian coordinates
1630 !-------------------------------------------------
1631  if (all(dtset%rfdir(:) == 1)) then
1632 
1633    if(usepaw.ne.1) then
1634      pelev=zero
1635    else
1636      call pawpolev(my_natom,natom,ntypat,pawrhoij,pawtab,pelev,&
1637 &     comm_atom=mpi_enreg%comm_atom)
1638 !    note that in the PAW case, the pelev contribution is already
1639 !    implicitly included in the electronic polarization, from the
1640 !    discretized derivative operator. In the NCPP case no such
1641 !    terms exist anyway. Actually in the PAW formulation
1642 !    such terms are included to all orders, unlike in USPP where only
1643 !    zeroth and first-order terms are. In USPP the first-order term
1644 !    is pelev. Here we compute pelev separately only for reporting
1645 !    purposes in polcart, it is not added into pel or used in the the
1646 !    PAW finite field code in make_grad_berry.F90
1647 !    13 June 2012 J Zwanziger
1648    end if
1649    call polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,3,&
1650 &   ptot_cart,rprimd,ucvol,unit_out,usepaw)
1651 
1652  end if
1653 
1654 !deallocate(pwind_k, dtm)
1655  ABI_DEALLOCATE(pwnsfac_k)
1656  ABI_DEALLOCATE(sflag_k)
1657  ABI_DEALLOCATE(cg1_k)
1658  if (ddkflag == 1)  then
1659    ABI_DEALLOCATE(cg1)
1660    ABI_DEALLOCATE(eig_dum)
1661    ABI_DEALLOCATE(occ_dum)
1662  end if
1663 
1664  if (usepaw == 1) then
1665    ABI_DEALLOCATE(dimlmn)
1666    call pawcprj_free(cprj_k)
1667    call pawcprj_free(cprj_kb)
1668    call pawcprj_free(cprj_gat)
1669    ABI_DATATYPE_DEALLOCATE(cprj_k)
1670    ABI_DATATYPE_DEALLOCATE(cprj_kb)
1671    ABI_DATATYPE_DEALLOCATE(cprj_gat)
1672    if (dtset%kptopt /= 3) then
1673      call pawcprj_free(cprj_ikn)
1674      call pawcprj_free(cprj_fkn)
1675      ABI_DATATYPE_DEALLOCATE(cprj_ikn)
1676      ABI_DATATYPE_DEALLOCATE(cprj_fkn)
1677    end if
1678    if (calc_epaw3_force) then
1679      ABI_DEALLOCATE(dsdr_sum)
1680      ABI_DEALLOCATE(epawf3_str)
1681    end if
1682    if (calc_epaw3_stress) then
1683      ABI_DEALLOCATE(dsds_sum)
1684      ABI_DEALLOCATE(epaws3_str)
1685    end if
1686 
1687    if (nproc>1) then
1688      call pawcprj_free(cprj_buf)
1689      ABI_DATATYPE_DEALLOCATE(cprj_buf)
1690    end if
1691 
1692  end if
1693 
1694  ABI_DEALLOCATE(ikpt3)
1695  ABI_DEALLOCATE(ikpt3i)
1696  ABI_DEALLOCATE(sflag_k_mult)
1697  ABI_DEALLOCATE(npw_k3)
1698  ABI_DEALLOCATE(pwind_k_mult)
1699  ABI_DEALLOCATE(itrs_mult)
1700  ABI_DEALLOCATE(coef)
1701  ABI_DEALLOCATE(polb_mult)
1702  ABI_DEALLOCATE(dtm_mult)
1703 
1704 !DEBUG
1705 !write(std_out,*)'berryphase_new exit'
1706 !END_DEBUG
1707 
1708 end subroutine berryphase_new

ABINIT/init_e_field_vars [ Functions ]

[ Top ] [ Functions ]

NAME

 init_e_field_vars

FUNCTION

 Initialization of variables and data structures used in polarization
 calculations

COPYRIGHT

 Copyright (C) 2004-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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  dtset <type(dataset_type)> = all input variables in this dataset
  gmet(3,3) = reciprocal space metric tensor in bohr**-2
  gprimd(3,3) = primitive translations in recip space
  kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
  mpi_enreg=information about MPI parallelization
  npwarr(nkpt) = number of planewaves in basis and boundary at this k point
  occ(mband*nkpt*nsppol) = occup number for each band at each k point
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3) = dimensional primitive vectors
  symrec(3,3,nsym) = symmetries in reciprocal space in terms of
    reciprocal space primitive translations
  xred(3,natom) = location of atoms in reduced units

OUTPUT

  dtefield <type(efield_type)> :: initialized polarization variables
  pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat
                         between k-points k and k +- dk where dk is
                         parallel to the direction idir
  pwind_alloc = first dimension of pwind and pwnsfac
  pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations

SIDE EFFECTS

 TO DO

NOTES

PARENTS

CHILDREN

      initberry

SOURCE

2725 subroutine init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,&
2726      &              mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,&
2727      &              pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred)
2728 
2729 
2730 !This section has been created automatically by the script Abilint (TD).
2731 !Do not modify the following lines by hand.
2732 #undef ABI_FUNC
2733 #define ABI_FUNC 'init_e_field_vars'
2734 !End of the abilint section
2735 
2736   implicit none
2737 
2738   !Arguments ------------------------------------
2739   !scalars
2740   integer,intent(out) :: pwind_alloc
2741   type(MPI_type),intent(inout) :: mpi_enreg
2742   type(dataset_type),intent(inout) :: dtset
2743   type(efield_type),intent(inout) :: dtefield !vz_i needs efield2
2744   type(pawang_type),intent(in) :: pawang
2745   type(pseudopotential_type),intent(in) :: psps
2746   !arrays
2747   integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
2748   integer,intent(in) :: symrec(3,3,dtset%nsym)
2749   integer,pointer :: pwind(:,:,:)
2750   real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
2751   real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom)
2752   real(dp),pointer :: pwnsfac(:,:)
2753   type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
2754   type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
2755 
2756   !Local variables-------------------------------
2757   logical :: initfield
2758   !scalars
2759 
2760   ! *************************************************************************
2761 
2762   initfield = .false.
2763 
2764   !initialization
2765   dtefield%has_qijb = 0
2766   dtefield%has_epawf3 = 0
2767   dtefield%has_epaws3 = 0
2768   dtefield%has_expibi = 0
2769   dtefield%has_rij = 0
2770   dtefield%usecprj = 0
2771   dtefield%berryopt = 0
2772 
2773   if ((dtset%berryopt < 0).or.(dtset%berryopt == 4) .or. (dtset%berryopt == 6) .or.(dtset%berryopt == 7) .or. &
2774        & (dtset%berryopt == 14) .or.(dtset%berryopt == 16) .or.(dtset%berryopt == 17)) then
2775      nullify(pwind,pwnsfac)
2776      call initberry(dtefield,dtset,gmet,gprimd,kg,&
2777           &   dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,&
2778           &   dtset%natom,dtset%nkpt,npwarr,dtset%nsppol,&
2779           &   dtset%nsym,dtset%ntypat,occ,pawang,pawrad,pawtab,&
2780           &   psps,pwind,pwind_alloc,pwnsfac,rprimd,symrec,&
2781           &   dtset%typat,psps%usepaw,xred)
2782      initfield = .true.
2783   end if
2784 
2785   if (.not. initfield .and. dtset%orbmag == 0) then
2786      ! initorbmag.F90 also allocates pwind and pwnsfac
2787      pwind_alloc = 1
2788      ABI_ALLOCATE(pwind,(pwind_alloc,2,3))
2789      ABI_ALLOCATE(pwnsfac,(2,pwind_alloc))
2790      pwind(:,:,:)=0
2791      pwnsfac(:,:)=zero
2792   end if
2793 
2794 end subroutine init_e_field_vars

ABINIT/initberry [ Functions ]

[ Top ] [ Functions ]

NAME

 initberry

FUNCTION

 Initialization of Berryphase calculation of the polarization, the
 ddk and the response of an insulator to a homogenous electric field.

COPYRIGHT

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

INPUTS

  dtset <type(dataset_type)> = all input variables in this dataset
  gmet(3,3) = reciprocal space metric tensor in bohr**-2
  gprimd(3,3) = primitive translations in recip space
  kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
  mband = maximum number of bands
  mkmem = maximum number of k-points in core memory
  mpw = maximum number of plane waves
  natom = number of atoms in unit cell
  nkpt = number of k points
  npwarr(nkpt) = number of planewaves in basis and boundary at this k point
  nsppol = 1 for unpolarized, 2 for spin-polarized
  nsym = number of symmetry operations
  ntypat = number of types of atoms in unit cell
  occ(mband*nkpt*nsppol) = occup number for each band at each k point
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3) = dimensional primitive vectors
  symrec(3,3,nsym) = symmetries in reciprocal space in terms of
    reciprocal space primitive translations
  typat = typat(natom) list of atom types
  usepaw = flag for PAW (1 PAW, 0 NCPP)
  xred(3,natom) = location of atoms in reduced units

OUTPUT

  dtefield <type(efield_type)> = variables related to Berry phase
      calculations
  pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat
                         between k-points k and k +- dk where dk is
                         parallel to the direction idir
    jpw = pwind(ipw,ifor,idir)
      * ipw = index of plane wave vector G for a given k-point k
      * ifor = 1: k + dk
               2: k - dk
      * idir = direction of the polarization/ddk calculation [dk(idir)
               is the only non-zero element of dk(:)]
      * jpw = index of plane wave vector G (+dG) at k +- dk
              where dG is a shift of one reciprocal lattice vector
              (required to close the strings of k-points using the
               periodic gauge condition)
    In case a G-vector of the basis sphere of plane waves at k
    does not belong to the basis sphere of plane waves at k+dk, jpw = 0.
   pwind_alloc = first dimension of pwind and pwnsfac
   pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations

SIDE EFFECTS

  mpi_enreg = informations about MPI parallelization
    kptdstrb(nproc,nneighbour,fmkmem_max*nsppol) : Array required
      by berryphase_new.f for MPI // over k-points. Defined
      for k-points in the fBZ
      but for k-points in the iBZ. Used by vtorho.f
           nproc = number of cpus
           nneighbour = number of neighbours for each k-point (= 6)

PARENTS

      init_e_field_vars

CHILDREN

      expibi,kpgsph,listkk,metric,pawcprj_alloc,pawcprj_getdim,qijb_kk
      setsym_ylm,smpbz,symatm,timab,wrtout,xmpi_max,xmpi_sum

SOURCE

2878 subroutine initberry(dtefield,dtset,gmet,gprimd,kg,mband,&
2879      &              mkmem,mpi_enreg,mpw,natom,nkpt,npwarr,nsppol,&
2880      &              nsym,ntypat,occ,pawang,pawrad,pawtab,psps,&
2881      &              pwind,pwind_alloc,pwnsfac,&
2882      &              rprimd,symrec,typat,usepaw,xred)
2883 
2884 
2885 !This section has been created automatically by the script Abilint (TD).
2886 !Do not modify the following lines by hand.
2887 #undef ABI_FUNC
2888 #define ABI_FUNC 'initberry'
2889 !End of the abilint section
2890 
2891   implicit none
2892 
2893   !Arguments ------------------------------------
2894   !scalars
2895   integer,intent(in) :: mband,mkmem,mpw,natom,nkpt,nsppol,nsym,ntypat,usepaw
2896   integer,intent(out) :: pwind_alloc
2897   type(MPI_type),intent(inout) :: mpi_enreg
2898   type(dataset_type),intent(inout) :: dtset
2899   type(efield_type),intent(inout) :: dtefield !vz_i
2900   type(pawang_type),intent(in) :: pawang
2901   type(pseudopotential_type),intent(in) :: psps
2902   !arrays
2903   integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
2904   integer,intent(in) :: symrec(3,3,nsym),typat(natom)
2905   integer,pointer :: pwind(:,:,:)
2906   real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol)
2907   real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
2908   real(dp),pointer :: pwnsfac(:,:)
2909   type(pawrad_type),intent(in) :: pawrad(ntypat)
2910   type(pawtab_type),intent(in) :: pawtab(ntypat)
2911 
2912   !Local variables-------------------------------
2913   !scalars
2914   integer :: exchn2n3d,flag,flag_kpt,fnkpt_computed,iband,icg,icprj
2915   integer :: idir,idum,idum1,ierr,ifor,ikg,ikg1,ikpt,ikpt1,ikpt1f
2916   integer :: ikpt1i,ikpt2,ikpt_loc,ikptf,ikpti,ikstr,index,ineigh,ipw,ipwnsfac
2917   integer :: isppol,istr,istwf_k,isym,isym1,itrs,itypat,iunmark,jpw,klmn,lmax,lmn2_size_max
2918   integer :: me,me_g0,mkmem_,my_nspinor,nband_k,mband_occ_k,ncpgr,nkstr,nproc,npw_k,npw_k1,spaceComm
2919   integer :: option, brav, mkpt, nkptlatt
2920   integer :: jstr,ii,jj,isign
2921   integer :: dk_flag, coord1, coord2
2922   integer :: mult
2923   real(dp) :: c1,ecut_eff,eg,eg_ev,rdum,diffk1,diffk2,diffk3
2924   real(dp) :: dist_, max_dist, last_dist, dist,kpt_shifted1,kpt_shifted2,kpt_shifted3
2925   real(dp) :: gprimdlc(3,3),rmetllc(3,3),gmetllc(3,3),ucvol_local
2926   ! gprimd(3,3) = inverse of rprimd
2927   ! rmetlcl(3,3)=real-space metric (same as rmet in metric.F90)
2928   ! gmetlcl(3,3)= same as gmet in metric.F90
2929   ! ucvol = volume of the unit cell in Bohr**3
2930 
2931   character(len=500) :: message
2932   logical :: calc_epaw3_force,calc_epaw3_stress,fieldflag
2933   !arrays
2934   integer :: dg(3),iadum(3),iadum1(3),neigh(6)
2935   integer,allocatable :: dimlmn(:),kg1_k(:,:),kpt_mark(:),nattyp_dum(:)
2936   real(dp) :: diffk(3),dk(3),dum33(3,3),eg_dir(3)
2937   real(dp) :: kpt1(3)
2938   real(dp) :: delta_str3(2), dstr(2),dk_str(2,2,3)
2939   real(dp) :: tsec(2)
2940   real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:),spkpt(:,:)
2941 
2942   ! *************************************************************************
2943 
2944   DBG_ENTER("COLL")
2945 
2946   call timab(1001,1,tsec)
2947   call timab(1002,1,tsec)
2948 
2949   !save the current value of berryopt
2950   dtefield%berryopt = dtset%berryopt
2951   !save the current value of nspinor
2952   dtefield%nspinor = dtset%nspinor
2953 
2954   !----------------------------------------------------------------------------
2955   !-------------------- Obtain k-point grid in the full BZ --------------------
2956   !----------------------------------------------------------------------------
2957 
2958   if(dtset%kptopt==1 .or. dtset%kptopt==2 .or. dtset%kptopt==4)then
2959      !  Compute the number of k points in the G-space unit cell
2960      nkptlatt=dtset%kptrlatt(1,1)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,3) &
2961           &   +dtset%kptrlatt(1,2)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,1) &
2962           &   +dtset%kptrlatt(1,3)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,2) &
2963           &   -dtset%kptrlatt(1,2)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,3) &
2964           &   -dtset%kptrlatt(1,3)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,1) &
2965           &   -dtset%kptrlatt(1,1)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,2)
2966 
2967      !  Call smpbz to obtain the list of k-point in the full BZ - without symmetry reduction
2968      option = 0
2969      brav = 1
2970      mkpt=nkptlatt*dtset%nshiftk
2971      ABI_ALLOCATE(spkpt,(3,mkpt))
2972      call smpbz(1,ab_out,dtset%kptrlatt,mkpt,fnkpt_computed,dtset%nshiftk,option,dtset%shiftk,spkpt)
2973      dtefield%fnkpt = fnkpt_computed
2974      ABI_ALLOCATE(dtefield%fkptns,(3,dtefield%fnkpt))
2975      dtefield%fkptns(:,:)=spkpt(:,1:dtefield%fnkpt)
2976      ABI_DEALLOCATE(spkpt)
2977   else if(dtset%kptopt==3.or.dtset%kptopt==0)then
2978      dtefield%fnkpt=nkpt
2979      ABI_ALLOCATE(dtefield%fkptns,(3,dtefield%fnkpt))
2980      dtefield%fkptns(1:3,1:dtefield%fnkpt)=dtset%kpt(1:3,1:dtefield%fnkpt)
2981      if(dtset%kptopt==0)then
2982         write(message,'(10a)') ch10,&
2983              &     ' initberry : WARNING -',ch10,&
2984              &     '  you have defined manually the k-point grid with kptopt = 0',ch10,&
2985              &     '  the berry phase calculation works only with a regular k-points grid,',ch10,&
2986              &     '  abinit doesn''t check if your grid is regular...'
2987         call wrtout(std_out,message,'PERS')
2988      end if
2989   end if
2990 
2991   !call listkk to get mapping from FBZ to IBZ
2992   rdum=1.0d-5  ! cutoff distance to decide when two k points match
2993   ABI_ALLOCATE(dtefield%indkk_f2ibz,(dtefield%fnkpt,6))
2994 
2995   my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2996 
2997   !ji: The following may need modification in the future
2998   !**** no spin-polarization doubling ; allow use of time reversal symmetry ****
2999 
3000   !Here is original call
3001   !
3002   !call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,&
3003   !& dtefield%fnkpt,dtset%nsym,1,dtset%symafm,dtset%symrel,1)
3004 
3005   call timab(1002,2,tsec)
3006   call timab(1003,1,tsec)
3007 
3008   call listkk(rdum,gmet,dtefield%indkk_f2ibz,dtset%kptns,dtefield%fkptns,nkpt,&
3009        & dtefield%fnkpt,dtset%nsym,1,dtset%symafm,symrec,1,use_symrec=.True.)
3010 
3011   call timab(1003,2,tsec)
3012   call timab(1004,1,tsec)
3013 
3014   !Construct i2fbz and f2ibz
3015   ABI_ALLOCATE(dtefield%i2fbz,(nkpt))
3016   idum=0
3017   do ikpt=1,dtefield%fnkpt
3018      if (dtefield%indkk_f2ibz(ikpt,2)==1 .and. &
3019           &   dtefield%indkk_f2ibz(ikpt,6) == 0 .and. &
3020           &   maxval(abs(dtefield%indkk_f2ibz(ikpt,3:5))) == 0 ) then
3021         dtefield%i2fbz(dtefield%indkk_f2ibz(ikpt,1))=ikpt
3022         idum=idum+1
3023      end if
3024   end do
3025   if (idum/=nkpt)then
3026      message = ' Found wrong number of k-points in IBZ'
3027      MSG_ERROR(message)
3028   end if
3029 
3030   !set flags for fields, forces, stresses
3031   fieldflag = ( (dtset%berryopt== 4) .or. (dtset%berryopt== 6) .or. (dtset%berryopt== 7)  &
3032        & .or. (dtset%berryopt==14) .or. (dtset%berryopt==16) .or. (dtset%berryopt==17) )
3033   ! following two flags activates computation of projector gradient contributions to force and
3034   ! stress in finite field PAW calculations
3035   calc_epaw3_force = (fieldflag .and. usepaw == 1 .and. dtset%optforces /= 0)
3036   calc_epaw3_stress = (fieldflag .and. usepaw == 1 .and. dtset%optstress /= 0)
3037 
3038 
3039 
3040   !----------------------------------------------------------------------------
3041   !------------- Allocate PAW space if necessary ------------------------------
3042   !----------------------------------------------------------------------------
3043 
3044   if (usepaw == 1) then
3045 
3046      dtefield%usepaw   = usepaw
3047      dtefield%natom    = natom
3048      dtefield%my_natom = mpi_enreg%my_natom
3049 
3050      ABI_ALLOCATE(dtefield%lmn_size,(ntypat))
3051      ABI_ALLOCATE(dtefield%lmn2_size,(ntypat))
3052      do itypat = 1, ntypat
3053         dtefield%lmn_size(itypat) = pawtab(itypat)%lmn_size
3054         dtefield%lmn2_size(itypat) = pawtab(itypat)%lmn2_size
3055      end do
3056 
3057      lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2
3058      dtefield%lmn2max = lmn2_size_max
3059 
3060      ! expibi and qijb_kk are NOT parallelized over atoms
3061      ! this may change in the future (JZwanziger 18 March 2014)
3062      ABI_ALLOCATE(dtefield%qijb_kk,(2,lmn2_size_max,dtefield%natom,3))
3063      ABI_ALLOCATE(dtefield%expibi,(2,dtefield%natom,3))
3064      dtefield%has_expibi = 1
3065      dtefield%has_qijb = 1
3066 
3067      if ( fieldflag .and. dtefield%has_rij==0) then
3068         lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2
3069         ABI_ALLOCATE(dtefield%rij,(lmn2_size_max,ntypat,3))
3070         dtefield%has_rij = 1
3071      end if
3072 
3073      ! additional F3-type force term for finite electric field with PAW. Same term
3074      ! might also apply for other displacement-type field calculations, but not sure yet
3075      ! JZwanziger 4 April 2014
3076      if ( calc_epaw3_force ) then
3077         ABI_ALLOCATE(dtefield%epawf3,(dtefield%natom,3,3))
3078         dtefield%has_epawf3 = 1
3079      end if
3080      if ( calc_epaw3_stress ) then
3081         ABI_ALLOCATE(dtefield%epaws3,(dtefield%natom,3,6))
3082         dtefield%has_epaws3 = 1
3083      end if
3084 
3085      ncpgr = 0
3086      if ( fieldflag .and. dtefield%usecprj == 0) then
3087         ABI_ALLOCATE(dimlmn,(natom))
3088         call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,typat,pawtab,'R')
3089         !    allocate space for cprj at kpts in BZ (IBZ or FBZ)
3090         ABI_DATATYPE_ALLOCATE(dtefield%cprj,(natom, mband*dtset%nspinor*dtset%nkpt*nsppol))
3091         !    write(std_out,*) "initberry alloc of cprj ", shape(dtefield%cprj)
3092         if (calc_epaw3_force .and. .not. calc_epaw3_stress) ncpgr = 3
3093         if (.not. calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 6
3094         if (calc_epaw3_force .and. calc_epaw3_stress) ncpgr = 9
3095         call pawcprj_alloc(dtefield%cprj,ncpgr,dimlmn)
3096         dtefield%usecprj = 1
3097         ABI_DEALLOCATE(dimlmn)
3098      end if
3099 
3100      ABI_ALLOCATE(dtefield%cprjindex,(nkpt,nsppol))
3101      dtefield%cprjindex(:,:) = 0
3102 
3103      if (dtset%kptopt /= 3) then
3104         ABI_ALLOCATE(dtefield%atom_indsym,(4,nsym,natom))
3105         call symatm(dtefield%atom_indsym,natom,nsym,symrec,dtset%tnons,tol8,typat,xred)
3106         lmax = psps%mpsang - 1
3107         ABI_ALLOCATE(dtefield%zarot,(2*lmax+1,2*lmax+1,lmax+1,nsym))
3108         call setsym_ylm(gprimd,lmax,nsym,1,rprimd,symrec,dtefield%zarot)
3109         dtefield%nsym = nsym
3110         dtefield%lmax = lmax
3111         dtefield%lmnmax = psps%lmnmax
3112      end if
3113 
3114   end if
3115 
3116   !------------------------------------------------------------------------------
3117   !------------------- Compute variables related to MPI // ----------------------
3118   !------------------------------------------------------------------------------
3119   spaceComm=mpi_enreg%comm_cell
3120   nproc=xmpi_comm_size(spaceComm)
3121   me=xmpi_comm_rank(spaceComm)
3122 
3123   if (nproc==1) then
3124      dtefield%fmkmem = dtefield%fnkpt
3125      dtefield%fmkmem_max = dtefield%fnkpt
3126      dtefield%mkmem_max = nkpt
3127   else
3128      dtefield%fmkmem = 0
3129      do ikpt = 1, dtefield%fnkpt
3130         ikpti = dtefield%indkk_f2ibz(ikpt,1)
3131         nband_k = dtset%nband(ikpti)
3132         if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,-1,me))) &
3133              &     dtefield%fmkmem = dtefield%fmkmem + 1
3134      end do
3135      !  Maximum value of mkmem and fmkmem
3136      call xmpi_max(dtefield%fmkmem,dtefield%fmkmem_max,spaceComm,ierr)
3137      !  I have to use the dummy variable mkmem_ because
3138      !  mkmem is declared as intent(in) while the first
3139      !  argument of xmpi_max must be intent(inout)
3140      mkmem_ = mkmem
3141      call xmpi_max(mkmem_,dtefield%mkmem_max,spaceComm,ierr)
3142   end if
3143 
3144   ABI_ALLOCATE(mpi_enreg%kpt_loc2fbz_sp,(0:nproc-1,1:dtefield%fmkmem_max*nsppol, 1:2))
3145   ABI_ALLOCATE(mpi_enreg%kpt_loc2ibz_sp,(0:nproc-1,1:dtefield%mkmem_max*nsppol, 1:2))
3146   ABI_ALLOCATE(mpi_enreg%kptdstrb,(nproc,6,dtefield%fmkmem_max*nsppol*2))
3147   ABI_ALLOCATE(mpi_enreg%mkmem,(0:nproc-1))
3148   mpi_enreg%kpt_loc2fbz_sp(:,:,:) = 0
3149   mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0
3150   mpi_enreg%kptdstrb(:,:,:)       = 0
3151   mpi_enreg%mkmem(:)              = 0
3152 
3153   if (fieldflag) then
3154      ABI_ALLOCATE(dtefield%cgqindex,(3,6,nkpt*nsppol))
3155      ABI_ALLOCATE(dtefield%nneigh,(nkpt))
3156      dtefield%cgqindex(:,:,:) = 0 ; dtefield%nneigh(:) = 0
3157   end if
3158 
3159   pwind_alloc = mpw*dtefield%fmkmem_max
3160   ABI_ALLOCATE(pwind,(pwind_alloc,2,3))
3161   ABI_ALLOCATE(pwnsfac,(2,pwind_alloc))
3162 
3163   !------------------------------------------------------------------------------
3164   !---------------------- Compute efield_type variables -------------------------
3165   !------------------------------------------------------------------------------
3166 
3167   !Initialization of efield_type variables
3168   mult=dtset%useria+1
3169   dtefield%efield_dot(:) = zero
3170   dtefield%dkvecs(:,:) = zero
3171   dtefield%maxnstr = 0    ; dtefield%maxnkstr  = 0
3172   dtefield%nstr(:) = 0    ; dtefield%nkstr(:) = 0
3173   ABI_ALLOCATE(dtefield%ikpt_dk,(dtefield%fnkpt,2,3))
3174   ABI_ALLOCATE(dtefield%cgindex,(nkpt,nsppol))
3175   ABI_ALLOCATE(dtefield%kgindex,(nkpt))
3176   ABI_ALLOCATE(dtefield%fkgindex,(dtefield%fnkpt))
3177   dtefield%ikpt_dk(:,:,:) = 0
3178   dtefield%cgindex(:,:) = 0
3179   dtefield%mband_occ = 0
3180   ABI_ALLOCATE(dtefield%nband_occ,(nsppol))
3181   dtefield%kgindex(:) = 0
3182   dtefield%fkgindex(:) = 0
3183 
3184   if (fieldflag) then
3185      dtset%rfdir(1:3) = 1
3186   end if
3187 
3188 
3189   !Compute spin degeneracy
3190   if (nsppol == 1 .and. dtset%nspinor == 1) then
3191      dtefield%sdeg = two
3192   else if (nsppol == 2 .or. my_nspinor == 2) then
3193      dtefield%sdeg = one
3194   end if
3195 
3196   !Compute the number of occupied bands and check that
3197   !it is the same for each k-point
3198 
3199   index = 0
3200   do isppol = 1, nsppol
3201      dtefield%nband_occ(isppol) = 0
3202      do ikpt = 1, nkpt
3203 
3204         mband_occ_k = 0
3205         nband_k = dtset%nband(ikpt + (isppol - 1)*nkpt)
3206 
3207         do iband = 1, nband_k
3208            index = index + 1
3209            if (abs(occ(index) - dtefield%sdeg) < tol8) mband_occ_k = mband_occ_k + 1
3210         end do
3211 
3212         if (fieldflag) then
3213            if (nband_k /= mband_occ_k) then
3214               write(message,'(a,a,a)')&
3215                    &         '  In a finite electric field, nband must be equal ',ch10,&
3216                    &         '  to the number of valence bands.'
3217               MSG_ERROR(message)
3218            end if
3219         end if
3220 
3221         if (ikpt > 1) then
3222            if (dtefield%nband_occ(isppol) /= mband_occ_k) then
3223               message = "The number of valence bands is not the same for every k-point of present spin channel"
3224               MSG_ERROR(message)
3225            end if
3226         else
3227            dtefield%mband_occ         = max(dtefield%mband_occ, mband_occ_k)
3228            dtefield%nband_occ(isppol) = mband_occ_k
3229         end if
3230 
3231      end do                ! close loop over ikpt
3232   end do                ! close loop over isppol
3233 
3234   if (fieldflag) then
3235      ABI_ALLOCATE(dtefield%smat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt*nsppol,2,3))
3236 
3237      dtefield%smat(:,:,:,:,:,:) = zero
3238   end if
3239 
3240   ABI_ALLOCATE(dtefield%sflag,(dtefield%mband_occ,nkpt*nsppol,2,3))
3241   dtefield%sflag(:,:,:,:) = 0
3242 
3243   !Compute the location of each wavefunction
3244 
3245   icg = 0
3246   icprj = 0
3247   !ikg = 0
3248   do isppol = 1, nsppol
3249      do ikpt = 1, nkpt
3250 
3251         nband_k = dtset%nband(ikpt + (isppol-1)*nkpt)
3252 
3253         if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
3254 
3255         dtefield%cgindex(ikpt,isppol) = icg
3256         npw_k = npwarr(ikpt)
3257         icg = icg + npw_k*dtefield%nspinor*nband_k
3258 
3259         if (usepaw == 1) then
3260            dtefield%cprjindex(ikpt,isppol) = icprj
3261            icprj = icprj + dtefield%nspinor*nband_k
3262         end if
3263 
3264      end do
3265   end do
3266 
3267   ikg = 0
3268   do ikpt = 1, nkpt
3269      if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)).and.&
3270           &   (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,nsppol,me))) cycle
3271 
3272      npw_k = npwarr(ikpt)
3273      dtefield%kgindex(ikpt) = ikg
3274      ikg = ikg + npw_k
3275   end do
3276 
3277   !Need to use dtset%red_efieldbar in the whole code
3278   !Compute the reciprocal lattice coordinates of the electric field
3279   if (fieldflag) then
3280 
3281      call  metric(gmetllc,gprimdlc,-1,rmetllc,rprimd,ucvol_local)
3282 
3283      if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7) then
3284 
3285         do ii=1,3
3286            dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii))
3287            dtefield%efield_dot(ii) =  dtset%red_efieldbar(ii)
3288         end do
3289 
3290         !    dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
3291         !    dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
3292         !    dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
3293 
3294         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3295              &     ' initberry: Reduced electric field (ebar)',ch10,&
3296              &     '  red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3297         call wrtout(std_out,message,'COLL')
3298 
3299      end if
3300 
3301      if (dtset%berryopt == 6 .or. dtset%berryopt ==7 ) then
3302 
3303         do ii=1,3
3304            dtset%red_dfield(ii)= (dot_product(dtset%dfield(:),gprimdlc(:,ii)))*ucvol_local/(4.d0*pi)
3305         end do
3306 
3307         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3308              &     ' initberry: Reduced electric displacement field',ch10,&
3309              &     '  red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10
3310         call wrtout(std_out,message,'COLL')
3311 
3312      end if
3313 
3314 
3315      if (  dtset%berryopt == 14 ) then
3316         !    transfer to unreduced electric field.
3317         do idir=1,3
3318            dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir))
3319            dtefield%efield_dot(idir) = dtset%red_efieldbar(idir)
3320            !      dtefield%efield2(idir)=dtset%red_efieldbar(idir)
3321         end do
3322 
3323         !    dtefield%efield_dot(1) = dtset%red_efieldbar(1)
3324         !    dtefield%efield_dot(2) = dtset%red_efieldbar(2)
3325         !    dtefield%efield_dot(3) = dtset%red_efieldbar(3)
3326 
3327         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3328              &     ' initberry: Unreduced electric field (a.u.)',ch10,&
3329              &     '  efield(1:3) = ',dtset%efield(1:3),ch10
3330         call wrtout(std_out,message,'COLL')
3331 
3332         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3333              &     ' initberry: Reduced electric field (ebar)',ch10,&
3334              &     '  red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3335         call wrtout(std_out,message,'COLL')
3336 
3337      end if
3338 
3339 
3340      if ( dtset%berryopt == 16 ) then
3341 
3342         !    to calculate D
3343         do ii=1,3
3344            dtset%dfield(ii)  =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,ii))
3345         end do
3346 
3347         do idir=1,3
3348            dtset%efield(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rprimd(:,idir))
3349         end do
3350 
3351         do idir=1,3
3352            dtset%red_efieldbar(idir)= (4*pi/ucvol_local)*dot_product(dtset%red_efield(:),rmetllc(:,idir))
3353            dtefield%efield_dot(idir) = dtset%red_efieldbar(idir)
3354         end do
3355 
3356         !    dtefield%efield_dot(1) = dtset%red_efieldbar(1)
3357         !    dtefield%efield_dot(2) = dtset%red_efieldbar(2)
3358         !    dtefield%efield_dot(3) = dtset%red_efieldbar(3)
3359 
3360 
3361         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3362              &     ' initberry: Unreduced electric displacement field (a.u.)',ch10,&
3363              &     '  dfield(1:3) = ',dtset%dfield(1:3),ch10
3364         call wrtout(std_out,message,'COLL')
3365 
3366         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3367              &     ' initberry: Unreduced electric field (a.u.)',ch10,&
3368              &     '  efield(1:3) = ',dtset%efield(1:3),ch10
3369         call wrtout(std_out,message,'COLL')
3370 
3371         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3372              &     ' initberry: Reduced electric field (ebar)',ch10,&
3373              &     '  red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3374         call wrtout(std_out,message,'COLL')
3375 
3376      end if
3377 
3378      if ( dtset%berryopt ==17) then
3379 
3380         !    to calculate D
3381 
3382         do idir=1,3
3383            dtset%efield(idir)= dot_product(dtset%red_efieldbar(:),gprimdlc(:,idir))  ! from ebar
3384            dtset%dfield(idir)  =(4*pi/ucvol_local)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
3385            !      dtset%red_efield(idir) = (ucvol_local/(4*pi))*dot_product(dtset%red_efieldbar(:),gmetllc(:,idir))
3386            dtefield%efield_dot(idir) = dtset%red_efieldbar(idir)
3387         end do
3388 
3389         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3390              &     ' initberry: Reduced electric field (ebar)',ch10,&
3391              &     '  red_efieldbar(1:3) = ',dtset%red_efieldbar(1:3),ch10
3392         call wrtout(std_out,message,'COLL')
3393 
3394 
3395         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3396              &     ' initberry: Unreduced electric field (a.u.)',ch10,&
3397              &     '  efield(1:3) = ',dtset%efield(1:3),ch10
3398         call wrtout(std_out,message,'COLL')
3399 
3400         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3401              &     ' initberry: Reduced electric displacement field (a.u.)',ch10,&
3402              &     '  red_dfield(1:3) = ',dtset%red_dfield(1:3),ch10
3403         call wrtout(std_out,message,'COLL')
3404 
3405         write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
3406              &     ' initberry: Unreduced electric displacement field (a.u.)',ch10,&
3407              &     '  dfield(1:3) = ',dtset%dfield(1:3),ch10
3408         call wrtout(std_out,message,'COLL')
3409 
3410 
3411      end if
3412 
3413 
3414 
3415   end if
3416 
3417   call timab(1004,2,tsec)
3418 
3419   !------------------------------------------------------------------------------
3420   !---------------------- Compute dk --------------------------------------------
3421   !------------------------------------------------------------------------------
3422 
3423   call timab(1005,1,tsec)
3424 
3425   do idir = 1, 3
3426 
3427      if (dtset%rfdir(idir) == 1) then
3428 
3429         !    Compute dk(:), the vector between a k-point and its nearest
3430         !    neighbour along the direction idir
3431 
3432         dk(:) = zero
3433         dk(idir) = 1._dp   ! 1 mean there is no other k-point un the direction idir
3434         do ikpt = 2, dtefield%fnkpt
3435            diffk(:) = abs(dtefield%fkptns(:,ikpt) - dtefield%fkptns(:,1))
3436            if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.&
3437                 &       (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:)
3438         end do
3439         dtefield%dkvecs(:,idir) = dk(:)
3440         !    DEBUG
3441         !    write(std_out,*)' initberry : idir, dk', idir, dk
3442         !    ENDDEBUG
3443 
3444         !    For each k point, find k_prim such that k_prim= k + dk mod(G)
3445         !    where G is a vector of the reciprocal lattice
3446 
3447         do ikpt = 1, dtefield%fnkpt
3448 
3449            !      First k+dk, then k-dk
3450            do isign=-1,1,2
3451               kpt_shifted1=dtefield%fkptns(1,ikpt)- isign*dk(1)
3452               kpt_shifted2=dtefield%fkptns(2,ikpt)- isign*dk(2)
3453               kpt_shifted3=dtefield%fkptns(3,ikpt)- isign*dk(3)
3454               !        Note that this is still a order fnkpt**2 algorithm.
3455               !        It is possible to implement a order fnkpt algorithm, see listkk.F90.
3456               do ikpt1 = 1, dtefield%fnkpt
3457                  diffk1=dtefield%fkptns(1,ikpt1) - kpt_shifted1
3458                  if(abs(diffk1-nint(diffk1))>tol8)cycle
3459                  diffk2=dtefield%fkptns(2,ikpt1) - kpt_shifted2
3460                  if(abs(diffk2-nint(diffk2))>tol8)cycle
3461                  diffk3=dtefield%fkptns(3,ikpt1) - kpt_shifted3
3462                  if(abs(diffk3-nint(diffk3))>tol8)cycle
3463                  dtefield%ikpt_dk(ikpt,(isign+3)/2,idir) = ikpt1
3464                  exit
3465               end do   ! ikpt1
3466            end do     ! isign
3467 
3468            !      OLD CODING
3469            !      First: k + dk
3470            !      do ikpt1 = 1, dtefield%fnkpt
3471            !      diffk(:) = abs(dtefield%fkptns(:,ikpt1) - &
3472            !      &         dtefield%fkptns(:,ikpt) - dk(:))
3473            !      if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
3474            !      dtefield%ikpt_dk(ikpt,1,idir) = ikpt1
3475            !      exit
3476            !      end if
3477            !      end do
3478 
3479            !      Second: k - dk
3480            !      do ikpt1 = 1, dtefield%fnkpt
3481            !      diffk(:) = abs(dtefield%fkptns(:,ikpt1) - &
3482            !      &         dtefield%fkptns(:,ikpt) + dk(:))
3483            !      if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
3484            !      dtefield%ikpt_dk(ikpt,2,idir) = ikpt1
3485            !      exit
3486            !      end if
3487            !      end do
3488 
3489         end do     ! ikpt
3490 
3491         !    Find the string length, starting from k point 1
3492         !    (all strings must have the same number of points)
3493 
3494         nkstr = 1
3495         ikpt1 = 1
3496         do ikpt = 1, dtefield%fnkpt
3497            ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
3498            if (ikpt1 == 1) exit
3499            nkstr = nkstr + 1
3500         end do
3501 
3502         !    Check that the string length is a divisor of nkpt
3503         if(mod(dtefield%fnkpt,nkstr) /= 0) then
3504            write(message,'(a,i5,a,i7)')&
3505                 &       ' The string length = ',nkstr,&
3506                 &       ', is not a divisor of fnkpt =',dtefield%fnkpt
3507            MSG_BUG(message)
3508         end if
3509 
3510         dtefield%nkstr(idir) = nkstr
3511         dtefield%nstr(idir)  = dtefield%fnkpt/nkstr
3512 
3513      end if      ! dtset%rfdir(idir) == 1
3514 
3515      write(message,'(a,i1,a,i3,a,i6)')&
3516           &   '  initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),&
3517           &   ', nstr = ',dtefield%nstr(idir)
3518      call wrtout(std_out,message,'COLL')
3519      call wrtout(ab_out,message,'COLL')
3520 
3521   end do     ! close loop over idir
3522 
3523   call timab(1005,2,tsec)
3524   call timab(1006,1,tsec)
3525 
3526   dtefield%maxnstr  = maxval(dtefield%nstr(:))
3527   dtefield%maxnkstr = maxval(dtefield%nkstr(:))
3528   ABI_ALLOCATE(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3))
3529   dtefield%idxkstr(:,:,:) = 0
3530 
3531   !for the geometry of the string space :
3532   ABI_ALLOCATE(dtefield%coord_str,(2,dtefield%maxnstr,3))
3533   ABI_ALLOCATE(dtefield%str_neigh,(-2:2,dtefield%maxnstr,3))
3534   ABI_ALLOCATE(dtefield%strg_neigh,(-2:2,dtefield%maxnstr,2,3))
3535   dtefield%coord_str(:,:,:) = 0.d0
3536   dtefield%str_neigh(:,:,:)=0
3537   dtefield%strg_neigh(:,:,:,:)=0
3538   dtefield%gmet_str(:,:,:)=0.d0
3539 
3540   !------------------------------------------------------------------------------
3541   !---------------------- Build the strings -------------------------------------
3542   !------------------------------------------------------------------------------
3543 
3544   ABI_ALLOCATE(kpt_mark,(dtefield%fnkpt))
3545   do idir = 1, 3
3546 
3547      if (dtset%rfdir(idir) == 1) then
3548 
3549         iunmark = 1
3550         kpt_mark(:) = 0
3551         do istr = 1, dtefield%nstr(idir)
3552 
3553            do while(kpt_mark(iunmark) /= 0)
3554               iunmark = iunmark + 1
3555            end do
3556            dtefield%idxkstr(1,istr,idir) = iunmark
3557            kpt_mark(iunmark) = 1
3558            do ikstr = 2, dtefield%nkstr(idir)
3559               ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir)
3560               ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir)
3561               dtefield%idxkstr(ikstr,istr,idir) = ikpt2
3562               kpt_mark(ikpt2) = 1
3563            end do
3564 
3565         end do    ! istr
3566 
3567         !    compute distance between strings
3568         !    compute the metric matrix of the strings space in the direction idir
3569         do ii = 1,3
3570            do jj = 1,3
3571               if (ii<idir.and.jj<idir) dtefield%gmet_str(ii  ,jj  ,idir) = &
3572                    &         gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3573               if (ii<idir.and.jj>idir) dtefield%gmet_str(ii  ,jj-1,idir) = &
3574                    &         gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3575               if (ii>idir.and.jj<idir) dtefield%gmet_str(ii-1,jj  ,idir) = &
3576                    &         gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3577               if (ii>idir.and.jj>idir) dtefield%gmet_str(ii-1,jj-1,idir) = &
3578                    &         gmet(ii,jj) - gmet(ii,idir)*gmet(jj,idir)/gmet(idir,idir)
3579            end do
3580         end do
3581         !    DEBUG
3582         !    write(std_out,*)'gmet'
3583         !    do ii=1,3
3584         !    write(std_out,*)gmet(ii,:)
3585         !    end do
3586         !    write(std_out,*)'gmet_str'
3587         !    do ii=1,2
3588         !    write(std_out,*)dtefield%gmet_str(ii,:,idir)
3589         !    end do
3590         !    ENDDEBUG
3591         do istr = 1, dtefield%nstr(idir)
3592            do ii = 1,3
3593               if (ii<idir) dtefield%coord_str(ii,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir))
3594               if (ii>idir) dtefield%coord_str(ii-1,istr,idir)=dtefield%fkptns(ii,dtefield%idxkstr(1,istr,idir))
3595            end do
3596         end do
3597 
3598         !    the following is very similar to getshell
3599         dist_ = 0._dp
3600         do ii = 1,2
3601            dist_ = dist_ + dtefield%gmet_str(ii,ii,idir)
3602         end do
3603         max_dist = 2._dp * dist_ * 2._dp
3604 
3605         dk_str(:,:,idir) = 0._dp
3606         last_dist = 0._dp
3607         !    ishell = 0
3608         !    dtefield%str_neigh(:,:,:) = 0
3609         dk_flag = 0
3610         do while (dk_flag /= 2)
3611            !      Advance shell counter
3612            !      ishell = ishell + 1
3613 
3614            !      Search the smallest distance between two strings
3615            dist = max_dist
3616            do istr = 1,dtefield%nstr(idir)
3617               delta_str3(:) = dtefield%coord_str(:,1,idir) - dtefield%coord_str(:,istr,idir)
3618               do coord1 = -1,1  !two loop to search also on the border of the BZ
3619                  do coord2 = -1,1
3620                     dist_ = 0._dp
3621                     dstr(:) = delta_str3(:) - nint(delta_str3(:))
3622                     dstr(1) = dstr(1) + real(coord1,dp)
3623                     dstr(2) = dstr(2) + real(coord2,dp)
3624                     do ii = 1,2
3625                        do jj = 1,2
3626                           dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj)
3627                        end do
3628                     end do
3629                     if ((dist_ < dist).and.(dist_ - last_dist > tol8)) then
3630                        dist = dist_
3631                     end if
3632                  end do
3633               end do
3634            end do
3635 
3636            last_dist = dist
3637 
3638            !      search the connecting vectors for that distance
3639            do istr = 1,dtefield%nstr(idir)
3640               delta_str3(:) = dtefield%coord_str(:,istr,idir) - dtefield%coord_str(:,1,idir)
3641               do coord1 = -1,1
3642                  do coord2 = -1,1
3643                     dist_ = 0._dp
3644                     dstr(:) = delta_str3(:) - nint(delta_str3(:))
3645                     dstr(1) = dstr(1) + real(coord1,dp)
3646                     dstr(2) = dstr(2) + real(coord2,dp)
3647                     do ii = 1,2
3648                        do jj = 1,2
3649                           dist_ = dist_ + dstr(ii)*dtefield%gmet_str(ii,jj,idir)*dstr(jj)
3650                        end do
3651                     end do
3652                     if (abs(dist_ - dist) < tol8) then
3653                        if (dk_flag == 0) then
3654                           dk_str(:,1,idir) = dstr(:)
3655                           dk_flag = 1
3656                           !                DEBUG
3657                           !                write(std_out,'(a,i4,2e15.4)')'1st connect', istr, dstr
3658                           !                ENDDEBUG
3659                        elseif (dk_str(1,1,idir)*dstr(2)-dk_str(2,1,idir)*dstr(1) > tol8) then
3660                           dk_str(:,2,idir) = dstr(:)
3661                           dk_flag = 2
3662                           !                DEBUG
3663                           !                write(std_out,'(a,i4,2e15.4)')'2nd connect', istr, dstr
3664                           !                ENDDEBUG
3665                           exit
3666                        end if
3667                     end if
3668                  end do
3669                  if (dk_flag == 2) exit
3670               end do
3671               if (dk_flag == 2) exit
3672            end do
3673 
3674         end do ! do while
3675 
3676         !    search the two neighbours for each string
3677         do istr = 1,dtefield%nstr(idir)
3678            dtefield%str_neigh(0,istr,idir) = istr
3679            dtefield%strg_neigh(0,istr,:,idir) = 0
3680            do jstr = 1,dtefield%nstr(idir)
3681               delta_str3(:) = dtefield%coord_str(:,jstr,idir) - dtefield%coord_str(:,istr,idir)
3682               do coord1 = -1,1
3683                  do coord2 = -1,1
3684                     dist_ = 0._dp
3685                     dstr(:) = delta_str3(:) - nint(delta_str3(:))
3686                     dstr(1) = dstr(1) + real(coord1,dp)
3687                     dstr(2) = dstr(2) + real(coord2,dp)
3688                     do ii = 1,2
3689                        if (sum(abs(dstr(:)-dk_str(:,ii,idir)))<tol8) then
3690                           dtefield%str_neigh(ii,istr,idir) = jstr
3691                           dtefield%strg_neigh(ii,istr,1,idir) = coord1
3692                           dtefield%strg_neigh(ii,istr,2,idir) = coord2
3693                        elseif (sum(abs(dstr(:)+dk_str(:,ii,idir)))<tol8) then
3694                           dtefield%str_neigh(-ii,istr,idir) = jstr
3695                           dtefield%strg_neigh(-ii,istr,1,idir) = coord1
3696                           dtefield%strg_neigh(-ii,istr,2,idir) = coord2
3697                        end if
3698                     end do
3699                  end do
3700               end do
3701            end do
3702         end do
3703 
3704         !    DEBUG
3705         !    write(std_out,'(a,e15.4,e15.4,e15.4,e15.4)')'dk_str',dk_str(1,1,idir),dk_str(2,1,idir),dk_str(1,2,idir),dk_str(2,2,idir)
3706         !    write(std_out,*)'istr, neigh1, strg(1,:), neigh2, strg(2,:),neigh-1, strg(-1,:), neigh-2, strg(-2,:)'
3707         !    do istr=1,dtefield%nstr(idir)
3708         !    write(std_out,'(13i4)')istr, &
3709         !    &       dtefield%str_neigh(1,istr,idir), dtefield%strg_neigh(1,istr,:,idir),&
3710         !    &       dtefield%str_neigh(2,istr,idir), dtefield%strg_neigh(2,istr,:,idir),&
3711         !    &       dtefield%str_neigh(-1,istr,idir), dtefield%strg_neigh(-1,istr,:,idir),&
3712         !    &       dtefield%str_neigh(-2,istr,idir), dtefield%strg_neigh(-2,istr,:,idir)
3713         !    end do
3714         !    ENDDEBUG
3715 
3716 
3717      end if         ! rfdir(idir) == 1
3718 
3719   end do           ! close loop over idir
3720 
3721   ABI_DEALLOCATE(kpt_mark)
3722 
3723   call timab(1006,2,tsec)
3724   call timab(1007,1,tsec)
3725 
3726   !------------------------------------------------------------------------------
3727   !------------ Compute PAW on-site terms if necessary --------------------------
3728   !------------------------------------------------------------------------------
3729 
3730   if (usepaw == 1 .and. dtefield%has_expibi == 1) then
3731      ABI_ALLOCATE(calc_expibi,(2,natom))
3732      do idir = 1, 3
3733         dk = dtefield%dkvecs(1:3,idir)
3734         calc_expibi = zero
3735         call expibi(calc_expibi,dk,natom,xred)
3736         dtefield%expibi(1:2,1:natom,idir) = calc_expibi
3737      end do
3738      !   call expibi(dtefield%expibi,dtefield%dkvecs,natom,xred)
3739      dtefield%has_expibi = 2
3740      ABI_DEALLOCATE(calc_expibi)
3741   end if
3742 
3743   if (usepaw == 1 .and. dtefield%has_qijb == 1) then
3744      ABI_ALLOCATE(calc_qijb,(2,dtefield%lmn2max,natom))
3745 
3746      do idir = 1, 3
3747         dk = dtefield%dkvecs(1:3,idir)
3748         calc_qijb = zero
3749         call qijb_kk(calc_qijb,dk,dtefield%expibi(1:2,1:natom,idir),&
3750              &     gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat)
3751         dtefield%qijb_kk(1:2,1:dtefield%lmn2max,1:natom,idir) = calc_qijb
3752         !    call qijb_kk(dtefield%qijb_kk,dtefield%dkvecs,dtefield%expibi,&
3753         ! &   gprimd,dtefield%lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat)
3754      end do
3755      dtefield%has_qijb = 2
3756      ABI_DEALLOCATE(calc_qijb)
3757   end if
3758 
3759   if (usepaw == 1 .and. dtefield%has_rij == 1) then
3760      c1=sqrt(four_pi/three)
3761      do itypat = 1, ntypat
3762         do klmn = 1, pawtab(itypat)%lmn2_size
3763            dtefield%rij(klmn,itypat,1) = c1*pawtab(itypat)%qijl(4,klmn) ! S_{1,1} ~ x
3764            dtefield%rij(klmn,itypat,2) = c1*pawtab(itypat)%qijl(2,klmn) ! S_{1,-1} ~ y
3765            dtefield%rij(klmn,itypat,3) = c1*pawtab(itypat)%qijl(3,klmn) ! S_{1,0} ~ z
3766         end do ! end loop over klmn
3767      end do ! end loop over itypat
3768      dtefield%has_rij = 2
3769   end if !
3770 
3771   call timab(1007,2,tsec)
3772   call timab(1008,1,tsec)
3773 
3774   !------------------------------------------------------------------------------
3775   !------------ Build the array pwind that is needed to compute the -------------
3776   !------------ overlap matrices at k +- dk                         -------------
3777   !------------------------------------------------------------------------------
3778 
3779   ecut_eff = dtset%ecut*(dtset%dilatmx)**2
3780   exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0
3781   pwind(:,:,:) = 0
3782   pwnsfac(1,:) = 1.0_dp
3783   pwnsfac(2,:) = 0.0_dp
3784   ABI_ALLOCATE(kg1_k,(3,mpw))
3785 
3786   ipwnsfac = 0
3787 
3788   do idir = 1, 3
3789 
3790      if (dtset%rfdir(idir) == 1) then
3791 
3792         dk(:) = dtefield%dkvecs(:,idir)
3793 
3794         do ifor = 1, 2
3795 
3796            if (ifor == 2) dk(:) = -1._dp*dk(:)
3797 
3798            !      Build pwind and kgindex
3799            !      NOTE: The array kgindex is important for parallel execution.
3800            !      In case nsppol = 2, it may happent that a particular processor
3801            !      treats k-points at different spin polarizations.
3802            !      In this case, it is not possible to address the elements of
3803            !      pwind correctly without making use of the kgindex array.
3804 
3805            ikg = 0 ; ikpt_loc = 0 ; isppol = 1
3806            do ikpt = 1, dtefield%fnkpt
3807 
3808               ikpti = dtefield%indkk_f2ibz(ikpt,1)
3809               nband_k = dtset%nband(ikpti)
3810               ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir)
3811               ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
3812 
3813               if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,1,me)).and.&
3814                    &         (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,nsppol,me))) cycle
3815 
3816               ikpt_loc = ikpt_loc + 1
3817 
3818               !        Build basis sphere of plane waves for the nearest neighbour of
3819               !        the k-point (important for MPI //)
3820 
3821               kg1_k(:,:) = 0
3822               kpt1(:) = dtset%kptns(:,ikpt1i)
3823               call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,&
3824                    &         1,mpi_enreg,mpw,npw_k1)
3825               me_g0=mpi_enreg%me_g0
3826 
3827 
3828               !        ji: fkgindex is defined here !
3829               dtefield%fkgindex(ikpt) = ikg
3830 
3831               !
3832               !        Deal with symmetry transformations
3833               !
3834 
3835               !        bra k-point k(b) and IBZ k-point kIBZ(b) related by
3836               !        k(b) = alpha(b) S(b)^t kIBZ(b) + G(b)
3837               !        where alpha(b), S(b) and G(b) are given by indkk_f2ibz
3838               !
3839               !        For the ket k-point:
3840               !        k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k)
3841               !        where GBZ(k) takes k(k) to the BZ
3842               !
3843 
3844               isym  = dtefield%indkk_f2ibz(ikpt,2)
3845               isym1 = dtefield%indkk_f2ibz(ikpt1f,2)
3846 
3847               !        Construct transformed G vector that enters the matching condition:
3848               !        alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) )
3849 
3850               dg(:) = -dtefield%indkk_f2ibz(ikpt,3:5) &
3851                    &         -nint(-dtefield%fkptns(:,ikpt) - dk(:) - tol10 + &
3852                    &         dtefield%fkptns(:,ikpt1f)) &
3853                    &         +dtefield%indkk_f2ibz(ikpt1f,3:5)
3854 
3855               !        old code
3856               !        iadum(:)=0
3857               !        do idum=1,3
3858               !        iadum(:)=iadum(:)+ symrec(:,idum,isym1)*dg(idum)
3859               !        end do
3860 
3861               !        new code
3862               iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:))
3863 
3864               dg(:) = iadum(:)
3865 
3866               if ( dtefield%indkk_f2ibz(ikpt1f,6) == 1 ) dg(:) = -dg(:)
3867 
3868               !        Construct S(k)^{t,-1} S(b)^{t}
3869 
3870               dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym))
3871 
3872               !        Construct alpha(k) alpha(b)
3873 
3874               if (dtefield%indkk_f2ibz(ikpt,6) == dtefield%indkk_f2ibz(ikpt1f,6)) then
3875                  itrs=0
3876               else
3877                  itrs=1
3878               end if
3879 
3880 
3881               npw_k  = npwarr(ikpti)
3882               !        npw_k1 = npwarr(ikpt1i)
3883 
3884               !        loop over bra G vectors
3885               do ipw = 1, npw_k
3886 
3887                  !          NOTE: the bra G vector is taken for the sym-related IBZ k point,
3888                  !          not for the FBZ k point
3889                  iadum(:) = kg(:,dtefield%kgindex(ikpti) + ipw)
3890 
3891                  !          Store non-symmorphic operation phase factor exp[i2\pi \alpha G \cdot t]
3892 
3893                  if ( ipwnsfac == 0 ) then
3894                     !            old code
3895                     rdum=0.0_dp
3896                     do idum=1,3
3897                        rdum=rdum+dble(iadum(idum))*dtset%tnons(idum,isym)
3898                     end do
3899                     rdum=two_pi*rdum
3900                     if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum
3901                     pwnsfac(1,ikg+ipw) = cos(rdum)
3902                     pwnsfac(2,ikg+ipw) = sin(rdum)
3903                     !
3904                     !            new code
3905                     !            rdum = DOT_PRODUCT(dble(iadum(:)),dtset%tnons(:,isym))
3906                     !            rdum= two_pi*rdum
3907                     !            if ( dtefield%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum
3908                     !            pwnsfac(1,ikg+ipw) = cos(rdum)
3909                     !            pwnsfac(2,ikg+ipw) = sin(rdum)
3910 
3911                  end if
3912 
3913                  !          to determine r.l.v. matchings, we transformed the bra vector
3914                  !          Rotation
3915                  iadum1(:)=0
3916                  do idum1=1,3
3917                     iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1)
3918                  end do
3919                  iadum(:)=iadum1(:)
3920                  !          Time reversal
3921                  if (itrs==1) iadum(:)=-iadum(:)
3922                  !          Translation
3923                  iadum(:) = iadum(:) + dg(:)
3924 
3925                  do jpw = 1, npw_k1
3926                     iadum1(1:3) = kg1_k(1:3,jpw)
3927                     if ( (iadum(1) == iadum1(1)).and. &
3928                          &             (iadum(2) == iadum1(2)).and. &
3929                          &             (iadum(3) == iadum1(3)) ) then
3930                        pwind(ikg + ipw,ifor,idir) = jpw
3931                        !              write(std_out,'(a,2x,3i4,2x,i4)') 'Found !:',iadum1(:),jpw
3932                        exit
3933                     end if
3934                  end do
3935               end do
3936 
3937               ikg  = ikg + npw_k
3938 
3939            end do    ! close loop over ikpt
3940 
3941            ipwnsfac = 1
3942 
3943         end do    ! close loop over ifor
3944 
3945      end if      ! rfdir(idir) == 1
3946 
3947   end do        ! close loop over idir
3948 
3949 
3950   call timab(1008,2,tsec)
3951   call timab(1009,1,tsec)
3952 
3953   !Build mpi_enreg%kptdstrb
3954   !array required to communicate the WFs between cpus in berryphase_new.f
3955   !(MPI // over k-points)
3956   if (nproc>1) then
3957      do idir = 1, 3
3958         if (dtset%rfdir(idir) == 1) then
3959            do ifor = 1, 2
3960 
3961               ikpt_loc = 0
3962               do isppol = 1, nsppol
3963 
3964                  do ikpt = 1, dtefield%fnkpt
3965 
3966                     ikpti = dtefield%indkk_f2ibz(ikpt,1)
3967                     nband_k = dtset%nband(ikpti)
3968                     ikpt1f = dtefield%ikpt_dk(ikpt,ifor,idir)
3969                     ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
3970 
3971                     if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle
3972 
3973                     ikpt_loc = ikpt_loc + 1
3974                     mpi_enreg%kptdstrb(me + 1,ifor+2*(idir-1),ikpt_loc) = &
3975                          &             ikpt1i + (isppol - 1)*nkpt
3976 
3977                     mpi_enreg%kptdstrb(me+1,ifor+2*(idir-1),&
3978                          &             ikpt_loc+dtefield%fmkmem_max*nsppol) = &
3979                          &             ikpt1f + (isppol - 1)*dtefield%fnkpt
3980 
3981                  end do   ! ikpt
3982               end do     ! isppol
3983            end do       ! ifor
3984         end if         ! dtset%rfdir(idir) == 1
3985      end do           ! idir
3986   end if             ! nproc>1
3987 
3988   !build mpi_enreg%kpt_loc2fbz_sp
3989   ikpt_loc = 0
3990   do isppol = 1, nsppol
3991      do ikpt = 1, dtefield%fnkpt
3992 
3993         ikpti = dtefield%indkk_f2ibz(ikpt,1)
3994         nband_k = dtset%nband(ikpti)
3995 
3996         if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle
3997 
3998         ikpt_loc = ikpt_loc + 1
3999 
4000         mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 1) = ikpt
4001         mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 2) = isppol
4002 
4003      end do
4004   end do
4005 
4006 
4007   !parallel case only :
4008   !build mpi_enreg%kpt_loc2ibz_sp, dtefield%cgqindex and dtefield%nneigh
4009   if ((fieldflag).and.(nproc>1)) then
4010      ikpt_loc = 0
4011      do isppol = 1, nsppol
4012         do ikpt = 1, nkpt
4013 
4014            ikptf = dtefield%i2fbz(ikpt)
4015            nband_k = dtset%nband(ikpti)
4016 
4017            neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0
4018            do idir=1, 3
4019 
4020               !        skip idir values for which efield_dot(idir) = 0
4021               if (abs(dtefield%efield_dot(idir)) < tol12) cycle
4022 
4023               do ifor = 1, 2
4024 
4025                  flag = 0
4026 
4027                  ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir)
4028                  ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
4029 
4030                  dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg
4031                  ikg = ikg + npwarr(ikpt1i)
4032 
4033                  !          check if this neighbour is also a previous neighbour
4034                  do ineigh = 1, (ifor+2*(idir-1))
4035                     if (neigh(ineigh) == ikpt1i) then
4036                        flag = 1
4037                        dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh
4038                        dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
4039                             &               dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt)
4040                        exit
4041                     end if
4042                  end do
4043                  !          create the cgqindex of the neighbour if necessary
4044                  if (flag == 0) then
4045                     neigh(ifor+2*(idir-1)) = ikpt1i
4046                     dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
4047                          &             ifor+2*(idir-1)
4048                     dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg
4049                     if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1
4050                     icg = icg + npwarr(ikpt1i)*dtefield%nspinor*nband_k
4051                  end if
4052               end do !ifor
4053            end do !idir
4054 
4055            if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me))) then
4056               !        ikpt is one of my kpt_loc
4057               ikpt_loc = ikpt_loc + 1
4058               mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 1) = ikpt
4059               mpi_enreg%kpt_loc2ibz_sp(me, ikpt_loc, 2) = isppol
4060            end if
4061 
4062         end do !ikpt
4063      end do !isppol
4064   end if !nproc>1
4065 
4066   !should be temporary
4067   !unassigned mpi_enreg%kpt_loc2fbz_sp are empty ; inform other cpu (there are better ways...)
4068   mpi_enreg%mkmem(me) = mkmem
4069   !do ii=ikpt_loc+1,dtefield%fmkmem_max
4070   !mpi_enreg%kpt_loc2fbz_sp(me, ii, 1) = -1
4071   !end do
4072 
4073 
4074   !(same as mpi_enreg%kptdstrb but for k-points in the iBZ),
4075   !dtefield%cgqindex and dtefield%nneigh
4076 
4077   if ((fieldflag).and.(nproc>1)) then
4078 
4079      ikpt_loc = 1
4080      do isppol = 1, nsppol
4081         do ikpt = 1, nkpt
4082 
4083            nband_k = dtset%nband(ikpt)
4084            ikptf = dtefield%i2fbz(ikpt)
4085 
4086            neigh(:) = 0 ; icg = 0 ; ikg = 0 ; flag_kpt = 0; icprj = 0
4087            do idir = 1, 3
4088 
4089               !        Skip idir values for which efield_dot(idir) = 0
4090               if (abs(dtefield%efield_dot(idir)) < tol12 .and. (fieldflag)) cycle
4091 
4092               do ifor = 1, 2
4093 
4094                  ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir)
4095                  ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
4096 
4097                  !          dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ikg
4098                  ikg = ikg + npwarr(ikpt1i)
4099 
4100                  flag = 0
4101                  do ineigh = 1, (ifor+2*(idir-1))
4102                     if (neigh(ineigh) == ikpt1i) then
4103                        flag = 1
4104                        !              dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = ineigh
4105                        !              dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
4106                        !              &               dtefield%cgqindex(2,ineigh,ikpt+(isppol-1)*nkpt)
4107                        exit
4108                     end if
4109                  end do
4110                  if (flag == 0) then
4111                     !            neigh(ifor+2*(idir-1)) = ikpt1i
4112                     !            dtefield%cgqindex(1,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = &
4113                     !            &             ifor+2*(idir-1)
4114                     !            dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) = icg
4115                     !            if (isppol == 1) dtefield%nneigh(ikpt) = dtefield%nneigh(ikpt) + 1
4116                     !            icg = icg + npwarr(ikpt1i)*dtset%nspinor*nband_k
4117                  end if
4118 
4119                  if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
4120 
4121                  flag_kpt = 1
4122 
4123                  !          MVeithen: the if condition allows to avoid that the same wavefunction
4124                  !          is send several times to a particular cpu
4125 
4126               end do    ! ifor
4127            end do    ! idir
4128 
4129            if (flag_kpt == 1) ikpt_loc = ikpt_loc + 1
4130 
4131         end do    ! ikpt
4132      end do    ! isppol
4133 
4134   end if   ! fieldflag
4135 
4136   call xmpi_sum(mpi_enreg%kptdstrb,spaceComm,ierr)
4137   call xmpi_sum(mpi_enreg%kpt_loc2fbz_sp,spaceComm,ierr)
4138   if (fieldflag) then
4139      call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr)
4140      call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr)
4141   end if
4142 
4143   !------------------------------------------------------------------------------
4144   !------------------------ Estimate critical field -----------------------------
4145   !------------------------------------------------------------------------------
4146 
4147   !Compute the minimal value of the bandgap required to be below
4148   !the critical field as defined by the relation
4149   !| E_i*a_i | < E_g/n_i
4150 
4151   if (fieldflag) then
4152 
4153      do idir = 1, 3
4154         !    eg_dir(idir) = abs(dtefield%efield_dot(idir))*dtefield%nkstr(idir)
4155         eg_dir(idir) = abs(dtset%red_efieldbar(idir))*dtefield%nkstr(idir)
4156      end do
4157 
4158 
4159      eg = maxval(eg_dir)
4160      eg_ev = eg*Ha_eV
4161 
4162      if (dtset%optcell ==0 .and. (dtset%berryopt == 4 .or. dtset%berryopt == 14)) then
4163         write(message,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,&
4164              &     ' initberry: COMMENT - ',ch10,&
4165              &     '  As a rough estimate,',ch10,&
4166              &     '  to be below the critical field, the bandgap of your system',ch10,&
4167              &     '  should be larger than ',eg_ev,' eV.',ch10
4168         call wrtout(ab_out,message,'COLL')
4169         call wrtout(std_out,message,'COLL')
4170 
4171      else
4172 
4173         write(message,'(a,a,a,a,a,a,a)') ch10,&
4174              &     ' initberry: COMMENT - ',ch10,&
4175              &     '  The estimation of critical electric field should be checked after calculation.',ch10,&
4176              &     '  It is printed out just after total energy.' ,ch10
4177 
4178         call wrtout(ab_out,message,'COLL')
4179         call wrtout(std_out,message,'COLL')
4180 
4181      end if
4182 
4183   end if
4184 
4185   ABI_DEALLOCATE(kg1_k)
4186 
4187   call timab(1009,2,tsec)
4188   call timab(1001,2,tsec)
4189 
4190   DBG_EXIT("COLL")
4191 
4192 end subroutine initberry

ABINIT/m_berryphase_new [ Modules ]

[ Top ] [ Modules ]

NAME

  m_berryphase_new

FUNCTION

COPYRIGHT

  Copyright (C) 2003-2018 ABINIT  group (MVeithen)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_berryphase_new
27 
28  use defs_abitypes
29  use defs_basis
30  use defs_datatypes
31  use defs_wvltypes
32  use m_efield
33  use m_errors
34  use m_abicore
35  use m_xmpi
36 
37  use m_berrytk,      only : smatrix, polcart
38  use m_cgprj,            only : ctocprj
39  use m_fftcore, only : kpgsph
40  use m_geometry,     only : xred2xcart, metric
41  use m_io_tools,     only : open_file
42  use m_iowf,         only : outwf
43  use m_kg,       only : getph
44  use m_kpts,    only : listkk, smpbz
45  use m_mpinfo, only : proc_distrb_cycle
46  use m_numeric_tools,only : rhophi
47  use m_pawang, only : pawang_type
48  use m_pawcprj,      only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, &
49                             pawcprj_put, pawcprj_copy, pawcprj_mpi_recv,  &
50                             pawcprj_mpi_send, pawcprj_free, pawcprj_getdim, pawcprj_symkn
51  use m_paw_dfpt,     only : dsdr_k_paw
52  use m_paw_efield,   only : pawpolev
53  use m_pawrad, only : pawrad_type
54  use m_pawrhoij,     only : pawrhoij_type
55  use m_paw_sphharm, only : setsym_ylm
56  use m_pawtab,       only : pawtab_type
57  use m_paw_overlap,  only : expibi,qijb_kk,smatrix_k_paw
58  use m_symtk,   only : symatm
59  use m_time,    only : timab
60 
61  implicit none
62 
63  private

ABINIT/prtefield [ Functions ]

[ Top ] [ Functions ]

NAME

 prtefield

FUNCTION

 Print components of electric field, displacement field and polarization in nice format

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LBoeri, MT)
 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

  dtset <type(dataset_type)>=all input variables in this dataset
   | berryopt
   | efield
   | dfield
   | red_efield
   | red_efieldbar
   | red_dfield
  dtefield <type(efield_type)>
   | efield2
   | red_ptot1
  iunit = unit number to which the data is printed
  rprimd

OUTPUT

  (only writing)

PARENTS

      gstate,update_e_field_vars

CHILDREN

      metric,wrtout

SOURCE

2315 subroutine prtefield(dtset,dtefield,iunit,rprimd)
2316 
2317 
2318 !This section has been created automatically by the script Abilint (TD).
2319 !Do not modify the following lines by hand.
2320 #undef ABI_FUNC
2321 #define ABI_FUNC 'prtefield'
2322 !End of the abilint section
2323 
2324   implicit none
2325 
2326   !Arguments ------------------------------------
2327   integer :: iunit
2328   real(dp),intent(in) :: rprimd(3,3)
2329   type(efield_type),intent(in) :: dtefield
2330   type(dataset_type),intent(inout) :: dtset
2331 
2332   !Local variables-------------------------------
2333   ! Do not modify the length of this string
2334   !scalars
2335   integer :: idir,ii
2336   character(len=1500) :: message
2337   character(len=7)   :: flag_field(3)
2338 
2339   real(dp) ::    ucvol
2340   ! arrays
2341   real(dp) :: ptot_cart(3),gmet(3,3),gprimd(3,3),rmet(3,3),red_pbar(3),red_dbar(3),red_dfieldbar(3)
2342   real(dp) :: red_efieldbar_lc(3),red_efield_lc(3)
2343 
2344 
2345   ! *************************************************************************
2346 
2347   !DEBUG
2348   !write(iout, '(a)') ' prtefield : enter '
2349   !ENDDEBUG
2350 
2351   !write here
2352 
2353   call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2354 
2355 
2356   ptot_cart(:)=zero
2357   do idir = 1,3
2358      ptot_cart(idir)=rprimd(idir,1)*dtefield%red_ptot1(1) + rprimd(idir,2)*dtefield%red_ptot1(2) + &
2359           &   rprimd(idir,3)*dtefield%red_ptot1(3)
2360   end do
2361   ptot_cart(:)=ptot_cart(:)/ucvol
2362 
2363   if (dtset%berryopt == 4) then
2364 
2365      !  to calculate e Eq.(25)
2366      do idir=1,3
2367         dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir))
2368      end do
2369 
2370      !  to calculate ebar Eq.(25)
2371      do idir=1,3
2372         dtset%red_efieldbar(idir)  =dot_product(dtset%efield(:),rprimd(:,idir))
2373      end do
2374 
2375 
2376      !  to calculate pbar
2377      do idir=1,3
2378         red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2379      end do
2380 
2381      !MGNAG: This message is too long and causes
2382      ! Runtime Error: wrtout_cpp.f90, line 893: Buffer overflow on output
2383      ! with NAG in test seq_tsv6_125 where we write to std_out!
2384      ! I cannot change the RECLEN of std_out!
2385 
2386      write(message,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2387           &   '       E:  ', (dtset%efield(ii), ii=1,3), ch10, &
2388           &   '       P:  ', (ptot_cart(ii), ii=1,3)
2389      call wrtout(iunit,message,'COLL')
2390 
2391      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,&
2392           &   '    ebar:  ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &   !!HONG need to change
2393           &  '    pbar:  ', (red_pbar(ii),ii=1,3)
2394      call wrtout(iunit,message,'COLL')
2395 
2396      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,&
2397           &   '       e:  ', (dtset%red_efield(ii),ii=1,3), ch10, &
2398           &   '       p:  ', (dtefield%red_ptot1(ii), ii=1,3)
2399      call wrtout(iunit,message,'COLL')
2400 
2401      write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10,&
2402           &   ' (S.I.), that is V/m for E, and C/m^2 for P', ch10, &
2403           &   '-      E:  ', (dtset%efield(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10, &    !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2404           &  '       P:  ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3),ch10
2405      call wrtout(iunit,message,'COLL')
2406 
2407   end if ! berryopt ==4
2408 
2409 
2410   if (dtset%berryopt == 6) then
2411      do idir=1,3
2412         dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir))
2413      end do
2414 
2415      !  to calculate ebar   !! Need to be changed
2416      do idir=1,3
2417         dtset%red_efieldbar(idir)  = dot_product(dtset%efield(:),rprimd(:,idir))
2418      end do
2419 
2420      !  to calculate red_pbar
2421      do idir=1,3
2422         red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2423      end do
2424 
2425      !  to calculate red_dbar
2426      do idir=1,3
2427         red_dbar(idir)  = dot_product(dtset%dfield(:),rprimd(:,idir))
2428      end do
2429 
2430      !  to calculate d
2431      do idir=1,3
2432         dtset%red_dfield(idir)  =(ucvol/(4*pi))*dot_product(dtset%dfield(:),gprimd(:,idir))
2433      end do
2434 
2435      write(message,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2436           &   '       e:  ', (dtset%red_efield(ii),ii=1,3), ch10, &
2437           &   '       p:  ', (dtefield%red_ptot1(ii), ii=1,3), ch10, &
2438           &   '       d:  ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2439           &   ' e  +  p:  ', (dtset%red_efield(ii)+dtefield%red_ptot1(ii),ii=1,3)
2440      call wrtout(iunit,message,'COLL')
2441 
2442      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10,&
2443           &   '    ebar:  ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &   !!HONG need to change
2444           &  '    pbar:  ', (red_pbar(ii),ii=1,3), ch10, &
2445           &   '    dbar:  ', (red_dbar(ii),ii=1,3), ch10, &
2446           &   ' eba+pba:  ', (dtset%red_efieldbar(ii)+red_pbar(ii),ii=1,3)
2447      call wrtout(iunit,message,'COLL')
2448 
2449      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2450           &   '       E:  ', (dtset%efield(ii), ii=1,3), ch10, &
2451           &   '       P:  ', (ptot_cart(ii), ii=1,3), ch10, &
2452           &   '       D:  ', (dtset%dfield(ii),ii = 1, 3), ch10, &
2453           &   'E+4*pi*P:  ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3)
2454      call wrtout(iunit,message,'COLL')
2455 
2456      write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10,&
2457 &     ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, &
2458 &     '-      E:  ', (dtset%efield(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10,& !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2459 &     '       P:  ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10,&
2460 &     '       D:  ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3),ch10,&
2461 &     'eps0*E+P:  ', (dtset%efield(ii)*eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))+ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii=1,3),ch10
2462         ! eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))=8.854187817620*5.14220652*1d-1
2463      call wrtout(iunit,message,'COLL')
2464 
2465      !MGNAG Runtime Error: wrtout_cpp.f90, line 896: Buffer overflow on output
2466 
2467   end if  ! berryopt ==6
2468 
2469 
2470   if (dtset%berryopt == 14)  then
2471 
2472      do idir=1,3   ! ebar local
2473         red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir))
2474      end do
2475 
2476 
2477      do idir=1,3
2478         dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2479      end do
2480 
2481      !  to calculate pbar
2482      do idir=1,3
2483         red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2484      end do
2485 
2486      write(message,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))')  ' (a. u.)', ch10,&
2487           &   '   ebar0:  ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &
2488           &   '    ebar:  ', (red_efieldbar_lc(ii),ii=1,3), ch10, &
2489           &   '    pbar:  ', (red_pbar(ii),ii=1,3)
2490      call wrtout(iunit,message,'COLL')
2491 
2492      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2493           &   '       e:  ', (dtset%red_efield(ii),ii=1,3), ch10, &
2494           &   '       p:  ', (dtefield%red_ptot1(ii), ii=1,3)
2495      call wrtout(iunit,message,'COLL')
2496 
2497      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2498           &   '       E:  ', (dtefield%efield2(ii), ii=1,3), ch10, &
2499           &   '       P:  ', (ptot_cart(ii), ii=1,3)
2500      call wrtout(iunit,message,'COLL')
2501 
2502      write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10, &
2503           &   ' (S.I.), that is V/m for E, and C/m^2 for P', ch10, &
2504           &   '-      E:  ', (dtefield%efield2(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10, &    !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2505           &  '       P:  ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3),ch10
2506      call wrtout(iunit,message,'COLL')
2507 
2508 
2509   end if  ! berryopt ==14
2510 
2511 
2512   if (dtset%berryopt == 16) then
2513 
2514      !  to calculate e Eq.(25)
2515      do idir=1,3
2516         dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,idir))
2517      end do
2518 
2519      !  to calculate ebar
2520      do idir=1,3
2521         dtset%red_efieldbar(idir)  = dot_product(dtset%efield(:),rprimd(:,idir))
2522      end do
2523 
2524      !  to calculate pbar
2525      do idir=1,3
2526         red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2527      end do
2528 
2529      !  to calculate dbar
2530      do idir=1,3
2531         red_dfieldbar(idir)  = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir))
2532      end do
2533 
2534      !  to calculate D
2535      do idir=1,3
2536         dtset%dfield(idir)  =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
2537      end do
2538 
2539      write(message,'(a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ' (a. u.)', ch10,&
2540           &   '       e:  ', (dtset%red_efield(ii),ii=1,3), ch10, &
2541           &   '       p:  ', (dtefield%red_ptot1(ii), ii=1,3), ch10, &
2542           &   '       d:  ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2543           &   ' e  +  p:  ', (dtset%red_efield(ii)+dtefield%red_ptot1(ii),ii=1,3)
2544      call wrtout(iunit,message,'COLL')
2545 
2546      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2547           &   '    ebar:  ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &
2548           &   '    pbar:  ', (red_pbar(ii),ii=1,3), ch10, &
2549           &   '    dbar:  ', (red_dfieldbar(ii),ii=1,3), ch10, &
2550           &   ' eba+pba:  ', (dtset%red_efieldbar(ii)+red_pbar(ii),ii=1,3)
2551      call wrtout(iunit,message,'COLL')
2552 
2553      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2554           &   '       E:  ', (dtset%efield(ii), ii=1,3), ch10, &
2555           &   '       P:  ', (ptot_cart(ii), ii=1,3), ch10, &
2556           &   '       D:  ', (dtset%dfield(ii),ii = 1, 3), ch10, &
2557           &   'E+4*pi*P:  ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3)
2558      call wrtout(iunit,message,'COLL')
2559 
2560      write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)') ch10, &
2561 &     ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, &
2562 &     '-      E:  ', (dtset%efield(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10, &    !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2563 &     '       P:  ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10, &
2564 &     '       D:  ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3), ch10, &
2565 &     'eps0*E+P:  ', (dtset%efield(ii)*eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))+ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii=1,3),ch10
2566      call wrtout(iunit,message,'COLL')
2567 
2568   end if  ! berryopt ==16
2569 
2570   if (dtset%berryopt == 17) then
2571 
2572      do idir=1,3   ! ebar local
2573         red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir))
2574      end do
2575 
2576      do idir=1,3
2577         dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2578      end do
2579 
2580      !  to calculate pbar
2581      do idir=1,3
2582         red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2583      end do
2584 
2585 
2586      !  do idir=1,3
2587      !  if (dtset%rfdir(idir)==1) then   ! fixed ebar
2588      !  red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir))  ! local efieldbar
2589      !  dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2590      !  red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2591      !  dtset%dfield(idir)=dtefield%efield2(idir)+4*pi*ptot_cart(idir)
2592      !  dtset%red_dfield(idir)=dtset%red_efield+dtefield%red_ptot1(idir)
2593      !  dtset%red_dfieldbar(idir)=red_efieldbar_lc(idir)+red_pbar(idir)
2594      !  E_lc(idir)=dtefield%efield2(idir)
2595      !  e_lc(idir)=red_efieldbar_lc(idir)
2596      !  ebar_lc(idir)=dtset%red_efieldbar(idir)
2597      !  else if (dtset%rfdir(idir)==2) then ! fixed d
2598      !  dtset%red_efield(idir)  =(ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2599      !  dtset%red_efieldbar(idir)  = dot_product(dtefield%efield2(:),rprimd(:,idir))
2600      !  red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2601      !  red_dfieldbar(idir)  = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir))
2602      !  dtset%dfield(idir)  =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
2603      !  E_lc(idir)=dtefield%efield2(idir)
2604      !  e_lc(idir)=dtset%red_efield(idir)
2605      !  ebar_lc(idir)=dtset%red_efieldbar(idir)
2606      !  end if
2607      !  enddo
2608 
2609 
2610 
2611      do idir=1,3
2612         red_efield_lc(idir)= (ucvol/(4*pi))*dot_product(dtefield%efield2(:),gprimd(:,idir))
2613         red_efieldbar_lc(idir)=dot_product(dtefield%efield2(:),rprimd(:,idir))  ! local efieldbar
2614         red_pbar(idir)  = (4*pi/ucvol)*dot_product(dtefield%red_ptot1(:),rmet(:,idir))
2615         red_dfieldbar(idir)  = (4*pi/ucvol)*dot_product(dtset%red_dfield(:),rmet(:,idir))
2616         dtset%dfield(idir)  =(4*pi/ucvol)*dot_product(dtset%red_dfield(:),rprimd(:,idir))
2617 
2618      end do
2619 
2620      do idir=1,3
2621         if(dtset%jfielddir(idir)==1) then
2622            flag_field(idir)="E-field"
2623         else
2624            flag_field(idir)="D-field"
2625         end if
2626      end do
2627 
2628      write(message,'(a,a,a,6x,a,11x,a,11x,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') &
2629           &   ' (a. u.)', ch10,&
2630           &   '           ', (flag_field(ii),ii=1,3),ch10, &
2631           &   '   ebar0:  ', (dtset%red_efieldbar(ii),ii=1,3), ch10, &
2632           &   '    ebar:  ', (red_efieldbar_lc(ii),ii=1,3), ch10, &
2633           &   '       d:  ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2634           &   ' e  +  p:  ', (red_efield_lc(ii)+dtefield%red_ptot1(ii),ii=1,3)
2635      call wrtout(iunit,message,'COLL')
2636 
2637      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2638           &   '       e:  ', (red_efield_lc(ii),ii=1,3), ch10, &
2639           &   '       p:  ', (dtefield%red_ptot1(ii), ii=1,3), ch10, &
2640           &   '       d:  ', (dtset%red_dfield(ii),ii = 1, 3), ch10, &
2641           &   ' e  +  p:  ', (red_efield_lc(ii)+dtefield%red_ptot1(ii),ii=1,3)
2642      call wrtout(iunit,message,'COLL')
2643 
2644      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2645           &   '    ebar:  ', (red_efieldbar_lc(ii),ii=1,3), ch10, &
2646           &   '    pbar:  ', (red_pbar(ii),ii=1,3), ch10, &
2647           &   '    dbar:  ', (red_dfieldbar(ii),ii=1,3), ch10, &
2648           &   ' eba+pba:  ', (red_efieldbar_lc(ii)+red_pbar(ii),ii=1,3)
2649      call wrtout(iunit,message,'COLL')
2650 
2651      write(message,'(a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x))') ch10, &
2652           &   '       E:  ', (dtefield%efield2(ii), ii=1,3), ch10, &
2653           &   '       P:  ', (ptot_cart(ii), ii=1,3), ch10, &
2654           &   '       D:  ', (dtset%dfield(ii),ii = 1, 3), ch10, &
2655           &   'E+4*pi*P:  ', (dtset%efield(ii)+4.0d0*pi*ptot_cart(ii),ii=1,3)
2656      call wrtout(iunit,message,'COLL')
2657 
2658      write(message,'(a,a,a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a,a,3(es16.9,2x),a)')  ch10, &
2659 &     ' (S.I.), that is V/m for E, and C/m^2 for P and D', ch10, &
2660 &     '       E:  ', (dtefield%efield2(ii)*(Ha_J/(e_Cb*Bohr_Ang*1d-10)), ii=1,3), ch10,& !(Ha_J/(e_Cb*Bohr_Ang*1d-10))= 5.14220652*1d+11
2661 &     '       P:  ', (ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2, ii=1,3), ch10, &
2662 &     '       D:  ', ((1.0d0/(4*pi))*dtset%dfield(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii = 1, 3), ch10, &
2663 &     'eps0*E+P:  ', (dtefield%efield2(ii)*eps0*(Ha_J/(e_Cb*Bohr_Ang*1d-10))+ptot_cart(ii)*(e_Cb)/(Bohr_Ang*1d-10)**2,ii=1,3),ch10
2664      call wrtout(iunit,message,'COLL')
2665 
2666   end if  ! berryopt ==17
2667 
2668 end subroutine prtefield

ABINIT/update_e_field_vars [ Functions ]

[ Top ] [ Functions ]

NAME

 update_e_field_vars

FUNCTION

 This routine updates E field variables

COPYRIGHT

 Copyright (C) 2003-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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

 atindx(natom)=index table for atoms, inverse of atindx (see gstate.f)
 atindx1(natom)=index table for atoms (see gstate.f)
 cg(2,mcg)=planewave coefficients of wavefunctions
 dimcprj(usepaw*natom)=lmn_size for each atom
 dtfil <type(datafiles_type)>=variables related to files
 gmet(3,3)=metric in reciprocal space
 gprimd(3,3)=reciprocal space dimensional primitive translations
 idir = determines directions for derivatives computed in ctocprj (0 for all)
 kg(3,mpw*mkmem)=reduced planewave coordinates
 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
 mkmem=number of k points treated by this node.
 mpw=maximum dimensioned size of npw
 my_natom=number of atoms treated by current processor
 natom=number of atoms in cell
 nattyp(ntypat)=number of atoms of each type
 ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft
 nkpt=number of k-points
 npwarr(nkpt)=number of planewaves in basis at this k point
 ntypat=number of types of atoms in unit cell
 pawrhoij(natom*usepaw) <type(pawrhoij_type)> atomic occupancies
 pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
 psps <type(pseudopotential_type)>=variables related to pseudopotentials
 pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
 pwind_alloc = first dimension of pwind
 pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
 rmet(3,3)=metric in real space
 rprimd(3,3)=dimensional real space primitive translations (bohr)
 scfcv_level= 0 if calling before scf loop, 1 if during
 scfcv_quit=signals whether calling during scf quit (see scfcv.F90)
 scfcv_step=istep value of loop counter from scfcv.F90
 ucvol=unit cell volume in bohr**3.
 unit_out= unit for output of the results (usually the .out file of ABINIT)
   The option unit_out = 0 is allowed. In this case, no information is written
   to the output file but only to the log file.
 usepaw= 1: use paw framework. 0:do not use paw.
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for
     each G and k point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real
     spherical harmonics

OUTPUT

 efield_old_cart(3)=updating cartesian values of efield (used in berryopt
                    6,16,17)
 pel_cg(3)=electronic polarization
 pelev(3)=leading order PAW contribution in pel_cg (for reporting purposes
          only)
 pion(3)=ionic part of polarization
 ptot(3)=total polarization
 red_efield2=updating efield used in berryopt 16,17
 red_efield2_old=updating efield used in berryopt 16.17
 red_ptot=updating efield used in berryopt 16.17

SIDE EFFECTS

 Input/Output
 dtset <type(dataset_type)>=all input variables in this dataset
 dtefield <type(efield_type)> = efield variables
 hdr <type(hdr_type)>=the header of wf, den and pot files
 mpi_enreg=information about MPI parallelization
 ptot_cart(3)=total polarization in cartesian coordinates
 xred(3,natom)=reduced atomic coordinates

TODO

NOTES

PARENTS

CHILDREN

SOURCE

1799 subroutine update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,&
1800 &  efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,&
1801 &  mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,ngfft,nkpt,npwarr,ntypat,&
1802 &  pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,pwind,&
1803 &  pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,&
1804 &  scfcv_level,scfcv_quit,scfcv_step,ucvol,unit_out,&
1805 &  usepaw,xred,ylm,ylmgr)
1806 
1807 
1808 !This section has been created automatically by the script Abilint (TD).
1809 !Do not modify the following lines by hand.
1810 #undef ABI_FUNC
1811 #define ABI_FUNC 'update_e_field_vars'
1812 !End of the abilint section
1813 
1814   implicit none
1815 
1816   !Arguments ------------------------------------
1817   integer, intent(in) :: idir,mcg,mkmem,mpw,my_natom,natom,nkpt,ntypat
1818   integer, intent(in) :: pwind_alloc,scfcv_level,scfcv_quit,scfcv_step,unit_out,usepaw
1819   real(dp), intent(in) :: ucvol
1820   type(datafiles_type), intent(in) :: dtfil
1821   type(pseudopotential_type),intent(in) :: psps
1822   type(dataset_type), intent(inout) :: dtset
1823   type(efield_type), intent(inout) :: dtefield
1824   type(hdr_type), intent(inout) :: hdr
1825   type(MPI_type), intent(inout) :: mpi_enreg
1826   !arrays
1827   integer, intent(in) :: atindx(natom),atindx1(natom),dimcprj(usepaw*natom)
1828   integer, intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat)
1829   integer, intent(in) :: ngfft(18),npwarr(nkpt),pwind(pwind_alloc,2,3)
1830   real(dp), intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3)
1831   real(dp), intent(in) :: pwnsfac(2,pwind_alloc)
1832   real(dp), intent(in) :: rmet(3,3),rprimd(3,3)
1833   real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
1834   real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
1835   real(dp), intent(inout) :: ptot_cart(3),xred(3,natom),efield_old_cart(3) !vz_i
1836   real(dp), intent(out) :: pel_cg(3),pelev(3),pion(3) !vz_i
1837   real(dp), intent(inout) :: red_efield2(3),red_efield2_old(3) !vz_i
1838   real(dp), intent(out) :: ptot(3),red_ptot(3) !vz_i
1839   type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw)
1840   type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
1841 
1842   !Local variables -------------------------
1843   !scalars
1844   character(len=500) :: message
1845   integer :: ctocprj_choice,iatom,ii,iorder_cprj,mcprj,my_nspinor,ncpgr
1846   integer :: optberry,usecprj
1847   logical :: calc_epaw3_force, calc_epaw3_stress, efield
1848   !arrays
1849   real(dp) :: efield_test_cart(3),red_efield1(3)
1850   real(dp),allocatable :: ph1d(:,:)
1851   type(pawcprj_type),allocatable :: cprj(:,:)
1852 
1853   ! *************************************************************************
1854 
1855   efield = .false.
1856 
1857   if ( dtset%berryopt == 4 .or. &
1858        & dtset%berryopt == 6 .or. &
1859        & dtset%berryopt == 7 .or. &
1860        & dtset%berryopt ==14 .or. &
1861        & dtset%berryopt ==16 .or. &
1862        & dtset%berryopt ==17 ) efield = .true.
1863   calc_epaw3_force = ( efield .and. dtset%optforces /= 0 .and. usepaw == 1 )
1864   calc_epaw3_stress = ( efield .and. dtset%optstress /= 0  .and. usepaw == 1 )
1865 
1866   usecprj=1; if (psps%usepaw==0)  usecprj = 0
1867   my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1868   mcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
1869 
1870   ncpgr = 0
1871   ctocprj_choice = 1 ! no derivs
1872   if ( efield .and. psps%usepaw == 1) then
1873      ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj))
1874      !  finite electric field may need gradients for forces, stress
1875      if (calc_epaw3_force .and. .not. calc_epaw3_stress) then
1876         ncpgr = 3; ctocprj_choice = 2 ! derivs w.r.t. position
1877      else if (.not. calc_epaw3_force .and. calc_epaw3_stress) then
1878         ncpgr = 6; ctocprj_choice = 3 ! derivs w.r.t strain
1879      else if (calc_epaw3_force .and. calc_epaw3_stress) then
1880         ncpgr = 9; ctocprj_choice = 23 ! derivs w.r.t. position and strain
1881      end if
1882      call pawcprj_alloc(cprj,ncpgr,dimcprj)
1883      iatom=0 ; iorder_cprj=1 ! retain ordering of input list
1884      !  all arguments to ctocprj are defined already except ph1d, do that here
1885      ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
1886      call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
1887      call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,iorder_cprj,&
1888           &   dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,&
1889           &   mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,nattyp,dtset%nband,&
1890           &   dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,&
1891           &   dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,&
1892           &   dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr)
1893      ABI_DEALLOCATE(ph1d)
1894   else
1895      ABI_DATATYPE_ALLOCATE(cprj,(0,0))
1896   end if ! end update of cprj
1897 
1898   if ( efield ) then ! compute polarization and if necessary store cprj in efield
1899      optberry=1
1900      pel_cg(:) = zero;pelev=zero
1901      call berryphase_new(atindx1,cg,cprj,dtefield,dtfil,dtset,psps,gprimd,hdr,psps%indlmn,kg,&
1902           &   psps%lmnmax,dtset%mband,mcg,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,&
1903           &   dtset%natom,npwarr,dtset%nsppol,psps%ntypat,dtset%nkpt,optberry,pawrhoij,pawtab,&
1904           &   pel_cg,pelev,pion,ptot,red_ptot,pwind,&
1905           &   pwind_alloc,pwnsfac,rprimd,dtset%typat,ucvol,&
1906           &   unit_out,usecprj,psps%usepaw,xred,psps%ziontypat)
1907 
1908      dtefield%red_ptot1(:)=red_ptot(:)
1909 
1910   end if ! end compute polarization and store cprj for efield
1911 
1912   if (efield .and. (scfcv_level == 0) ) then ! do this before scfcv loop
1913 
1914      efield_old_cart(:)=dtset%efield(:)   !!HONG
1915 
1916      !  save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar
1917      dtefield%efield2(:)=dtset%efield(:)
1918 
1919      if ( dtset%berryopt ==16 .or. dtset%berryopt ==17) then   !!HONG
1920         do ii=1,3
1921            red_efield2(ii)=zero
1922            red_efield2_old(ii)  =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii))
1923         end do
1924      end if
1925 
1926      if (dtset%berryopt == 14 .and. scfcv_quit /=1) then
1927         !    ! Convert polarization to cartesian coords
1928 
1929         ptot_cart(:)=zero
1930         do ii = 1,3
1931            ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + &
1932                 &       rprimd(ii,3)*red_ptot(3)
1933         end do
1934         ptot_cart(:)=ptot_cart(:)/ucvol
1935 
1936         do ii=1,3
1937            dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii))
1938         end do
1939 
1940         !    !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
1941         write(message,'(a,a)')   ch10, 'scfcv: Constant reduced ebar-field:'
1942 
1943         call wrtout(std_out,message,'COLL')
1944         call prtefield(dtset,dtefield,std_out,rprimd)
1945 
1946         if(dtset%prtvol>=10)then
1947            call wrtout(ab_out,message,'COLL')
1948            call prtefield(dtset,dtefield,ab_out,rprimd)
1949         end if
1950 
1951         !    updating E field
1952         do ii =1,3   ! desired E field
1953            efield_test_cart(ii)=gprimd(ii,1)*dtset%red_efieldbar(1) + &
1954                 &       gprimd(ii,2)*dtset%red_efieldbar(2)+gprimd(ii,3)*dtset%red_efieldbar(3)
1955         end do
1956 
1957         !    if not convergence well, need to add some code here to make sure efield_test_cart(:) not change much
1958         dtset%efield(:) = efield_test_cart(:)
1959 
1960      end if  ! berryopt ==14
1961 
1962   end if ! end efield .and. scfcv_level 0 tasks
1963 
1964 !!!
1965 !!! Various printing and update steps for the different efield options
1966 !!!
1967 
1968   if (efield .and. (scfcv_level == 1) ) then ! do this each scf step
1969 
1970      if (dtset%prtvol >= 10)then
1971         write(message,'(6(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,&
1972              &     ' scfcv: New value of the polarization:',ch10,&
1973              &     ' (reduced coordinates, a. u.)',ch10,&
1974              &     '     Electronic berry phase:       ', (pel_cg(ii), ii = 1, 3)
1975         call wrtout(ab_out,message,'COLL')
1976         call wrtout(std_out,message,'COLL')
1977         if(psps%usepaw==1) then
1978            write(message,'(a,3(e16.9,2x))')&
1979                 &       '     ...includes PAW on-site term: ', (pelev(ii), ii = 1, 3)
1980            call wrtout(ab_out,message,'COLL')
1981            call wrtout(std_out,message,'COLL')
1982         end if
1983         write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')&
1984              &     '     Ionic:                        ', (pion(ii), ii = 1, 3), ch10, &
1985              &     '     Total:                        ', (red_ptot(ii), ii = 1, 3) !!REC
1986         call wrtout(ab_out,message,'COLL')
1987         call wrtout(std_out,message,'COLL')
1988      end if ! end prtvol >= 10 output
1989 
1990      ptot_cart(:)=zero
1991      do ii = 1,3
1992         ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + &
1993              &     rprimd(ii,3)*red_ptot(3)
1994      end do
1995      ptot_cart(:)=ptot_cart(:)/ucvol
1996 
1997      !  !===================================================================================================
1998      !  !                                       OUTPUT  for fixed E
1999      !  !===================================================================================================
2000 
2001      if (dtset%berryopt == 4) then
2002 
2003         !    !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
2004         write(message,'(a,a)')   ch10, 'scfcv: Constant unreduced E-field:'
2005         call wrtout(std_out,message,'COLL')
2006         call prtefield(dtset,dtefield,std_out,rprimd)
2007         if(dtset%prtvol>=10)then
2008            call wrtout(ab_out,message,'COLL')
2009            call prtefield(dtset,dtefield,ab_out,rprimd)
2010         end if
2011      end if ! end berryopt 4 output
2012 
2013      !  =====================================================================================
2014      !  !                                      fixed D calculation
2015      !  !====================================================================================
2016      if (dtset%berryopt == 6) then
2017         if (scfcv_step > 1) then
2018 
2019            !      ! update efield taking damping into account dfield is in cartesian in dtset structure (contains input value)
2020            !      ! same goes for efield - update the dtset%efield value
2021            efield_test_cart(:)=dtset%ddamp*(dtset%dfield(:)-4.0d0*pi*ptot_cart(:))+&
2022                 &       (1.0d0-dtset%ddamp)*efield_old_cart(:)
2023 
2024            !      ! test whether change in efield in any direction exceed maxestep, if so, set the
2025            !      ! change to maxestep instead   ! need optimized !
2026            do ii = 1,3
2027 
2028              if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then
2029 
2030                write(std_out,'(a,a,i5)') "JH - ","  E-field component:",ii
2031                write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), &
2032 &               ",    E(n-1)=",efield_old_cart(ii), ",    E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), &
2033 &               ",    maxestep=",dtset%maxestep
2034 
2035 
2036                if (efield_test_cart(ii) > efield_old_cart(ii)) then
2037                  efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep)
2038                else
2039                  efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep)
2040                end if
2041              end if
2042            end do
2043 
2044            dtset%efield(:) = efield_test_cart(:)
2045 
2046            !      !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
2047            write(message,'(a,a)')   ch10, 'scfcv: Constant unreduced D-field  - updating E-field:'
2048            call wrtout(std_out,message,'COLL')
2049            call prtefield(dtset,dtefield,std_out,rprimd)
2050            if(dtset%prtvol>=10)then
2051               call wrtout(ab_out,message,'COLL')
2052               call prtefield(dtset,dtefield,ab_out,rprimd)
2053            end if
2054 
2055            !      ! need to update dtset%efield_dot(:) with new value
2056            dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
2057            dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
2058            dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
2059 
2060         else
2061 
2062            write(message,'(a,a)')   ch10, 'scfcv: Constant unreduced D-field  - Pre E-field:'
2063            call wrtout(std_out,message,'COLL')
2064            call prtefield(dtset,dtefield,std_out,rprimd)
2065            if(dtset%prtvol>=10)then
2066               call wrtout(ab_out,message,'COLL')
2067               call prtefield(dtset,dtefield,ab_out,rprimd)
2068            end if
2069 
2070         end if  ! scfcv_step >1
2071 
2072         efield_old_cart(:)=dtset%efield(:)
2073      end if  ! berryopt ==6
2074      !  !===================================================================================================
2075      !  !                                      fixed reduced d calculation
2076      !  !===================================================================================================
2077      if (dtset%berryopt == 16) then
2078 
2079         if (scfcv_step > 1) then
2080            !      ! update efield taking damping into account reduced red_dfield
2081            !      red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]]
2082 
2083            red_efield2(:)=dtset%ddamp*(dtset%red_dfield(:)-red_ptot(:))+ (1.0d0-dtset%ddamp)*red_efield2_old(:)
2084 
2085            !      to calculate unreduced E
2086            efield_test_cart(:)=(4*pi/ucvol)*(rprimd(:,1)*red_efield2(1)+rprimd(:,2)*red_efield2(2)+rprimd(:,3)*red_efield2(3))
2087 
2088            !      ! test whether change in efield in any direction exceed maxestep, if so, set the
2089            !      ! change to maxestep instead   ! need optimized !
2090            do ii = 1,3
2091 
2092              if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then
2093 
2094                write(std_out,'(a,a,i5)') "JH - ","  E-field component:",ii
2095                write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), &
2096 &               ",    E(n-1)=",efield_old_cart(ii), ",    E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), &
2097 &               ",    maxestep=",dtset%maxestep
2098 
2099                if (efield_test_cart(ii) > efield_old_cart(ii)) then
2100                  efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep)
2101                else
2102                  efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep)
2103                end if
2104              end if
2105            end do
2106 
2107            dtset%efield(:) = efield_test_cart(:)
2108 
2109            !      !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
2110            write(message,'(a,a)')   ch10, 'scfcv: Constant reduced d-field  - updating E-field:'
2111            call wrtout(std_out,message,'COLL')
2112            call prtefield(dtset,dtefield,std_out,rprimd)
2113            if(dtset%prtvol>=10)then
2114               call wrtout(ab_out,message,'COLL')
2115               call prtefield(dtset,dtefield,ab_out,rprimd)
2116            end if
2117 
2118            !      ! need to update dtset%efield_dot(:) with new value
2119            !      ! This needs to be deleted  when efield_dot is deleted
2120            dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
2121            dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
2122            dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
2123 
2124         else
2125 
2126            write(message,'(a,a)')   ch10, 'scfcv: Constant reduced d-field  - Pre E-field:'
2127            call wrtout(std_out,message,'COLL')
2128            call prtefield(dtset,dtefield,std_out,rprimd)
2129            if(dtset%prtvol>=10)then
2130               call wrtout(ab_out,message,'COLL')
2131               call prtefield(dtset,dtefield,ab_out,rprimd)
2132            end if
2133 
2134         end if  ! scfcv_step > 1
2135 
2136         efield_old_cart(:)=dtset%efield(:)
2137         red_efield2_old(:)=red_efield2(:)
2138      end if  ! berryopt ==16
2139 
2140 
2141      !  !===================================================================================================
2142      !  !                                      fixed reduced d and ebar calculation (mixed BC)
2143      !  !===================================================================================================
2144      if (dtset%berryopt == 17) then
2145 
2146         if (scfcv_step > 1) then
2147            !      ! update efield taking damping into account reduced red_dfield
2148            !      red_efield1 and red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel1999]]
2149            !      red_efield1 for fixed ebar, red_efield2 for fixed d calculation
2150 
2151            !      save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar
2152            dtefield%efield2(:)=dtset%efield(:)
2153 
2154            !      write(*,'(a,3i4)') "jfielddir=", (dtset%jfielddir(ii),ii=1,3)
2155 
2156            do ii=1,3
2157               if (dtset%jfielddir(ii) ==2 ) then    ! direction under fixed d
2158                  dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) !  update ebar which is not fixed
2159                  dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii))
2160                  red_efield2(ii)=dtset%ddamp*(dtset%red_dfield(ii) - red_ptot(ii)) +  &
2161                       &           (1.0d0-dtset%ddamp)*red_efield2_old(ii)         ! d(ii) is fixed, update e(ii)  may need ddamping here
2162 
2163                  !          write(message,'(a,a,i5,a,i5)')   ch10, 'direction  ', ii,'   for fixed d, value is (2)  ', dtset%jfielddir(ii)
2164                  !          call wrtout(ab_out,message,'COLL')
2165                  !          call wrtout(std_out,message,'COLL')
2166 
2167               else if (dtset%jfielddir(ii) ==1 ) then   ! direction under fixed ebar
2168                  red_efield2(ii)= (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) !  update e which is not fixed
2169                  dtset%red_dfield(ii)=red_ptot(ii) +  (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii))  ! update d
2170 
2171                  !          write(message,'(a,a,i5,a,i5)')   ch10, 'direction  ', ii,'   for fixed ebar, value is (1)  ', dtset%jfielddir(ii)
2172                  !          call wrtout(ab_out,message,'COLL')
2173                  !          call wrtout(std_out,message,'COLL')
2174 
2175               end if
2176            end do
2177 
2178            do ii=1,3
2179               red_efield1(ii)  =(ucvol/(4*pi))*dot_product(dtset%red_efieldbar(:),gmet(:,ii))
2180            end do
2181 
2182 
2183            dtset%red_efield(:)=(red_efield1(:) + red_efield2(:))/2.0d0 ! average reduced efield,
2184            !      one is from fixed ebar part,
2185            !      the other is from fixed d part.
2186            !      This may need to be optimized !!
2187 
2188            write(message,'(a,a,a,a,3(es16.9,2x),a)')   ch10, 'Reduced efield from fixed ebar:', ch10, &
2189                 &       '       e:  ', (red_efield1(ii),ii=1,3), ch10
2190 
2191            !      call wrtout(ab_out,message,'COLL')
2192            call wrtout(std_out,message,'COLL')
2193 
2194            write(message,'(a,a,a,a,3(es16.9,2x),a)')   ch10, 'Reduced efield from fixed d:', ch10, &
2195                 &       '       e:  ', (red_efield2(ii),ii=1,3), ch10
2196 
2197            !      call wrtout(ab_out,message,'COLL')
2198            call wrtout(std_out,message,'COLL')
2199 
2200            write(message,'(a,a,a,a,3(es16.9,2x),a)')   ch10, 'Average reduced efield:', ch10, &
2201                 &       '       e:  ', (dtset%red_efield(ii),ii=1,3), ch10
2202 
2203            !      call wrtout(ab_out,message,'COLL')
2204            call wrtout(std_out,message,'COLL')
2205 
2206            !      to calculate unreduced E
2207            do ii=1,3
2208               efield_test_cart(ii)  = (4*pi/ucvol)* dot_product(dtset%red_efield(:),rprimd(:,ii))
2209            end do
2210 
2211            !      ! test whether change in efield in any direction exceed maxestep, if so, set the
2212            !      ! change to maxestep instead   ! need optimized !
2213            do ii = 1,3
2214              if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then
2215 
2216                write(std_out,'(a,a,i5)') "JH - ","  E-field component:",ii
2217                write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), &
2218 &               ",    E(n-1)=",efield_old_cart(ii), ",    E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), &
2219 &               ",    maxestep=",dtset%maxestep
2220 
2221                if (efield_test_cart(ii) > efield_old_cart(ii)) then
2222                  efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep)
2223                else
2224                  efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep)
2225                end if
2226              end if
2227            end do
2228 
2229            dtset%efield(:) = efield_test_cart(:)
2230 
2231            !      !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar
2232            write(message,'(a,a)')   ch10, 'scfcv: Constant reduced ebar and d-field  - updating E-field:'
2233            call wrtout(std_out,message,'COLL')
2234            call prtefield(dtset,dtefield,std_out,rprimd)
2235            if(dtset%prtvol>=10)then
2236               call wrtout(ab_out,message,'COLL')
2237               call prtefield(dtset,dtefield,ab_out,rprimd)
2238            end if
2239 
2240 
2241            !      ! need to update dtset%efield_dot(:) with new value
2242            !      ! This needs to be deleted  when efield_dot is deleted
2243            dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
2244            dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
2245            dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
2246 
2247         else
2248 
2249            write(message,'(a,a)')   ch10, 'scfcv: Constant reduced ebar and d-field  - Pre E-field:'
2250            call wrtout(std_out,message,'COLL')
2251            call prtefield(dtset,dtefield,std_out,rprimd)
2252            if(dtset%prtvol>=10)then
2253               call wrtout(ab_out,message,'COLL')
2254               call prtefield(dtset,dtefield,ab_out,rprimd)
2255            end if
2256 
2257         end if  ! scfcv_step > 1
2258 
2259         efield_old_cart(:)=dtset%efield(:)
2260         red_efield2_old(:)=red_efield2(:)
2261 
2262      end if  ! berryopt ==17
2263 
2264   end if ! end efield .and. scfcv_level 1 tasks
2265 
2266   !deallocate cprj
2267   if ( efield .and. psps%usepaw == 1) then
2268      call pawcprj_free(cprj)
2269   end if
2270   ABI_DATATYPE_DEALLOCATE(cprj)
2271 
2272 end subroutine update_e_field_vars